hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
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);