lcsim/src/org/lcsim/contrib/tracking
diff -u -r1.1 -r1.2
--- StandaloneAxialBarrelTrack1.java 9 Jul 2007 06:09:40 -0000 1.1
+++ StandaloneAxialBarrelTrack1.java 9 Jul 2007 19:58:54 -0000 1.2
@@ -8,15 +8,18 @@
* Open. You can then make changes to the template in the Source Editor.
*/
-package org.lcsim.contrib.tracking;
+package org.lcsim.recon.tracking.standalone;
-import hep.physics.matrix.SymmetricMatrix;
import org.lcsim.event.Track;
import org.lcsim.event.TrackerHit;
import org.lcsim.event.SimTrackerHit;
import org.lcsim.fit.circle.CircleFit;
-
+import org.lcsim.event.MCParticle;
import java.util.*;
+import java.io.*;
+import org.lcsim.util.aida.AIDA;
+import hep.aida.*;
+import hep.physics.matrix.SymmetricMatrix;
/**
* Implements the Track interface for tracks found using only outer tracker information by TrackFinder
@@ -44,14 +47,15 @@
double _module_length = 100.0;
List<Track> _tracks = new ArrayList<Track>();
boolean _isreferencepoint_pca = false;
-
+ MCParticle _majority_particle;
+ double _purity = 0.0;
// Static
- static int _type = 2; // ???
+ static int _type = 2; // ???
// Constructors
//=============
-
+
/**
* Create a default instance of StandaloneAxialBarrelTrack
*/
@@ -61,12 +65,12 @@
/**
* Create a StandaloneAxialBarrelTrack using information from TrackFinder
- *
+ *
* @param b_field Magnetic field: required to generate momentum from curvature
* @param module_length Length of readout modules: required to estimate tan(lambda)
* @param fit The CircleFit for this track
* @param hits The list of SimTrackerHits used for this track
- *
+ *
* @see org.lcsim.fit.circle.CircleFit
*/
public StandaloneAxialBarrelTrack1(double b_field, double module_length, CircleFit fit, List<SimTrackerHit> hits)
@@ -82,20 +86,50 @@
_reference_point[1] = fit.yref();
_reference_point[2] = 0.0;
- _sim_tracker_hits = hits;
+ final List<SimTrackerHit> _sim_tracker_hits = hits;
_track_parameters[0] = fit.dca();
_track_parameters[1] = fit.phi()-Math.PI; // who ordered this?
_track_parameters[2] = 10.0*fit.curvature(); // convert 1/mm to 1/cm
_track_parameters[3] = 0.0;
_track_parameters[4] = getTanLambda();
-
+
_isreferencepoint_pca = true;
_tracks.add(this);
}
-
+
+ public StandaloneAxialBarrelTrack1(double b_field, double module_length, CircleFit fit, List<SimTrackerHit> hits, MCParticle majority_particle, double purity)
+ {
+ _b_field = b_field;
+ _module_length = module_length;
+
+ _fit_success = true;
+ _chi2 = fit.chisq();
+ _ndof = hits.size();
+
+ _reference_point[0] = fit.xref();
+ _reference_point[1] = fit.yref();
+ _reference_point[2] = 0.0;
+
+ _sim_tracker_hits = hits;
+
+ _track_parameters[0] = fit.dca();
+ _track_parameters[1] = fit.phi()-Math.PI; // who ordered this?
+ _track_parameters[2] = fit.curvature(); // convert 1/mm to 1/cm
+ _track_parameters[3] = 0.0;
+ _track_parameters[4] = getTanLambda();
+
+ _isreferencepoint_pca = true;
+
+ _majority_particle = majority_particle;
+ _purity = purity;
+
+ _tracks.add(this);
+
+ }
+
// Methods needed to further initialize a StandaloneAxialBarrelTrack
//=======================================================
@@ -149,21 +183,21 @@
}
/**
- * The y coordinate of the reference point for the track parameters
- */
+ * The y coordinate of the reference point for the track parameters
+ */
public double getReferencePointY()
{
return _reference_point[1];
}
-
- /**
- * The z coordinate of the reference point for the track parameters
- */
+
+ /**
+ * The z coordinate of the reference point for the track parameters
+ */
public double getReferencePointZ()
{
return _reference_point[2];
}
-
+
/**
* Is this reference point the point of closest approach used to define the track parameters?
*/
@@ -187,27 +221,27 @@
}
/**
- * The x-component of the momentum of the particle corresponding to this track: depends upon
+ * The x-component of the momentum of the particle corresponding to this track: depends upon
* b_field being appropriately when constructing the StandaloneAxialBarrelTrack
*/
public double getPX()
{
return getMomentum()[0];
}
-
+
/**
- * The y-component of the momentum of the particle corresponding to this track: depends upon
+ * The y-component of the momentum of the particle corresponding to this track: depends upon
* b_field being appropriately when constructing the StandaloneAxialBarrelTrack
- */
+ */
public double getPY()
{
return getMomentum()[1];
}
/**
- * The z-component of the momentum of the particle corresponding to this track: depends upon
+ * The z-component of the momentum of the particle corresponding to this track: depends upon
* b_field being appropriately when constructing the StandaloneAxialBarrelTrack
- */
+ */
public double getPZ()
{
return getMomentum()[2];
@@ -245,7 +279,7 @@
return _track_parameters;
}
- /**
+ /**
* Retreive the error matrix of the track fit
* <p>
* <b>N.B. Currently returns an empty matrix</b>
@@ -255,7 +289,7 @@
return _error_matrix;
}
- /**
+ /**
* Retreive the chisquared of the track fit
*/
public double getChi2()
@@ -281,7 +315,7 @@
return _dedx;
}
- /**
+ /**
* Retrieve the dEdx error of the track: zero unless set explicitly
*<p>
*<b>N.B. Currently zero unless set explicitly</b>
@@ -296,10 +330,14 @@
*/
public double getRadiusOfInnermostHit()
{
- double min_radius = 1000000.0; // 1km
- for (SimTrackerHit hit : _sim_tracker_hits){
- min_radius = Math.min(Math.sqrt(hit.getPoint()[0]*hit.getPoint()[0] +
- hit.getPoint()[1]*hit.getPoint()[1]), min_radius);
+ double min_radius = 1000000.0; // 1km
+ for (SimTrackerHit hit : _sim_tracker_hits)
+ {
+ if(hit != null)
+ {
+ min_radius = Math.min(Math.sqrt(hit.getPoint()[0]*hit.getPoint()[0] +
+ hit.getPoint()[1]*hit.getPoint()[1]), min_radius);
+ }
}
return min_radius;
}
@@ -315,7 +353,7 @@
}
/**
- * Get tracks associated with this track.
+ * Get tracks associated with this track.
* <p>
* <b>Currently returns only this unless set explicitly</b>
*/
@@ -355,7 +393,8 @@
/**
* Calculate tan(lambda) using outermost hit radius and module length
*/
- private double getTanLambda() {
+ private double getTanLambda()
+ {
SimTrackerHit hit = getOutermostHit();
double x_outer = hit.getPoint()[0];
double y_outer = hit.getPoint()[1];
@@ -380,17 +419,119 @@
{
double max_radius = 0.0;
SimTrackerHit outermost_hit = null;
- for (SimTrackerHit hit : _sim_tracker_hits){
+ for (SimTrackerHit hit : _sim_tracker_hits)
+ {
if (hit == null) continue;
double radius = Math.sqrt(hit.getPoint()[0]*hit.getPoint()[0] +
- hit.getPoint()[1]*hit.getPoint()[1]);
+ hit.getPoint()[1]*hit.getPoint()[1]);
if ( radius > max_radius)
{
max_radius = radius;
- outermost_hit = hit;
+ outermost_hit = hit;
}
}
return outermost_hit;
}
+ /**
+ *Get majority MCParticle associated with track
+ */
+ public MCParticle getMCParticle()
+ {
+ return _majority_particle;
+ }
+
+ /**
+ *Return percentage of hits from majority particle
+ */
+ public double purity()
+ {
+ return _purity;
+ }
+ public double chiSquared()
+ {
+ double resolution = .007; //arclength resolution in mm
+
+ //Get track parameters
+ double dca = _track_parameters[0];
+ double phi = _track_parameters[1];
+
+ // double rtrack = -1.0/_track_parameters[2]; //radius of track in mm, I hope // **GOLD**
+
+ double rtrack = -1.0/_track_parameters[2];
+
+ // int sign = (int)Math.signum(rtrack); // **GOLD**
+
+ int sign = (int)Math.signum(rtrack);
+
+ double charge = getCharge();
+ double tp2 = _track_parameters[2];
+
+ //System.out.println( "Sign is: " +sign+ " and rtrack is: " +rtrack);
+ //System.out.println( "TP2: " +tp2+ "Charge: " +charge);
+
+ /* GOLD
+ //Calculate center of circle defined by track, in polar coordinates
+ double rcenter = Math.abs(rtrack) - (sign*dca);
+ double phicenter = phi-sign*0.5*Math.PI;
+ END GOLD! */
+
+ double rcenter = Math.abs(rtrack) - (sign*dca);
+ double phicenter = sign*phi-0.5*Math.PI;
+
+ /*
+ if (sign > 0){
+ // phicenter = phi - (sign*0.5*Math.PI); // **GOLD**
+
+ //phicenter = sign*phi - (0.5*Math.PI);
+
+ }
+ if (sign < 0){
+ // phicenter = -1*phi+(sign*0.5*Math.PI); //**GOLD**
+
+ //phicenter = sign*phi-(0.5*Math.PI);
+
+ }
+ */
+
+ double wssr = 0.0; //sum of squares of residuals
+
+ //Calculate residuals
+ for(SimTrackerHit hit : _sim_tracker_hits)
+ {
+ if(hit != null)
+ {
+ //Find polar coordinates of hit
+ double rhit = Math.sqrt(hit.getPoint()[0]*hit.getPoint()[0] + hit.getPoint()[1]*hit.getPoint()[1]);
+ double phihit = Math.atan2(hit.getPoint()[1],hit.getPoint()[0]);
+ if(phihit > Math.PI) phihit -= 2*Math.PI;
+ if(phihit < -1.0 * Math.PI) phihit += 2*Math.PI;
+
+ //Find phi coordinate of the point on the track with same radius as hit.
+ //This equation comes from solving for the phi coordinate of the intersection of two circles.
+ //There are, of course, two solutions. I hope I did the right thing to pick which one is good.
+ double deltaPhi = Math.acos((rhit*rhit + rcenter*rcenter - rtrack*rtrack)/(2.0 * rhit * rcenter));
+ double phitrack = 0.0;
+ double phitrack1 = sign*phicenter + deltaPhi;
+ double phitrack2 = sign*phicenter - deltaPhi;
+ if(Math.abs(phihit - phitrack1) < Math.abs(phihit - phitrack2)) phitrack = phitrack1;
+ else phitrack = phitrack2;
+ if(Math.abs(phitrack) > Math.PI)
+ {
+ if(phitrack>0) phitrack -= 2*Math.PI;
+ if(phitrack<0) phitrack += 2*Math.PI;
+ }
+ double dPhi = phihit - phitrack;
+ //Find weighted squared residual
+ double residual = (dPhi * dPhi) * (rhit * rhit)/(resolution * resolution);
+ wssr += residual;
+ //System.out.println("R Hit is: " +rhit+ " dPhi is: " +dPhi);
+ //double logchisquared = Math.log10(wssr);
+ //aida.cloud2D("Log ChiSquared vs. Phi").fill(logchisquared, phi);
+ }
+
+
+ }
+ return wssr;
+ }
}
\ No newline at end of file