 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/recon/tracking on MAIN | |||
| gbl/GBLFileIO.java | +2 | -2 | 1.7 -> 1.8 |
| /GBLOutput.java | +136 | -296 | 1.10 -> 1.11 |
| /GBLOutputDriver.java | +19 | -32 | 1.5 -> 1.6 |
| TrackUtils.java | +28 | -78 | 1.23 -> 1.24 |
| +185 | -408 | ||
Cleaning up.
diff -u -r1.7 -r1.8 --- GBLFileIO.java 12 Sep 2013 17:57:38 -0000 1.7 +++ GBLFileIO.java 15 Sep 2013 19:24:49 -0000 1.8 @@ -34,8 +34,8 @@
openFile();
}
- public void printEventNr(int evtnr) {
- addLine(String.format("New Event %d", evtnr));
+ public void printEventInfo(int evtnr, double Bz) {
+ addLine(String.format("New Event %d %.12f", evtnr, Bz));
}
protected void addLine(String line) {
diff -u -r1.10 -r1.11 --- GBLOutput.java 12 Sep 2013 17:57:38 -0000 1.10 +++ GBLOutput.java 15 Sep 2013 19:24:49 -0000 1.11 @@ -9,6 +9,7 @@
import hep.physics.matrix.MatrixOp; import hep.physics.matrix.SymmetricMatrix; import hep.physics.vec.*;
+
import java.io.FileWriter; import java.io.PrintWriter; import java.util.ArrayList;
@@ -16,6 +17,7 @@
import java.util.HashMap; import java.util.List; import java.util.Map;
+
import org.lcsim.constants.Constants; import org.lcsim.event.*; import org.lcsim.fit.helicaltrack.*;
@@ -39,11 +41,13 @@
private int _debug;
private GBLFileIO file;
private Hep3Vector _B;
- private TrackerHitUtils _trackerHitUtils = new TrackerHitUtils();
+ private TrackerHitUtils _trackerHitUtils = new TrackerHitUtils();
private MaterialSupervisor _materialmanager;
private MultipleScattering _scattering;
private HPSTransformations _hpstrans = new HPSTransformations();
private double _beamEnergy = 2.2; //GeV
+ private boolean AprimeEvent = false; + private boolean hasXPlanes = false;
@@ -66,8 +70,8 @@
public void buildModel(Detector detector) {
_materialmanager.buildModel(detector);
}
- void printNewEvent(int eventNumber) {
- file.printEventNr(eventNumber);
+ void printNewEvent(int eventNumber,double Bz) {
+ file.printEventInfo(eventNumber,Bz);
}
void printTrackID(int iTrack) {
file.printTrackID(iTrack);
@@ -75,145 +79,106 @@
void close() {
file.closeFile();
}
+ void setAPrimeEventFlag(boolean flag) {
+ this.AprimeEvent = flag;
+ }
+ void setXPlaneFlag(boolean flag) {
+ this.hasXPlanes = flag;
+ }
+ public Hep3Vector get_B() {
+ return _B;
+ }
+ public void set_B(Hep3Vector _B) {
+ this._B = _B;
+ }
+
void printGBL(Track trk, List<MCParticle> mcParticles, List<SimTrackerHit> simTrackerHits) {
SeedTrack st = (SeedTrack)trk;
SeedCandidate seed = st.getSeedCandidate();
- HelicalTrackFit htf = seed.getHelix(); - double[] htf_params = new double[5]; - System.arraycopy(htf.parameters(), 0, htf_params, 0, 5); - htf_params[0] *= -1; - HelicalTrackFit htf_dca_flip = new HelicalTrackFit(htf_params, htf.covariance(), htf.chisq(), htf.ndf(), htf.PathMap(), htf.ScatterMap()); - //htf = htf_dca_flip;
+ HelicalTrackFit htf = seed.getHelix(); + _scattering.setDebug(this._debug>0?true:false); + ScatterPoints scatters = _scattering.FindHPSScatterPoints(htf); + +
+ // Hits on track
List<HelicalTrackHit> hits = seed.getHits();
// Find the truth particle of the track
MCParticle mcp = getMatchedTruthParticle(trk);
+
+ if(mcp==null) {
+ System.out.printf("%s: no truth particle found!\n",this.getClass().getSimpleName());
+ System.exit(1);
+ }
- // Get track paramteres from MC particle
+ if(AprimeEvent ) {
+ checkAprimeTruth(mcp,mcParticles);
+ }
+
+
+ // Get track parameters from MC particle
HelicalTrackFit htfTruth = getHTF(mcp);
// Use the truth helix as the initial track for GBL?
//htf = htfTruth;
+
-
- if (this._debug>1) {
- System.out.printf("%s: find scatters\n",this.getClass().getSimpleName());
- }
- //find the scattering angle from material manager
- _scattering.setDebug(this._debug>0?true:false);
- ScatterPoints scatters = _scattering.FindHPSScatterPoints(htf);
-
-
- // Convert from perigee (L3 parameters) to curvilinear frame
+ // Get perigee parameters to curvilinear frame
PerigeeParams perPar = new PerigeeParams(htf);
+ PerigeeParams perParTruth = new PerigeeParams(htfTruth); + file.printPerTrackParam(perPar); + file.printPerTrackParamTruth(perParTruth); + + // Get curvilinear parameters
ClParams clPar = new ClParams(htf);
- //BasicMatrix jacPerToCl = getJacPerToCl(htf); - //BasicMatrix clPar = (BasicMatrix) MatrixOp.mult(jacPerToCl, MatrixOp.transposed(perPar)); - //clPar = (BasicMatrix) MatrixOp.transposed(clPar); - - // print chi2 of fit - file.printChi2(htf.chisq(),htf.ndf());
+ ClParams clParTruth = new ClParams(htfTruth); + file.printClTrackParam(clPar); + file.printClTrackParamTruth(clParTruth);
- // print track parameters for easier debugging later - file.printOldPerTrackParam(htf);
if(_debug>0) {
System.out.printf("%s\n",file.getClTrackParamStr(clPar));
System.out.printf("%s\n",file.getPerTrackParamStr(perPar));
}
- file.printPerTrackParam(perPar);
- file.printClTrackParam(clPar);
-
-
-
-
- // Use the sim tracker hits to find the truth particle
- // this would give me the particle kinematics responsible for the track
-
- List<MCParticle> mcp_pair = getAprimeDecayProducts(mcParticles);
- if(_debug>0) {
- double invMassTruth = Math.sqrt( Math.pow(mcp_pair.get(0).getEnergy()+mcp_pair.get(1).getEnergy(),2) - VecOp.add(mcp_pair.get(0).getMomentum(), mcp_pair.get(1).getMomentum()).magnitudeSquared() );
- double invMassTruthTrks = getInvMassTracks(getHTF(mcp_pair.get(0)),getHTF(mcp_pair.get(1)));
- System.out.printf("%s: invM = %f\n",this.getClass().getSimpleName(),invMassTruth);
- System.out.printf("%s: invMTracks = %f\n",this.getClass().getSimpleName(),invMassTruthTrks);
- }
+
-
+ // find the projection from the I,J,K to U,V,T curvilinear coordinates + Hep3Matrix perToClPrj = getPerToClPrj(htf); + file.printPerToClPrj(perToClPrj); + + // print chi2 of fit + file.printChi2(htf.chisq(),htf.ndf());
- // cross-check
- if(!mcp_pair.contains(mcp)) {
- boolean hasBeamElectronParent = false;
- for(MCParticle parent : mcp.getParents()) {
- if(parent.getGeneratorStatus()!=MCParticle.FINAL_STATE && parent.getPDGID()==11 && parent.getMomentum().y()==0.0 && Math.abs(parent.getMomentum().magnitude() - _beamEnergy)<0.01) {
- hasBeamElectronParent = true;
+ // build map of layer to SimTrackerHits that belongs to the MC particle
+ Map<Integer, SimTrackerHit> simHitsLayerMap = new HashMap<Integer, SimTrackerHit >();
+ for (SimTrackerHit sh : simTrackerHits) {
+ if(sh.getMCParticle()==mcp) {
+ int layer = sh.getIdentifierFieldValue("layer");
+ if(!simHitsLayerMap.containsKey(layer) || (sh.getPathLength() < simHitsLayerMap.get(layer).getPathLength()) ) {
+ simHitsLayerMap.put(layer, sh);
}
}
- if(!hasBeamElectronParent) {
- System.out.printf("%s: the matched MC particle is not an A' daughter and not a the recoil electrons!?\n",this.getClass().getSimpleName());
- System.out.printf("%s: %s %d p %s org %s\n",this.getClass().getSimpleName(),mcp.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",mcp.getPDGID(),mcp.getMomentum().toString(),mcp.getOrigin().toString());
- printMCParticles(mcParticles);
- System.exit(1);
- } else {
- if(_debug>0) System.out.printf("%s: the matched MC particle is the recoil electron\n",this.getClass().getSimpleName());
- }
- }
-
- List<SimTrackerHit> simHits = new ArrayList<SimTrackerHit>();
- for (SimTrackerHit hit : simTrackerHits) {
- if(hit.getMCParticle()==mcp) simHits.add(hit);
- }
-
- if(_debug>0) {
- System.out.printf("%s: %d sim tracker hits to MC particle with p = %f (trk p = %f q=%f):\n",this.getClass().getSimpleName(),simHits.size(),mcp.getMomentum().magnitude(),htf.p(this._B.magnitude()),Math.signum(htf.R()));
- for(SimTrackerHit sh : simHits) System.out.printf("%s: sim hit at layer %d and pos %s \n",this.getClass().getSimpleName(),sh.getIdentifierFieldValue("layer"),sh.getPositionVec().toString());
- }
-
- if(mcp==null) {
- System.out.printf("%s: no truth particle found!\n",this.getClass().getSimpleName());
- System.exit(1);
- }
-
- // map layer and sim hits
- Map<Integer, SimTrackerHit> simHitsLayerMap = new HashMap<Integer, SimTrackerHit >();
- for(SimTrackerHit sh : simHits) {
- int layer = sh.getIdentifierFieldValue("layer");
- if(!simHitsLayerMap.containsKey(layer) || (sh.getPathLength() < simHitsLayerMap.get(layer).getPathLength()) ) {
- simHitsLayerMap.put(layer, sh);
- }
}
- - PerigeeParams perParTruth = new PerigeeParams(htfTruth); - file.printPerTrackParamTruth(perParTruth); - ClParams clParTruth = new ClParams(htfTruth); - file.printClTrackParamTruth(clParTruth); - - - // find the projection from the I,J,K to U,V,T coordinates - Hep3Matrix perToClPrj = this.getPerToClPrj(htf); - - file.printPerToClPrj(perToClPrj); - - // covariance matrix from the helix fit
+ // covariance matrix from the fit
file.printPerTrackCov(htf);
- // calculate truth chi2 - double chi2 = truthTrackFitChi2(perPar,perParTruth,htf.covariance());
+ // dummy cov matrix for CL parameters + BasicMatrix clCov = new BasicMatrix(5,5); + initUnit(clCov); + clCov = (BasicMatrix) MatrixOp.mult(0.1*0.1,clCov); + file.printCLTrackCov(clCov); +
if(_debug>0) {
System.out.printf("%s: perPar covariance matrix\n%s\n",this.getClass().getSimpleName(),htf.covariance().toString());
- System.out.printf("%s: truth perPar chi2 %f\n",this.getClass().getSimpleName(),chi2);
+ double chi2_truth = truthTrackFitChi2(perPar,perParTruth,htf.covariance());
+ System.out.printf("%s: truth perPar chi2 %f\n",this.getClass().getSimpleName(),chi2_truth);
}
- - - BasicMatrix clCov = new BasicMatrix(5,5); - this.initUnit(clCov); - clCov = (BasicMatrix) MatrixOp.mult(0.1*0.1,clCov); - file.printCLTrackCov(clCov);
if(_debug>0) System.out.printf("%d hits\n",hits.size());
@@ -232,8 +197,7 @@
file.printStrip(istrip,strip.layer());
//Center of the sensor
- Hep3Vector origin = strip.origin(); -
+ Hep3Vector origin = strip.origin();
file.printOrigin(origin);
// associated 3D position of the crossing of this and it's stereo partner sensor
@@ -241,109 +205,39 @@
//Find intercept point with sensor in tracking frame
Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(htf, strip, Math.abs(_B.z()));
- Hep3Vector trkposXPlane = TrackUtils.getHelixXPlaneIntercept(htf, strip.w(), origin);
- Hep3Vector trkposXPlaneIter = TrackUtils.getHelixPlanePositionIter(htf, strip.origin(), strip.w(), 1.0e-8);
- Hep3Vector trkposDiff = VecOp.sub(trkposXPlaneIter, trkpos);
- if(trkposDiff.magnitude() > 1.0e-7) {
- System.out.printf("WARNING trkposDiff mag = %.10f [%.10f %.10f %.10f]\n",trkposDiff.magnitude(),trkposDiff.x(),trkposDiff.y(),trkposDiff.z());
- System.exit(1);
- }
-
+ Hep3Vector trkposTruth = TrackUtils.getHelixPlaneIntercept(htfTruth, strip, Math.abs(_B.z())); + file.printStripTrackPos(trkpos); +
if(_debug>0) {
System.out.printf("trkpos at intercept [%.10f %.10f %.10f]\n",trkpos.x(),trkpos.y(),trkpos.z());
- //System.out.printf("trkposXPlane at intercept %s\n",trkposXPlane.toString());
- System.out.printf("trkposXPlaneIter at intercept [%.10f %.10f %.10f]\n",trkposXPlaneIter.x(),trkposXPlaneIter.y(),trkposXPlaneIter.z());
- // is it on the plane?
- //double trkpos_d = VecOp.dot(strip.w(),VecOp.unit(VecOp.sub(trkpos, strip.origin())));
- //double trkposXPlaneIter_d = VecOp.dot(strip.w(),VecOp.unit(VecOp.sub(trkposXPlaneIter, strip.origin())));
- //System.out.printf("trkpos_d %f with diff %s trkpos %s org %s \n",trkpos_d,VecOp.unit(VecOp.sub(trkpos, strip.origin())), trkpos.toString(),strip.origin().toString());
- //System.out.printf("trkposXPlane_d %f with diff %s trkpos %s org %s \n",trkposXPlaneIter_d,VecOp.unit(VecOp.sub(trkposXPlaneIter, strip.origin())), trkposXPlaneIter.toString(),strip.origin().toString());
+ System.out.printf("trkposTruth at intercept %s\n",trkposTruth.toString());
}
- //trkpos = trkposXPlane;
- //trkpos = trkposXPlaneDCAFlip;
-
-// if(strip.layer()==1) trkpos = new BasicHep3Vector(88.8078567539 ,0.213035827282, 1.5292560945);
-// if(strip.layer()==2) trkpos = new BasicHep3Vector(96.2602646374, 0.378914765887 , 1.686096828);
-// if(strip.layer()==3) trkpos = new BasicHep3Vector(188.751065257, 3.59937204783, 3.63336621837);
-// if(strip.layer()==4) trkpos = new BasicHep3Vector(196.197763346 ,3.95227914386, 3.7902238256);
-// if(strip.layer()==5) trkpos = new BasicHep3Vector(288.617639967, 9.49588626783 ,5.7383167412);
-// if(strip.layer()==6) trkpos = new BasicHep3Vector(296.05862198 ,10.0360590169 ,5.89529022351 );
-// if(strip.layer()==7) trkpos = new BasicHep3Vector(488.64433778 ,28.9100728781, 9.96718984117);
-// if(strip.layer()==8) trkpos = new BasicHep3Vector(496.071890128, 29.8276792157, 10.1246568417);
-// if(strip.layer()==9) trkpos = new BasicHep3Vector(687.835519286, 58.4555308451, 14.2045621589);
-// if(strip.layer()==10) trkpos = new BasicHep3Vector(695.251412131, 59.7550500287, 14.3629733178);
-// if(strip.layer()==11) trkpos = new BasicHep3Vector(886.710664575, 98.3529112595, 18.472815436);
-// if(strip.layer()==12) trkpos = new BasicHep3Vector(894.114639235, 100.042816527 ,18.6326045072 );
-//
-
-
- if(_debug>0) {
- //calculate track position for a given number of x-plane positions
- double[] x_tmp = {88.8078567539,96.2602646374,188.751065257,196.197763346,288.617639967,296.05862198,488.64433778,496.071890128,687.835519286,695.251412131,886.710664575,894.114639235};
- for(int i=0;i<x_tmp.length;++i) {
- double s_tmp = HelixUtils.PathToXPlane(htf, x_tmp[i], 0., 0).get(0);
- Hep3Vector trkpos_x_tmp = HelixUtils.PointOnHelix(htf, s_tmp);
-
- //Calculate path length by hand for x plane
- double dca = htf.dca();
- double xc = (htf.R() - dca) * Math.sin(htf.phi0());
- double sinPhi = (xc - x_tmp[i])/htf.R();
- double phi_at_x = Math.asin(sinPhi);
- //double phi_at_x = Math.asin( 1/htf.R()* ( -1*dca*Math.sin(htf.phi0()) + htf.R()*Math.sin(htf.phi0()) - trkpos.x() ) );
- double dphi_at_x = phi_at_x - htf.phi0();
- if(dphi_at_x > Math.PI) dphi_at_x -= 2.0*Math.PI;
- if(dphi_at_x < -Math.PI) dphi_at_x += 2.0*Math.PI;
- double s_at_x = -1.0*dphi_at_x*htf.R();
- double y_at_x = dca*Math.cos(htf.phi0()) - htf.R()*Math.cos(htf.phi0()) + htf.R()*Math.cos(phi_at_x);
- double z_at_x = htf.z0() + s_at_x*htf.slope();
- Hep3Vector trkposmanual_x_tmp = new BasicHep3Vector(x_tmp[i],y_at_x,z_at_x);
- //System.out.printf("s = %f trkpos %s (initial manual xplane)\n", s_at_x , trkposmanual_at_x.toString());
- //System.out.printf("params %f %f %f %f %f\n", dca,htf.z0(),htf.phi0(),htf.slope(),htf.curvature());
- System.out.printf("x %f gives trkpos at %s for s %f s3D %f\n", x_tmp[i],trkpos_x_tmp.toString(),s_tmp,s_tmp/Math.cos(Math.atan(htf.slope())));
- System.out.printf("x %f gives trkpos at %s for s %f s3D %f (manual)\n", x_tmp[i],trkposmanual_x_tmp.toString(),s_at_x,s_at_x/Math.cos(Math.atan(htf.slope())));
-
- }
+ // cross-check intercept point
+ if(hasXPlanes ) {
+ Hep3Vector trkposXPlaneIter = TrackUtils.getHelixPlanePositionIter(htf, strip.origin(), strip.w(), 1.0e-8);
+ Hep3Vector trkposDiff = VecOp.sub(trkposXPlaneIter, trkpos);
+ if(trkposDiff.magnitude() > 1.0e-7) {
+ System.out.printf("WARNING trkposDiff mag = %.10f [%.10f %.10f %.10f]\n",trkposDiff.magnitude(),trkposDiff.x(),trkposDiff.y(),trkposDiff.z());
+ System.exit(1);
+ }
+ if(_debug>0) System.out.printf("trkposXPlaneIter at intercept [%.10f %.10f %.10f]\n",trkposXPlaneIter.x(),trkposXPlaneIter.y(),trkposXPlaneIter.z());
}
-
-
- Hep3Vector trkposTruth = TrackUtils.getHelixPlaneIntercept(htfTruth, strip, Math.abs(_B.z()));
-
- if(_debug>0) {
- System.out.printf("trkposTruth at intercept %s\n",trkposTruth.toString());
- }
-
// Find the sim tracker hit for this layer
SimTrackerHit simHit = simHitsLayerMap.get(strip.layer());
-
if(simHit==null) {
System.out.printf("%s: no sim hit for strip hit at layer %d\n",this.getClass().getSimpleName(),strip.layer());
System.out.printf("%s: it as %d mc particles associated with it:\n",this.getClass().getSimpleName(),hit.getMCParticles().size());
- for (MCParticle particle : hit.getMCParticles()) {
- System.out.printf("%s: %d p %s \n",this.getClass().getSimpleName(),particle.getPDGID(),particle.getMomentum().toString());
- }
- System.out.printf("%s: these are sim hits in the layer map:\n",this.getClass().getSimpleName());
- for(Map.Entry<Integer,SimTrackerHit> entry : simHitsLayerMap.entrySet()) {
- SimTrackerHit simhit = entry.getValue();
- System.out.printf("%s layer %d sim hit at %s with MC particle pdgid %d with p %s \n",this.getClass().getSimpleName(),entry.getKey(),simhit.getPositionVec().toString(),simhit.getMCParticle().getPDGID(),simhit.getMCParticle().getMomentum().toString());
- }
+ for (MCParticle particle : hit.getMCParticles()) System.out.printf("%s: %d p %s \n",this.getClass().getSimpleName(),particle.getPDGID(),particle.getMomentum().toString());
System.out.printf("%s: these are sim hits in the event:\n",this.getClass().getSimpleName());
- for (SimTrackerHit simhit : simTrackerHits) {
- System.out.printf("%s sim hit at %s with MC particle pdgid %d with p %s \n",this.getClass().getSimpleName(),simhit.getPositionVec().toString(),simhit.getMCParticle().getPDGID(),simhit.getMCParticle().getMomentum().toString());
- }
+ for (SimTrackerHit simhit : simTrackerHits) System.out.printf("%s sim hit at %s with MC particle pdgid %d with p %s \n",this.getClass().getSimpleName(),simhit.getPositionVec().toString(),simhit.getMCParticle().getPDGID(),simhit.getMCParticle().getMomentum().toString());
System.out.printf("%s: these are all the MC particles in the event:\n",this.getClass().getSimpleName());
- for(MCParticle loop_p : mcParticles) {
- System.out.printf("%s: %d p=%s parents ",this.getClass().getSimpleName(),loop_p.getPDGID(),loop_p.getMomentum().toString());
- for(MCParticle parent : loop_p.getParents()) {
- System.out.printf(" [%d p=%s]",parent.getPDGID(),parent.getMomentum().toString());
- }
- System.out.printf("\n");
- }
- //System.exit(1);
+ System.exit(1);
}
+
if(_debug>0) {
double s_truthSimHit = HelixUtils.PathToXPlane(htfTruth, simHit.getPositionVec().z(), 0, 0).get(0);
Hep3Vector trkposTruthSimHit = HelixUtils.PointOnHelix(htfTruth, s_truthSimHit);
@@ -354,35 +248,10 @@
//path length to intercept
double s = HelixUtils.PathToXPlane(htf,trkpos.x(),0,0).get(0);
double s3D = s / Math.cos(Math.atan(htf.slope()));
-
-
-
- if(_debug>0) {
-
-
- double s_XPlane = HelixUtils.PathToXPlane(htf,trkposXPlane.x(),0,0).get(0);
- double s_truth = HelixUtils.PathToXPlane(htfTruth, trkposTruth.x(), 0, 0).get(0);
- double s3D_truth = s_truth / Math.cos(Math.atan(htfTruth.slope()));
- double s3D_XPlane = s_XPlane / Math.cos(Math.atan(htf.slope()));
- System.out.printf("s = %f (s_3D=%f) trkpos %s (initial)\n", s, s3D , trkpos.toString());
- System.out.printf("s = %f (s_3D=%f) trkpos %s (truth) \n", s_truth, s3D_truth, trkposTruth.toString());
- System.out.printf("s = %f (s_3D=%f) trkpos %s (XPlane) \n", s_XPlane, s3D_XPlane, trkposXPlane.toString());
- double s_test = 88.8207081021 * Math.cos(Math.atan(htf.slope()));
- Hep3Vector trkpos_test = HelixUtils.PointOnHelix(htf, s_test);
- System.out.printf("s = %f trkpos %s (initial)\n", s_test, trkpos_test.toString());
- s_test = 88.8103791522 * Math.cos(Math.atan(htfTruth.slope()));
- trkpos_test = HelixUtils.PointOnHelix(htfTruth, s_test);
- System.out.printf("s = %f trkpos %s (truth)\n", s_test, trkpos_test.toString());
- }
-
- //print the path length
file.printStripPathLen(s);
file.printStripPathLen3D(s3D);
- //file.addLine(String.format("trkpos at s=%f is %s",88.7866518242,HelixUtils.PointOnHelix(htf, 88.7866518242).toString()));
-
//print stereo angle in YZ plane
- if(_debug>0) System.out.printf("u = %s v = %s w = %s \n",strip.u().toString(),strip.v().toString(), strip.w().toString());
file.printMeasDir(strip.u());
file.printNonMeasDir(strip.v());
file.printNormalDir(strip.w());
@@ -393,82 +262,56 @@
double lambda = Math.atan(htf.slope());
file.printStripTrackDir(Math.sin(phi),Math.sin(lambda));
file.printStripTrackDirFull(tDir);
- file.printStripTrackPos(trkpos);
- //Print residual in measurment system
+ //Print residual in measurement system
// start by find the distance vector between the center and the track position
Hep3Vector vdiffTrk = VecOp.sub(trkpos, origin);
Hep3Vector vdiffTrkTruth = VecOp.sub(trkposTruth, origin);
- Hep3Vector vdiffTrkIter = VecOp.sub(trkposXPlaneIter, origin);
// then find the rotation from tracking to measurement frame
Hep3Matrix trkToStripRot = _trackerHitUtils.getTrackToStripRotation(strip);
+
// then rotate that vector into the measurement frame to get the predicted measurement position
Hep3Vector trkpos_meas = VecOp.mult(trkToStripRot, vdiffTrk);
Hep3Vector trkposTruth_meas = VecOp.mult(trkToStripRot, vdiffTrkTruth);
+ // hit measurement and uncertainty in measurement frame + Hep3Vector m_meas = new BasicHep3Vector(strip.umeas(),0.,0.); + Hep3Vector res_err_meas = new BasicHep3Vector(strip.du(),(strip.vmax() - strip.vmin()) / Math.sqrt(12),10.0/Math.sqrt(12));
- - // hit uncertainty in measurement frame - double umeas = strip.umeas(); - double uError = strip.du(); - double vmeas = 0; - double vError = (strip.vmax() - strip.vmin()) / Math.sqrt(12); - double wmeas = 0; - double wError = 10.0/Math.sqrt(12); //0.001; - - // measurement in measurement frame - Hep3Vector m_meas = new BasicHep3Vector(umeas,vmeas,wmeas); - - file.printStripMeas(strip.umeas());
+ file.printStripMeas(m_meas.x());
// residual in measurement frame
Hep3Vector res_meas = VecOp.sub(m_meas, trkpos_meas);
Hep3Vector resTruth_meas = VecOp.sub(m_meas, trkposTruth_meas);
- double pred_meas_iter = VecOp.dot(strip.u(),vdiffTrkIter);
- double res_meas_iter = umeas - pred_meas_iter;
-
- if(_debug>0) {
- System.out.printf("layer %d uRes %.10f uRes_alt %.10f\n",strip.layer(),res_meas.x(),res_meas_iter);
- }
+ file.printStripMeasRes(res_meas.x(),res_err_meas.x());
+ file.printStripMeasResTruth(resTruth_meas.x(),res_err_meas.x());
+
+ if(_debug>0) System.out.printf("layer %d uRes %.10f\n",strip.layer(),res_meas.x());
- // residual uncertainty in measurement frame - Hep3Vector res_err_meas = new BasicHep3Vector(uError,vError,wError); - Hep3Vector resTruth_err_meas = new BasicHep3Vector(uError,vError,wError); - Hep3Vector resSimHit_err_meas = new BasicHep3Vector(uError,vError,wError);
+ // sim hit residual
- // print measurement to file - file.printStripMeasRes(res_meas.x(),res_err_meas.x()); - file.printStripMeasResTruth(resTruth_meas.x(),resTruth_err_meas.x());
if(simHit!=null) {
- Hep3Vector simHitPos = null;simHitPos = this._hpstrans.transformVectorToTracking(simHit.getPositionVec());
- if(_debug>0) {
- System.out.printf("simHitPos %s\n",simHitPos.toString());
- }
+ Hep3Vector simHitPos = this._hpstrans.transformVectorToTracking(simHit.getPositionVec());
+ if(_debug>0) System.out.printf("simHitPos %s\n",simHitPos.toString());
Hep3Vector vdiffSimHit = VecOp.sub(simHitPos, trkpos);
Hep3Vector simHitPos_meas = VecOp.mult(trkToStripRot, vdiffSimHit);
- Hep3Vector resSimHit_meas = simHitPos_meas; //VecOp.sub(m_meas, simHitPos_meas); - file.printStripMeasResSimHit(resSimHit_meas.x(),resSimHit_err_meas.x());
+ file.printStripMeasResSimHit(simHitPos_meas.x(),res_err_meas.x());
} else {
file.printStripMeasResSimHit(-999999.9,-999999.9);
}
- - // print scattering angle
+ // find scattering angle
ScatterPoint scatter = scatters.getScatterPoint(((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement());
double scatAngle;
if(scatter != null) {
-
scatAngle = scatter.getScatterAngle().Angle();
-
}
else {
-
System.out.printf("%s: WARNING cannot find scatter for detector %s with strip cluster at %s\n",this.getClass(),((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement().getName(),strip.origin().toString());
-
//can be edge case where helix is outside, but close to sensor, so use hack with the sensor origin closest to hit -> FIX THIS!
DetectorPlane closest = null;
double dx = 999999.9;
@@ -515,18 +358,36 @@
}
- List<MCParticle> getTruthParticle(int pdgId, List<MCParticle> mcParticles) {
-
- if(mcParticles==null) return null;
+ private void checkAprimeTruth(MCParticle mcp, List<MCParticle> mcParticles) {
+ List<MCParticle> mcp_pair = getAprimeDecayProducts(mcParticles);
- List<MCParticle> fsParticles = new ArrayList<MCParticle>();
- for (MCParticle mcparticle : mcParticles) {
- if (mcparticle.getGeneratorStatus() == MCParticle.FINAL_STATE && pdgId==mcparticle.getPDGID()) {
- fsParticles.add(mcparticle);
- }
- }
- return fsParticles;
- }
+ if(_debug>0) {
+ double invMassTruth = Math.sqrt( Math.pow(mcp_pair.get(0).getEnergy()+mcp_pair.get(1).getEnergy(),2) - VecOp.add(mcp_pair.get(0).getMomentum(), mcp_pair.get(1).getMomentum()).magnitudeSquared() );
+ double invMassTruthTrks = getInvMassTracks(getHTF(mcp_pair.get(0)),getHTF(mcp_pair.get(1)));
+ System.out.printf("%s: invM = %f\n",this.getClass().getSimpleName(),invMassTruth);
+ System.out.printf("%s: invMTracks = %f\n",this.getClass().getSimpleName(),invMassTruthTrks);
+ }
+
+ // cross-check
+ if(!mcp_pair.contains(mcp)) {
+ boolean hasBeamElectronParent = false;
+ for(MCParticle parent : mcp.getParents()) {
+ if(parent.getGeneratorStatus()!=MCParticle.FINAL_STATE && parent.getPDGID()==11 && parent.getMomentum().y()==0.0 && Math.abs(parent.getMomentum().magnitude() - _beamEnergy)<0.01) {
+ hasBeamElectronParent = true;
+ }
+ }
+ if(!hasBeamElectronParent) {
+ System.out.printf("%s: the matched MC particle is not an A' daughter and not a the recoil electrons!?\n",this.getClass().getSimpleName());
+ System.out.printf("%s: %s %d p %s org %s\n",this.getClass().getSimpleName(),mcp.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",mcp.getPDGID(),mcp.getMomentum().toString(),mcp.getOrigin().toString());
+ printMCParticles(mcParticles);
+ System.exit(1);
+ } else {
+ if(_debug>0) System.out.printf("%s: the matched MC particle is the recoil electron\n",this.getClass().getSimpleName());
+ }
+ }
+
+ }
+
MCParticle getMatchedTruthParticle(Track track) {
boolean debug = false;
@@ -576,8 +437,6 @@
private BasicMatrix getPerParVector(HelicalTrackFit htf) {
BasicMatrix perPar = new BasicMatrix(1,5);
- Hep3Vector T = HelixUtils.Direction(htf, 0.); - Hep3Vector p = VecOp.mult(htf.p(Math.abs(_B.z())), T);
double kappa = -1.0*Math.signum(htf.R())*Constants.fieldConversion*this._B.z()/htf.pT(Math.abs(_B.z()));
double theta = Math.PI/2.0 - Math.atan(htf.slope());
perPar.setElement(0,0,kappa);
@@ -589,21 +448,7 @@
}
-
- private Hep3Vector getMomentum(Track trk) {
- double pvec[] = trk.getTrackStates().get(0).getMomentum();
- //Hep3Vector p = VecOp.mult(1.0e3,new BasicHep3Vector(pvec[0],pvec[1],pvec[2]));
- Hep3Vector p = new BasicHep3Vector(pvec[0],pvec[1],pvec[2]);
-
-
- return p;
- }
-
- private double getLambda(Track trk) {
- return Math.atan(trk.getTrackStates().get(0).getTanLambda());
- }
-
private BasicMatrix getJacPerToCl(HelicalTrackFit htf) {
@@ -627,7 +472,6 @@
double alpha = VecOp.cross(H,T).magnitude();
Hep3Vector N = VecOp.mult(1./alpha,VecOp.cross(H, T));
Hep3Vector K = Z;
- double gamma = VecOp.dot(H, T);
double Q = -Bnorm.magnitude()*q/p.magnitude();
double kappa = -1.0*q*Bnorm.z()/pT;
@@ -748,12 +592,6 @@
pt.setElement(0, 4, Math.tan(Math.PI/2.0-perParTruth.getTheta()));
Matrix error_matrix = MatrixOp.inverse(covariance);
BasicMatrix res = (BasicMatrix) MatrixOp.sub(p, pt);
-// for(int i=0;i<res.getNRows();++i) {
-// for(int j=0;j<res.getNColumns();++j) {
-// double v = res.e(i, j)*res.e(i, j);
-// res.setElement(i, j, v);
-// }
-// }
BasicMatrix chi2 = (BasicMatrix) MatrixOp.mult(res,MatrixOp.mult(error_matrix, MatrixOp.transposed(res)));
if(chi2.getNColumns()!=1 ||chi2.getNRows()!=1) {
throw new RuntimeException("matrix dim is screwed up!");
@@ -885,7 +723,9 @@
}
- public class ClParams {
+
+
+ public class ClParams {
private BasicMatrix _params = new BasicMatrix(1,5);
private ClParams(HelicalTrackFit htf) {
@@ -900,7 +740,7 @@
//System.out.printf("%s: vecCl=%s\n",this.getClass().getSimpleName(),vecCl.toString());
double xT = vecCl.x();
double yT = vecCl.y();
- double zT = vecCl.z();
+ //double zT = vecCl.z();
Hep3Vector T = HelixUtils.Direction(htf, 0.);
Hep3Vector p = VecOp.mult(htf.p(Math.abs(_B.z())), T);
diff -u -r1.5 -r1.6 --- GBLOutputDriver.java 12 Sep 2013 17:57:38 -0000 1.5 +++ GBLOutputDriver.java 15 Sep 2013 19:24:49 -0000 1.6 @@ -2,10 +2,13 @@
import hep.physics.vec.BasicHep3Vector; import hep.physics.vec.Hep3Vector;
+
import java.io.IOException;
+import java.util.ArrayList;
import java.util.List; import java.util.logging.Level; import java.util.logging.Logger;
+
import org.lcsim.event.EventHeader; import org.lcsim.event.MCParticle; import org.lcsim.event.SimTrackerHit;
@@ -13,7 +16,6 @@
import org.lcsim.geometry.Detector; import org.lcsim.hps.recon.tracking.EventQuality; import org.lcsim.hps.recon.tracking.TrackUtils;
-import org.lcsim.hps.users.phansson.TrigRateDriver;
import org.lcsim.util.Driver; import org.lcsim.util.aida.AIDA;
@@ -34,7 +36,6 @@
private String outputPlotFileName="";
private boolean hideFrame = false;
private int _debug = 0;
- private String simTrackerHitCollectionName = "TrackerHits";
private String MCParticleCollectionName = "MCParticle";
private int iTrack = 0;
private int iEvent = 0;
@@ -91,50 +92,36 @@
}
List<SimTrackerHit> simTrackerHits = event.getSimTrackerHits("TrackerHits");
- truthRes.processSim(mcParticles, simTrackerHits);
+ //truthRes.processSim(mcParticles, simTrackerHits);
- //gbl.printNewEvent(event.getEventNumber()); - gbl.printNewEvent(iEvent); - - iTrack = 0;
+ + + List<Track> selected_tracks = new ArrayList<Track>();
for (Track trk : tracklist) {
-
- //if(trk.getCharge()>0) continue;
- //if(trk.getTrackStates().get(0).getMomentum()[0]>0.8) continue;
-
- totalTracks++;
-
- TrackUtils _trackUtils = new TrackUtils();
-
- if(!TrackUtils.isGoodTrack(trk, tracklist, EventQuality.Quality.MEDIUM)) {
- if(_debug>0) {
- System.out.printf("%s: Track failed selection\n", this.getClass().getSimpleName());
- }
- continue;
+ totalTracks++;
+ if(TrackUtils.isGoodTrack(trk, tracklist, EventQuality.Quality.MEDIUM)) {
+ if(_debug>0) System.out.printf("%s: Track failed selection\n", this.getClass().getSimpleName());
+ selected_tracks.add(trk);
}
+ }
- if(_debug>0) {
- System.out.printf("%s: Print GBL output for this track\n", this.getClass().getSimpleName());
- }
-
- totalTracksProcessed++;
-
- gbl.printTrackID(iTrack);
+ //gbl.printNewEvent(event.getEventNumber()); + gbl.printNewEvent(iEvent,gbl.get_B().z());
-
+ iTrack = 0;
+ for (Track trk : selected_tracks) {
+ if(_debug>0) System.out.printf("%s: Print GBL output for this track\n", this.getClass().getSimpleName());
+ gbl.printTrackID(iTrack);
gbl.printGBL(trk,mcParticles,simTrackerHits);
- -
+ totalTracksProcessed++;
++iTrack;
}
++iEvent;
- -
}
@Override
diff -u -r1.23 -r1.24 --- TrackUtils.java 12 Sep 2013 17:57:54 -0000 1.23 +++ TrackUtils.java 15 Sep 2013 19:24:49 -0000 1.24 @@ -30,7 +30,7 @@
/** * * @author Omar Moreno <[log in to unmask]>
- * @version $Id: TrackUtils.java,v 1.23 2013/09/12 17:57:54 phansson Exp $
+ * @version $Id: TrackUtils.java,v 1.24 2013/09/15 19:24:49 phansson Exp $
* TODO: Switch to JLab coordinates */
@@ -289,6 +289,29 @@
}
+ /*
+ * Calculates the point on the helix in the x-y plane at the intercept with plane
+ * The normal of the plane is in the same x-y plane as the circle.
+ *
+ * @param helix
+ * @param vector normal to plane
+ * @param origin of plane
+ * @return point in the x-y plane of the intercept
+ *
+ */
+ public Hep3Vector getHelixXPlaneIntercept(HelicalTrackFit helix, Hep3Vector w, Hep3Vector origin) {
+ throw new RuntimeException("this function is not working properly; don't use it");
+
+ // FInd the intercept point x_int,y_int, between the circle and sensor, which becomes a line in the x-y plane in this case.
+ // y_int = k*x_int + m
+ // R^2 = (y_int-y_c)^2 + (x_int-x_c)^2
+ // solve for x_int
+
+ }
+
+
+
+
/**
* @param helix input helix object
* @param origin of the plane to intercept
@@ -353,69 +376,7 @@
return pos;
}
- /*
- * Calculates the point on the helix in the x-y plane at the intercept with plane
- * The normal of the plane is in the same x-y plane as the circle.
- *
- * @param helix
- * @param vector normal to plane
- * @param origin of plane
- * @return point in the x-y plane of the intercept
- *
- */
- public static Hep3Vector getHelixXPlaneIntercept(HelicalTrackFit helix, Hep3Vector w, Hep3Vector origin) {
- // FInd the intercept point x_int,y_int, between the circle and sensor, which becomes a line in the x-y plane in this case.
- // y_int = k*x_int + m
- // R^2 = (y_int-y_c)^2 + (x_int-x_c)^2
- // solve for x_int
- double R = helix.R();
- double d0 = helix.dca();
- double sinPhi0 = Math.sin(helix.phi0());
- double cosPhi0 = Math.cos(helix.phi0());
- //double x_0 = -d0*sinPhi0; // POCA
- //double y_0 = d0*cosPhi0; // POCA
- double x_c = (R-d0)*sinPhi0; // center of circle
- double y_c = -1*(R-d0)*cosPhi0; // center of circle
- double k = -1*w.x()/w.y(); // dy/dx slope of line, note change of sign
- double m = origin.y() - k * origin.x(); // intercept
-
- double k2 = k*k;
- double R2 = R*R;
- double m2 = m*m;
- double A = k2*(m+y_c);
- double B = Math.sqrt( k2*( -m2 + (1+k2)*R2 +x_c*x_c + k2*x_c*x_c - 2*m*y_c - y_c*y_c ) );
- double C = k+k*k*k;
-
- double xint1 = (A - B)/C;
- double xint2 = (A + B)/C;
-
- xint1 *= -1.;
- xint2 *= -1.;
-
- //if(Double.isInfinite(xint1) || Double.isNaN(xint1)) {
- //System.out.printf("xint1 = %f xint2 = %f\n",xint1,xint2);
- double xint;
- if(xint1>0 && xint2>0) {
- xint = xint1<xint2 ? xint1 : xint2; // chose smaller one
- } else if(xint1<0 && xint2<0) {
- System.out.printf("problemxint1 = %f xint2 = %f\n",xint1,xint2);
- xint = xint1;
- } else {
- xint = xint1>0 ? xint1 : xint2;
- }
- double s = HelixUtils.PathToXPlane(helix, xint, 0, 0).get(0);
- Hep3Vector pointOnHelix = HelixUtils.PointOnHelix(helix, s);
- //double[] point = {xint,k*xint+m};
- return pointOnHelix;
- }
-
-
-
-
-
- /**
- *
- */
+
public double getPhi(Hep3Vector position){
// Check if a track has been set
@@ -554,13 +515,9 @@
double s = HelixUtils.PathToXPlane(track, hth.x(), 0, 0).get(0);
//System.out.printf("x %f s %f smap %f\n",hth.x(),s,s_wrong);
if(Double.isNaN(s)) {
- double x0 = track.x0(); - double y0 = track.y0(); -
double xc=track.xc();
double yc=track.yc();
double RC = track.R();
- double y=yc+Math.signum(RC)*Math.sqrt(RC*RC-Math.pow(hth.x()-xc,2));
System.out.printf("calculateTrackHitResidual: s is NaN. p=%.3f RC=%.3f, x=%.3f, xc=%.3f\n",track.p(-0.491),RC,hth.x(),xc);
}
@@ -624,16 +581,10 @@
Hep3Vector u = strip.u();
- Hep3Vector v = strip.v(); - Hep3Vector w = strip.w(); - Hep3Vector eta = VecOp.cross(u,v);
Hep3Vector corigin = strip.origin();
- String side = SvtUtils.getInstance().isTopLayer((SiSensor)((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom"; - - double bfield = bFieldInZ; //Math.abs(this._bfield.z());
//Find interception with plane that the strips belongs to
- Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(_trk, strip, bfield);
+ Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(_trk, strip, bFieldInZ);
if(debug) {
System.out.printf("calculateLocalTrackHitResiduals: found interception point at %s \n",trkpos.toString());
@@ -643,8 +594,8 @@
if(Double.isNaN(trkpos.x()) || Double.isNaN(trkpos.y()) || Double.isNaN(trkpos.z())) {
System.out.printf("calculateLocalTrackHitResiduals: failed to get interception point (%s) \n",trkpos.toString());
System.out.printf("calculateLocalTrackHitResiduals: track params\n%s\n",_trk.toString());
- System.out.printf("calculateLocalTrackHitResiduals: track pT=%.3f chi2=[%.3f][%.3f] \n",_trk.pT(bfield),_trk.chisq()[0],_trk.chisq()[1]);
- trkpos = TrackUtils.getHelixPlaneIntercept(_trk, strip, bfield);
+ System.out.printf("calculateLocalTrackHitResiduals: track pT=%.3f chi2=[%.3f][%.3f] \n",_trk.pT(bFieldInZ),_trk.chisq()[0],_trk.chisq()[1]);
+ trkpos = TrackUtils.getHelixPlaneIntercept(_trk, strip, bFieldInZ);
System.exit(1);
}
@@ -731,7 +682,6 @@
HelicalTrackHit hth = (HelicalTrackHit) hit;
for(Track track : othertracks) {
List<TrackerHit> hitsOnTrack = track.getTrackerHits();
- boolean shared = false;
for(TrackerHit loop_hit : hitsOnTrack) {
HelicalTrackHit loop_hth = (HelicalTrackHit) loop_hit;
if(hth.equals(loop_hth)) {
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1