Author: [log in to unmask]
Date: Fri Aug 28 16:36:48 2015
New Revision: 3452
Log:
clean up before I change stuff
Modified:
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java Fri Aug 28 16:36:48 2015
@@ -1,4 +1,5 @@
package org.hps.recon.tracking.gbl;
+
import hep.physics.matrix.BasicMatrix;
import hep.physics.matrix.Matrix;
import hep.physics.matrix.MatrixOp;
@@ -41,37 +42,36 @@
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
-
/**
* Calculate the input needed for Millepede minimization.
+ *
* @author Per Hansson Adrian <[log in to unmask]>
*
*/
public class GBLOutput {
-
+
private int _debug = 0;
private GBLFileIO textFile = null;
private Hep3Vector _B;
- private final TrackerHitUtils _trackerHitUtils = new TrackerHitUtils();
+ private final TrackerHitUtils _trackerHitUtils = new TrackerHitUtils();
private final MaterialSupervisor _materialmanager;
private final MultipleScattering _scattering;
private final double _beamEnergy = 1.1; //GeV
- private boolean AprimeEvent = false; // do extra checks
- private boolean hasXPlanes = false;
-
-
-
-
+ private boolean AprimeEvent = false; // do extra checks
+ private boolean hasXPlanes = false;
+
/**
* Constructor
- * @param outputFileName is the filename given to the text-based output file. If empty no output file is written
+ *
+ * @param outputFileName is the filename given to the text-based output
+ * file. If empty no output file is written
* @param bfield magnetic field in Tesla
*/
- GBLOutput(String outputFileName,Hep3Vector bfield) {
- //System.out.printf("name \"%s\" \n", outputFileName);
- if(!outputFileName.equalsIgnoreCase("")) {
- textFile = new GBLFileIO(outputFileName);
- }
+ GBLOutput(String outputFileName, Hep3Vector bfield) {
+ //System.out.printf("name \"%s\" \n", outputFileName);
+ if (!outputFileName.equalsIgnoreCase("")) {
+ textFile = new GBLFileIO(outputFileName);
+ }
_materialmanager = new MaterialSupervisor();
_scattering = new MultipleScattering(_materialmanager);
_B = CoordinateTransformations.transformVectorToTracking(bfield);
@@ -80,661 +80,623 @@
public void setDebug(int debug) {
_debug = debug;
- _scattering.setDebug((_debug>0));
- }
+ _scattering.setDebug((_debug > 0));
+ }
+
public void buildModel(Detector detector) {
_materialmanager.buildModel(detector);
}
- void printNewEvent(int eventNumber,double Bz) {
- if(textFile != null) {
- textFile.printEventInfo(eventNumber,Bz);
- }
- }
+
+ void printNewEvent(int eventNumber, double Bz) {
+ if (textFile != null) {
+ textFile.printEventInfo(eventNumber, Bz);
+ }
+ }
+
void printTrackID(int iTrack) {
- if(textFile != null) {
- textFile.printTrackID(iTrack);
- }
- }
+ if (textFile != null) {
+ textFile.printTrackID(iTrack);
+ }
+ }
+
void close() {
- if(textFile != null) {
- textFile.closeFile();
- }
- }
+ if (textFile != null) {
+ textFile.closeFile();
+ }
+ }
+
void setAPrimeEventFlag(boolean flag) {
- this.AprimeEvent = flag;
- }
+ this.AprimeEvent = flag;
+ }
+
void setXPlaneFlag(boolean flag) {
- this.hasXPlanes = flag;
- }
+ this.hasXPlanes = flag;
+ }
+
public Hep3Vector get_B() {
- return _B;
- }
- public void set_B(Hep3Vector _B) {
- this._B = _B;
- }
-
-
-
+ return _B;
+ }
+
+ public void set_B(Hep3Vector _B) {
+ this._B = _B;
+ }
+
void printGBL(Track trk, List<SiTrackerHitStrip1D> stripHits, GBLTrackData gtd, List<GBLStripClusterData> stripClusterDataList, List<MCParticle> mcParticles, List<SimTrackerHit> simTrackerHits, boolean isMC) {
- SeedTrack st = (SeedTrack)trk;
+ SeedTrack st = (SeedTrack) trk;
SeedCandidate seed = st.getSeedCandidate();
- HelicalTrackFit htf = seed.getHelix();
+ HelicalTrackFit htf = seed.getHelix();
// Find scatter points along the path
ScatterPoints scatters = _scattering.FindHPSScatterPoints(htf);
-
+
// Hits on track
List<HelicalTrackHit> hits = seed.getHits();
// Find the truth particle of the track
MCParticle mcp = null;
-
- if(isMC) {
- mcp = getMatchedTruthParticle(trk);
-
- if(mcp==null) {
- System.out.printf("%s: WARNING!! no truth particle found in event!\n",this.getClass().getSimpleName());
- this.printMCParticles(mcParticles);
- //System.exit(1);
- } else {
- if(_debug>0) System.out.printf("%s: truth particle (pdgif %d ) found in event!\n",this.getClass().getSimpleName(),mcp.getPDGID());
- }
-
- if(AprimeEvent ) {
- checkAprimeTruth(mcp,mcParticles);
- }
- }
-
+
+ if (isMC) {
+ mcp = getMatchedTruthParticle(trk);
+
+ if (mcp == null) {
+ System.out.printf("%s: WARNING!! no truth particle found in event!\n", this.getClass().getSimpleName());
+ this.printMCParticles(mcParticles);
+ //System.exit(1);
+ } else {
+ if (_debug > 0) {
+ System.out.printf("%s: truth particle (pdgif %d ) found in event!\n", this.getClass().getSimpleName(), mcp.getPDGID());
+ }
+ }
+
+ if (AprimeEvent) {
+ checkAprimeTruth(mcp, mcParticles);
+ }
+ }
+
// Get track parameters from MC particle
- HelicalTrackFit htfTruth = (isMC&&mcp!=null) ? TrackUtils.getHTF(mcp,-1.0*this._B.z()) : null;
-
+ HelicalTrackFit htfTruth = (isMC && mcp != null) ? TrackUtils.getHTF(mcp, -1.0 * this._B.z()) : null;
+
// Use the truth helix as the initial track for GBL?
//htf = htfTruth;
-
-
// Get perigee parameters to curvilinear frame
PerigeeParams perPar = new PerigeeParams(htf);
PerigeeParams perParTruth = new PerigeeParams(htfTruth);
- if(textFile != null) {
- textFile.printPerTrackParam(perPar);
- textFile.printPerTrackParamTruth(perParTruth);
- }
-
+ if (textFile != null) {
+ textFile.printPerTrackParam(perPar);
+ textFile.printPerTrackParamTruth(perParTruth);
+ }
+
//GBLDATA
gtd.setPerigeeTrackParameters(perPar);
// Get curvilinear parameters
ClParams clPar = new ClParams(htf);
ClParams clParTruth = new ClParams(htfTruth);
- if(textFile != null) {
- textFile.printClTrackParam(clPar);
- textFile.printClTrackParamTruth(clParTruth);
-
- if(_debug>0) {
- System.out.printf("%s\n",textFile.getClTrackParamStr(clPar));
- System.out.printf("%s\n",textFile.getPerTrackParamStr(perPar));
- }
- }
-
-
+ if (textFile != null) {
+ textFile.printClTrackParam(clPar);
+ textFile.printClTrackParamTruth(clParTruth);
+
+ if (_debug > 0) {
+ System.out.printf("%s\n", textFile.getClTrackParamStr(clPar));
+ System.out.printf("%s\n", textFile.getPerTrackParamStr(perPar));
+ }
+ }
+
// find the projection from the I,J,K to U,V,T curvilinear coordinates
Hep3Matrix perToClPrj = getPerToClPrj(htf);
-
- if(textFile != null) {
- textFile.printPerToClPrj(perToClPrj);
- }
-
+
+ if (textFile != null) {
+ textFile.printPerToClPrj(perToClPrj);
+ }
+
//GBLDATA
- for(int row=0; row<perToClPrj.getNRows();++row) {
- for(int col=0; col<perToClPrj.getNColumns();++col) {
- gtd.setPrjPerToCl(row, col, perToClPrj.e(row, col));
- }
- }
-
-
-
+ for (int row = 0; row < perToClPrj.getNRows(); ++row) {
+ for (int col = 0; col < perToClPrj.getNColumns(); ++col) {
+ gtd.setPrjPerToCl(row, col, perToClPrj.e(row, col));
+ }
+ }
+
//GBLDATA
- for(int row=0; row<perToClPrj.getNRows();++row) {
- for(int col=0; col<perToClPrj.getNColumns();++col) {
- gtd.setPrjPerToCl(row, col, perToClPrj.e(row, col));
- }
- }
-
-
+ for (int row = 0; row < perToClPrj.getNRows(); ++row) {
+ for (int col = 0; col < perToClPrj.getNColumns(); ++col) {
+ gtd.setPrjPerToCl(row, col, perToClPrj.e(row, col));
+ }
+ }
+
// print chi2 of fit
- if(textFile != null) {
- textFile.printChi2(htf.chisq(),htf.ndf());
- }
-
+ if (textFile != null) {
+ textFile.printChi2(htf.chisq(), htf.ndf());
+ }
+
// build map of layer to SimTrackerHits that belongs to the MC particle
- Map<Integer, SimTrackerHit> simHitsLayerMap = new HashMap<Integer, SimTrackerHit >();
+ 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()) ) {
+ if (sh.getMCParticle() == mcp) {
+ int layer = sh.getIdentifierFieldValue("layer");
+ if (!simHitsLayerMap.containsKey(layer) || (sh.getPathLength() < simHitsLayerMap.get(layer).getPathLength())) {
simHitsLayerMap.put(layer, sh);
}
}
}
-
-
+
// covariance matrix from the fit
- if(textFile != null) {
- textFile.printPerTrackCov(htf);
- }
-
+ if (textFile != null) {
+ textFile.printPerTrackCov(htf);
+ }
+
// dummy cov matrix for CL parameters
- BasicMatrix clCov = new BasicMatrix(5,5);
+ BasicMatrix clCov = new BasicMatrix(5, 5);
initUnit(clCov);
- clCov = (BasicMatrix) MatrixOp.mult(0.1*0.1,clCov);
-
- if(textFile != null) {
- textFile.printCLTrackCov(clCov);
- }
-
- if(_debug>0) {
- System.out.printf("%s: perPar covariance matrix\n%s\n",this.getClass().getSimpleName(),htf.covariance().toString());
- double chi2_truth = truthTrackFitChi2(perPar,perParTruth,htf.covariance());
- System.out.printf("%s: truth perPar chi2 %f\n",this.getClass().getSimpleName(),chi2_truth);
- }
-
- if(_debug>0) System.out.printf("%d hits\n",hits.size());
-
+ clCov = (BasicMatrix) MatrixOp.mult(0.1 * 0.1, clCov);
+
+ if (textFile != null) {
+ textFile.printCLTrackCov(clCov);
+ }
+
+ if (_debug > 0) {
+ System.out.printf("%s: perPar covariance matrix\n%s\n", this.getClass().getSimpleName(), htf.covariance().toString());
+ double chi2_truth = truthTrackFitChi2(perPar, perParTruth, htf.covariance());
+ System.out.printf("%s: truth perPar chi2 %f\n", this.getClass().getSimpleName(), chi2_truth);
+ }
+
+ if (_debug > 0) {
+ System.out.printf("%d hits\n", hits.size());
+ }
int istrip = 0;
- for(int ihit=0;ihit!=hits.size();++ihit) {
-
+ for (int ihit = 0; ihit != hits.size(); ++ihit) {
+
HelicalTrackHit hit = hits.get(ihit);
HelicalTrackCross htc = (HelicalTrackCross) hit;
List<HelicalTrackStrip> strips = htc.getStrips();
-
- for(HelicalTrackStrip stripOld : strips) {
+
+ for (HelicalTrackStrip stripOld : strips) {
HelicalTrackStripGbl strip = new HelicalTrackStripGbl(stripOld, true);
-
+
// find Millepede layer definition from DetectorElement
IDetectorElement de = ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement();
HpsSiSensor sensor;
- if( de instanceof HpsSiSensor) {
+ if (de instanceof HpsSiSensor) {
sensor = (HpsSiSensor) de;
} else {
throw new ClassCastException("Detector element " + de.getName() + " couldn't be casted to " + HpsSiSensor.class.getName());
}
int millepedeId = sensor.getMillepedeId();
-
- if(_debug>0) System.out.printf("%s: layer %d millepede %d (DE=\"%s\", origin %s) \n",this.getClass().getSimpleName(),strip.layer(), millepedeId, sensor.getName(), strip.origin().toString());
-
- if(textFile != null) {
- textFile.printStrip(istrip,millepedeId,de.getName());
- }
-
+
+ if (_debug > 0) {
+ System.out.printf("%s: layer %d millepede %d (DE=\"%s\", origin %s) \n", this.getClass().getSimpleName(), strip.layer(), millepedeId, sensor.getName(), strip.origin().toString());
+ }
+
+ if (textFile != null) {
+ textFile.printStrip(istrip, millepedeId, de.getName());
+ }
+
//GBLDATA
GBLStripClusterData stripData = new GBLStripClusterData(millepedeId);
//Add to output list
stripClusterDataList.add(stripData);
-
-
-
+
//Center of the sensor
- Hep3Vector origin = strip.origin();
-
- if(textFile != null) {
- textFile.printOrigin(origin);
- }
-
+ Hep3Vector origin = strip.origin();
+
+ if (textFile != null) {
+ textFile.printOrigin(origin);
+ }
+
// associated 3D position of the crossing of this and it's stereo partner sensor
- if(textFile != null) {
- textFile.printHitPos3D(hit.getCorrectedPosition());
- }
-
+ if (textFile != null) {
+ textFile.printHitPos3D(hit.getCorrectedPosition());
+ }
+
//Find intercept point with sensor in tracking frame
Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(htf, strip, Math.abs(_B.z()));
- Hep3Vector trkposTruth = htfTruth!=null ? TrackUtils.getHelixPlaneIntercept(htfTruth, strip, Math.abs(_B.z())) : new BasicHep3Vector(-999999.9,-999999.9,-999999.9);
- if(textFile != null) {
- textFile.printStripTrackPos(trkpos);
- }
- if(_debug>0) {
- System.out.printf("trkpos at intercept [%.10f %.10f %.10f]\n",trkpos.x(),trkpos.y(),trkpos.z());
- System.out.printf("trkposTruth at intercept %s\n",trkposTruth!=null?trkposTruth.toString():"no truth track");
- }
-
+ Hep3Vector trkposTruth = htfTruth != null ? TrackUtils.getHelixPlaneIntercept(htfTruth, strip, Math.abs(_B.z())) : new BasicHep3Vector(-999999.9, -999999.9, -999999.9);
+ if (textFile != null) {
+ textFile.printStripTrackPos(trkpos);
+ }
+ if (_debug > 0) {
+ System.out.printf("trkpos at intercept [%.10f %.10f %.10f]\n", trkpos.x(), trkpos.y(), trkpos.z());
+ System.out.printf("trkposTruth at intercept %s\n", trkposTruth != null ? trkposTruth.toString() : "no truth track");
+ }
+
// 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());
- }
-
-
+ 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());
+ }
+ }
+
// Find the sim tracker hit for this layer
SimTrackerHit simHit = simHitsLayerMap.get(strip.layer());
- if( isMC ) {
- 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 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());
- //System.out.printf("%s: these are all the MC particles in the event:\n",this.getClass().getSimpleName());
- //System.exit(1);
- }
-
- if(_debug>0) {
- if(htfTruth!=null && simHit!=null) {
- double s_truthSimHit = HelixUtils.PathToXPlane(htfTruth, simHit.getPositionVec().z(), 0, 0).get(0);
- Hep3Vector trkposTruthSimHit = HelixUtils.PointOnHelix(htfTruth, s_truthSimHit);
- Hep3Vector resTruthSimHit = VecOp.sub(CoordinateTransformations.transformVectorToTracking(simHit.getPositionVec()),trkposTruthSimHit);
- System.out.printf("TruthSimHit residual %s for layer %d\n",resTruthSimHit.toString(),strip.layer());
- }
- }
- }
-
+ if (isMC) {
+ 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 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());
+ }
+ //System.out.printf("%s: these are all the MC particles in the event:\n",this.getClass().getSimpleName());
+ //System.exit(1);
+ }
+
+ if (_debug > 0) {
+ if (htfTruth != null && simHit != null) {
+ double s_truthSimHit = HelixUtils.PathToXPlane(htfTruth, simHit.getPositionVec().z(), 0, 0).get(0);
+ Hep3Vector trkposTruthSimHit = HelixUtils.PointOnHelix(htfTruth, s_truthSimHit);
+ Hep3Vector resTruthSimHit = VecOp.sub(CoordinateTransformations.transformVectorToTracking(simHit.getPositionVec()), trkposTruthSimHit);
+ System.out.printf("TruthSimHit residual %s for layer %d\n", resTruthSimHit.toString(), strip.layer());
+ }
+ }
+ }
+
//path length to intercept
- double s = HelixUtils.PathToXPlane(htf,trkpos.x(),0,0).get(0);
+ double s = HelixUtils.PathToXPlane(htf, trkpos.x(), 0, 0).get(0);
double s3D = s / Math.cos(Math.atan(htf.slope()));
- if(textFile != null) {
- textFile.printStripPathLen(s);
- textFile.printStripPathLen3D(s3D);
- }
-
+ if (textFile != null) {
+ textFile.printStripPathLen(s);
+ textFile.printStripPathLen3D(s3D);
+ }
+
//GBLDATA
stripData.setPath(s);
stripData.setPath3D(s3D);
-
-
-
+
//print stereo angle in YZ plane
- if(textFile != null) {
- textFile.printMeasDir(strip.u());
- textFile.printNonMeasDir(strip.v());
- textFile.printNormalDir(strip.w());
- }
-
+ if (textFile != null) {
+ textFile.printMeasDir(strip.u());
+ textFile.printNonMeasDir(strip.v());
+ textFile.printNormalDir(strip.w());
+ }
+
//GBLDATA
stripData.setU(strip.u());
stripData.setV(strip.v());
stripData.setW(strip.w());
-
-
+
//Print track direction at intercept
Hep3Vector tDir = HelixUtils.Direction(htf, s);
- double phi = htf.phi0() - s/htf.R();
+ double phi = htf.phi0() - s / htf.R();
double lambda = Math.atan(htf.slope());
- if(textFile != null) {
- textFile.printStripTrackDir(Math.sin(phi),Math.sin(lambda));
- textFile.printStripTrackDirFull(tDir);
- }
-
+ if (textFile != null) {
+ textFile.printStripTrackDir(Math.sin(phi), Math.sin(lambda));
+ textFile.printStripTrackDirFull(tDir);
+ }
+
//GBLDATA
stripData.setTrackDir(tDir);
stripData.setTrackPhi(phi);
stripData.setTrackLambda(lambda);
-
-
+
// calculate isolation to other strip clusters
double stripIsoMin = 9999.9;
for (SiTrackerHitStrip1D stripHit : stripHits) {
- if(stripHit.getRawHits().get(0).getDetectorElement().getName().equals(de.getName())) {
+ if (stripHit.getRawHits().get(0).getDetectorElement().getName().equals(de.getName())) {
SiTrackerHitStrip1D local = stripHit.getTransformedHit(TrackerHitType.CoordinateSystem.SENSOR);
double d = Math.abs(strip.umeas() - local.getPosition()[0]);
- if (d<stripIsoMin && d>0) {
+ if (d < stripIsoMin && d > 0) {
stripIsoMin = d;
}
}
}
-
- if(_debug>0) System.out.printf("%s: stripIsoMin = %f \n", this.getClass().getSimpleName(), stripIsoMin);
-
+
+ if (_debug > 0) {
+ System.out.printf("%s: stripIsoMin = %f \n", this.getClass().getSimpleName(), stripIsoMin);
+ }
+
// Add isolation to text file output
- if(textFile != null) {
+ if (textFile != null) {
textFile.printStripIso(stripIsoMin);
}
-
-
+
//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 = htfTruth!=null ? VecOp.sub(trkposTruth, origin): null;
-
+ Hep3Vector vdiffTrkTruth = htfTruth != null ? VecOp.sub(trkposTruth, origin) : null;
+
// then find the rotation from tracking to measurement frame
Hep3Matrix trkToStripRot = _trackerHitUtils.getTrackToStripRotation(strip.getStrip());
-
+
// then rotate that vector into the measurement frame to get the predicted measurement position
Hep3Vector trkpos_meas = VecOp.mult(trkToStripRot, vdiffTrk);
- Hep3Vector trkposTruth_meas = vdiffTrkTruth!=null ? VecOp.mult(trkToStripRot, vdiffTrkTruth) : null;
-
-
+ Hep3Vector trkposTruth_meas = vdiffTrkTruth != null ? VecOp.mult(trkToStripRot, vdiffTrkTruth) : null;
+
// 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));
-
- if(textFile != null) {
- textFile.printStripMeas(m_meas.x());
+ 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));
+
+ if (textFile != null) {
+ textFile.printStripMeas(m_meas.x());
}
//if(textFile != null) {
// textFile.printStripTrackPosMeasFrame(trkpos_meas);
//}
-
-
//GBLDATA
stripData.setMeas(strip.umeas());
stripData.setTrackPos(trkpos_meas);
-
- if(_debug>1) {
- System.out.printf("%s: rotation matrix to meas frame\n%s\n", getClass().getSimpleName(), VecOp.toString(trkToStripRot));
- System.out.printf("%s: tPosGlobal %s origin %s\n", getClass().getSimpleName(), trkpos.toString(), origin.toString());
- System.out.printf("%s: tDiff %s\n", getClass().getSimpleName(), vdiffTrk.toString());
- System.out.printf("%s: tPosMeas %s\n", getClass().getSimpleName(), trkpos_meas.toString());
- }
-
-
-
+
+ if (_debug > 1) {
+ System.out.printf("%s: rotation matrix to meas frame\n%s\n", getClass().getSimpleName(), VecOp.toString(trkToStripRot));
+ System.out.printf("%s: tPosGlobal %s origin %s\n", getClass().getSimpleName(), trkpos.toString(), origin.toString());
+ System.out.printf("%s: tDiff %s\n", getClass().getSimpleName(), vdiffTrk.toString());
+ System.out.printf("%s: tPosMeas %s\n", getClass().getSimpleName(), trkpos_meas.toString());
+ }
+
// residual in measurement frame
Hep3Vector res_meas = VecOp.sub(m_meas, trkpos_meas);
Hep3Vector resTruth_meas = trkposTruth_meas != null ? VecOp.sub(m_meas, trkposTruth_meas) : null;
- if(textFile != null) {
- textFile.printStripMeasRes(res_meas.x(),res_err_meas.x());
- textFile.printStripMeasResTruth(resTruth_meas!=null ? resTruth_meas.x() : -9999999.9,res_err_meas.x());
- }
-
+ if (textFile != null) {
+ textFile.printStripMeasRes(res_meas.x(), res_err_meas.x());
+ textFile.printStripMeasResTruth(resTruth_meas != null ? resTruth_meas.x() : -9999999.9, res_err_meas.x());
+ }
+
//GBLDATA
stripData.setMeasErr(res_err_meas.x());
-
-
-
- if(_debug>0) System.out.printf("layer %d millePedeId %d uRes %.10f\n",strip.layer(), millepedeId ,res_meas.x());
-
+
+ if (_debug > 0) {
+ System.out.printf("layer %d millePedeId %d uRes %.10f\n", strip.layer(), millepedeId, res_meas.x());
+ }
+
// sim hit residual
-
- if(simHit!=null) {
+ if (simHit != null) {
Hep3Vector simHitPos = CoordinateTransformations.transformVectorToTracking(simHit.getPositionVec());
- if(_debug>0) System.out.printf("simHitPos %s\n",simHitPos.toString());
+ if (_debug > 0) {
+ System.out.printf("simHitPos %s\n", simHitPos.toString());
+ }
Hep3Vector vdiffSimHit = VecOp.sub(simHitPos, trkpos);
Hep3Vector simHitPos_meas = VecOp.mult(trkToStripRot, vdiffSimHit);
- if(textFile != null) {
- textFile.printStripMeasResSimHit(simHitPos_meas.x(),res_err_meas.x());
+ if (textFile != null) {
+ textFile.printStripMeasResSimHit(simHitPos_meas.x(), res_err_meas.x());
}
} else {
- if(textFile != null) {
- textFile.printStripMeasResSimHit(-999999.9,-999999.9);
- }
- }
-
+ if (textFile != null) {
+ textFile.printStripMeasResSimHit(-999999.9, -999999.9);
+ }
+ }
+
// find scattering angle
- ScatterPoint scatter = scatters.getScatterPoint(((RawTrackerHit)strip.getStrip().rawhits().get(0)).getDetectorElement());
+ ScatterPoint scatter = scatters.getScatterPoint(((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement());
double scatAngle;
-
- if(scatter != null) {
+
+ if (scatter != null) {
scatAngle = scatter.getScatterAngle().Angle();
- }
- else {
+ } else {
if (_debug > 0) {
System.out.printf("%s: WARNING cannot find scatter for detector %s with strip cluster at %s\n", this.getClass(), ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement().getName(), strip.origin().toString());
}
//can be edge case where helix is outside, but close to sensor, so make a new scatter point assuming the helix does pass through the sensor
DetectorPlane hitPlane = null;
- if(MaterialSupervisor.class.isInstance(_scattering.getMaterialManager())) {
- MaterialSupervisor matSup = (MaterialSupervisor)_scattering.getMaterialManager();
- IDetectorElement hitElement = ((RawTrackerHit)strip.getStrip().rawhits().get(0)).getDetectorElement();
- for(ScatteringDetectorVolume vol : matSup.getMaterialVolumes()) {
- if (vol.getDetectorElement()==hitElement) {
+ if (MaterialSupervisor.class.isInstance(_scattering.getMaterialManager())) {
+ MaterialSupervisor matSup = (MaterialSupervisor) _scattering.getMaterialManager();
+ IDetectorElement hitElement = ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement();
+ for (ScatteringDetectorVolume vol : matSup.getMaterialVolumes()) {
+ if (vol.getDetectorElement() == hitElement) {
hitPlane = (DetectorPlane) vol;
break;
}
}
- if(hitPlane==null) {
+ if (hitPlane == null) {
throw new RuntimeException("cannot find plane for hit!");
} else {
// find scatterlength
- double s_closest = HelixUtils.PathToXPlane(htf,hitPlane.origin().x(), 0., 0).get(0);
+ double s_closest = HelixUtils.PathToXPlane(htf, hitPlane.origin().x(), 0., 0).get(0);
double X0 = hitPlane.getMaterialTraversedInRL(HelixUtils.Direction(htf, s_closest));
- ScatterAngle scatterAngle = new ScatterAngle(s_closest, _scattering.msangle(htf.p(this._B.magnitude()),X0));
+ ScatterAngle scatterAngle = new ScatterAngle(s_closest, _scattering.msangle(htf.p(this._B.magnitude()), X0));
scatAngle = scatterAngle.Angle();
}
- }
- else {
+ } else {
throw new UnsupportedOperationException("Should not happen. This problem is only solved with the MaterialSupervisor.");
}
- }
-
-
+ }
+
//print scatterer to file
- if(textFile != null) {
- textFile.printStripScat(scatAngle);
- }
-
+ if (textFile != null) {
+ textFile.printStripScat(scatAngle);
+ }
+
//GBLDATA
stripData.setScatterAngle(scatAngle);
-
+
++istrip;
-
-
-
- }
-
- }
-
-
-
- }
-
-
+ }
+ }
+ }
+
private void checkAprimeTruth(MCParticle mcp, List<MCParticle> mcParticles) {
- 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(TrackUtils.getHTF(mcp_pair.get(0),-1.0*this._B.z()),TrackUtils.getHTF(mcp_pair.get(1),-1.0*this._B.z()));
- 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());
- }
- }
-
- }
-
+ 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(TrackUtils.getHTF(mcp_pair.get(0), -1.0 * this._B.z()), TrackUtils.getHTF(mcp_pair.get(1), -1.0 * this._B.z()));
+ 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;
-
- Map<MCParticle,Integer> particlesOnTrack = new HashMap<MCParticle,Integer>();
-
- if(debug) System.out.printf("getmatched mc particle from %d tracker hits on the track \n",track.getTrackerHits().size());
-
-
- for(TrackerHit hit : track.getTrackerHits()) {
- List<MCParticle> mcps = ((HelicalTrackHit)hit).getMCParticles();
- if(mcps==null) {
- System.out.printf("%s: warning, this hit (layer %d pos=%s) has no mc particles.\n",this.getClass().getSimpleName(),((HelicalTrackHit)hit).Layer(),((HelicalTrackHit)hit).getCorrectedPosition().toString());
- }
- else {
- if( debug ) System.out.printf("%s: this hit (layer %d pos=%s) has %d mc particles.\n",this.getClass().getSimpleName(),((HelicalTrackHit)hit).Layer(),((HelicalTrackHit)hit).getCorrectedPosition().toString(),mcps.size());
- for(MCParticle mcp : mcps) {
- if( !particlesOnTrack.containsKey(mcp) ) {
+
+ Map<MCParticle, Integer> particlesOnTrack = new HashMap<MCParticle, Integer>();
+
+ if (debug) {
+ System.out.printf("getmatched mc particle from %d tracker hits on the track \n", track.getTrackerHits().size());
+ }
+
+ for (TrackerHit hit : track.getTrackerHits()) {
+ List<MCParticle> mcps = ((HelicalTrackHit) hit).getMCParticles();
+ if (mcps == null) {
+ System.out.printf("%s: warning, this hit (layer %d pos=%s) has no mc particles.\n", this.getClass().getSimpleName(), ((HelicalTrackHit) hit).Layer(), ((HelicalTrackHit) hit).getCorrectedPosition().toString());
+ } else {
+ if (debug) {
+ System.out.printf("%s: this hit (layer %d pos=%s) has %d mc particles.\n", this.getClass().getSimpleName(), ((HelicalTrackHit) hit).Layer(), ((HelicalTrackHit) hit).getCorrectedPosition().toString(), mcps.size());
+ }
+ for (MCParticle mcp : mcps) {
+ if (!particlesOnTrack.containsKey(mcp)) {
particlesOnTrack.put(mcp, 0);
}
int c = particlesOnTrack.get(mcp);
- particlesOnTrack.put(mcp, c+1);
- }
- }
- }
- if(debug) {
- System.out.printf("Track p=[ %f, %f, %f] \n",track.getTrackStates().get(0).getMomentum()[0],track.getTrackStates().get(0).getMomentum()[1],track.getTrackStates().get(0).getMomentum()[1]);
- System.out.printf("FOund %d particles\n",particlesOnTrack.size());
- for(Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet()) {
- System.out.printf("%d hits assigned to %d p=%s \n",entry.getValue(),entry.getKey().getPDGID(),entry.getKey().getMomentum().toString());
- }
- }
- Map.Entry<MCParticle,Integer> maxEntry = null;
- for(Map.Entry<MCParticle,Integer> entry : particlesOnTrack.entrySet()) {
- if ( maxEntry == null || entry.getValue().compareTo(maxEntry.getValue()) > 0 ) maxEntry = entry;
+ particlesOnTrack.put(mcp, c + 1);
+ }
+ }
+ }
+ if (debug) {
+ System.out.printf("Track p=[ %f, %f, %f] \n", track.getTrackStates().get(0).getMomentum()[0], track.getTrackStates().get(0).getMomentum()[1], track.getTrackStates().get(0).getMomentum()[1]);
+ System.out.printf("FOund %d particles\n", particlesOnTrack.size());
+ for (Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet()) {
+ System.out.printf("%d hits assigned to %d p=%s \n", entry.getValue(), entry.getKey().getPDGID(), entry.getKey().getMomentum().toString());
+ }
+ }
+ Map.Entry<MCParticle, Integer> maxEntry = null;
+ for (Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet()) {
+ if (maxEntry == null || entry.getValue().compareTo(maxEntry.getValue()) > 0) {
+ maxEntry = entry;
+ }
//if ( maxEntry != null ) {
// if(entry.getValue().compareTo(maxEntry.getValue()) < 0) continue;
//}
//maxEntry = entry;
}
- if(debug) {
- if (maxEntry != null ) {
- System.out.printf("Matched particle with pdgId=%d and mom %s to track with charge %d and momentum [%f %f %f]\n",
- maxEntry.getKey().getPDGID(),maxEntry.getKey().getMomentum().toString(),
- track.getCharge(),track.getTrackStates().get(0).getMomentum()[0],track.getTrackStates().get(0).getMomentum()[1],track.getTrackStates().get(0).getMomentum()[2]);
- } else {
- System.out.printf("No truth particle found on this track\n");
- }
+ if (debug) {
+ if (maxEntry != null) {
+ System.out.printf("Matched particle with pdgId=%d and mom %s to track with charge %d and momentum [%f %f %f]\n",
+ maxEntry.getKey().getPDGID(), maxEntry.getKey().getMomentum().toString(),
+ track.getCharge(), track.getTrackStates().get(0).getMomentum()[0], track.getTrackStates().get(0).getMomentum()[1], track.getTrackStates().get(0).getMomentum()[2]);
+ } else {
+ System.out.printf("No truth particle found on this track\n");
+ }
}
return maxEntry == null ? null : maxEntry.getKey();
}
-
-
- private BasicMatrix getPerParVector(HelicalTrackFit htf) {
- BasicMatrix perPar = new BasicMatrix(1,5);
- if( htf != null) {
- 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);
- perPar.setElement(0,1,theta);
- perPar.setElement(0,2,htf.phi0());
- perPar.setElement(0,3,htf.dca());
- perPar.setElement(0,4,htf.z0());
- }
- return perPar;
-
- }
-
-
-
-
- private BasicMatrix getJacPerToCl(HelicalTrackFit htf) {
- System.out.printf("%s: getJacPerToCl\n",this.getClass().getSimpleName());
- //use propoerly normalized B-field
- Hep3Vector Bnorm = VecOp.mult(Constants.fieldConversion, _B);
- //init jacobian to zero
- BasicMatrix j = new BasicMatrix(5,5);
- initZero(j);
- double lambda = Math.atan(htf.slope());
- double q = Math.signum(htf.R());
- double theta = Math.PI/2.0 - lambda;
- Hep3Vector T = HelixUtils.Direction(htf, 0.);
- Hep3Vector p = VecOp.mult(htf.p(Math.abs(_B.z())), T);
- double pT = htf.pT(Math.abs(_B.z()));
- Hep3Vector H = VecOp.mult(1./(Bnorm.magnitude()), Bnorm);
- Hep3Vector Z = new BasicHep3Vector(0,0,1);
- Hep3Vector J = VecOp.mult(1./VecOp.cross(T,Z).magnitude(), VecOp.cross(T, Z));
- Hep3Vector U = VecOp.mult(-1, J);
- Hep3Vector V = VecOp.cross(T, U);
- double alpha = VecOp.cross(H,T).magnitude();
- Hep3Vector N = VecOp.mult(1./alpha,VecOp.cross(H, T));
- Hep3Vector K = Z;
- double Q = -Bnorm.magnitude()*q/p.magnitude();
- double kappa = -1.0*q*Bnorm.z()/pT;
-
- if (this._debug!=0) {
- System.out.printf("%s: Bnorm=%s mag(Bnorm)=%f\n",this.getClass().getSimpleName(),Bnorm.toString(),Bnorm.magnitude());
- System.out.printf("%s: p=%s |p|=%f pT=%f\n",this.getClass().getSimpleName(),p.toString(),p.magnitude(),pT);
- System.out.printf("%s: q=%f\n",this.getClass().getSimpleName(),q);
- System.out.printf("%s: q/p=%f\n",this.getClass().getSimpleName(),q/p.magnitude());
- System.out.printf("%s: T=%s\n",this.getClass().getSimpleName(),T.toString());
- System.out.printf("%s: H=%s\n",this.getClass().getSimpleName(),H.toString());
- System.out.printf("%s: kappa=%f\n",this.getClass().getSimpleName(),kappa);
- System.out.printf("%s: alpha=%f Q=%f \n",this.getClass().getSimpleName(),alpha,Q);
- System.out.printf("%s: J=%s \n",this.getClass().getSimpleName(),J.toString());
- System.out.printf("%s: V=%s \n",this.getClass().getSimpleName(),V.toString());
- System.out.printf("%s: N=%s \n",this.getClass().getSimpleName(),N.toString());
- System.out.printf("%s: TdotJ=%f \n",this.getClass().getSimpleName(),VecOp.dot(T, J));
- System.out.printf("%s: VdotN=%f \n",this.getClass().getSimpleName(),VecOp.dot(V, N));
- System.out.printf("%s: TdotK=%f \n",this.getClass().getSimpleName(),VecOp.dot(T, K));
- System.out.printf("%s: UdotN=%f \n",this.getClass().getSimpleName(),VecOp.dot(U, N));
- }
-
-
-
-
-
- j.setElement(0,0,-1.0*Math.sin(theta)/Bnorm.z());
-
- j.setElement(0,1,q/(p.magnitude()*Math.tan(theta)));
-
- j.setElement(1,1,-1);
-
- j.setElement(1,3,-alpha*Q*VecOp.dot(T, J)*VecOp.dot(V, N));
-
- j.setElement(1,4,-alpha*Q*VecOp.dot(T, K)*VecOp.dot(V, N));
-
- j.setElement(2,2,1);
-
- j.setElement(2,3,-alpha*Q*VecOp.dot(T,J)*VecOp.dot(U, N)/Math.cos(lambda));
-
- j.setElement(2,4,-alpha*Q*VecOp.dot(T,K)*VecOp.dot(U, N)/Math.cos(lambda));
-
- j.setElement(3,3,-1);
-
- j.setElement(4,4,VecOp.dot(V, K));
-
-
- if(_debug>0) {
- System.out.printf("%s: lambda= J(1,1)=%f * theta + J(1,3)=%f * eps + J(1,4)=%f * z0 \n",
- this.getClass().getSimpleName(),
- j.e(1, 1),j.e(1,3),j.e(1,4));
-
- }
-
-
- return j;
-
- }
-
-
-
-
-
-
+
+// private BasicMatrix getJacPerToCl(HelicalTrackFit htf) {
+// System.out.printf("%s: getJacPerToCl\n", this.getClass().getSimpleName());
+// //use propoerly normalized B-field
+// Hep3Vector Bnorm = VecOp.mult(Constants.fieldConversion, _B);
+// //init jacobian to zero
+// BasicMatrix j = new BasicMatrix(5, 5);
+// initZero(j);
+// double lambda = Math.atan(htf.slope());
+// double q = Math.signum(htf.R());
+// double theta = Math.PI / 2.0 - lambda;
+// Hep3Vector T = HelixUtils.Direction(htf, 0.);
+// Hep3Vector p = VecOp.mult(htf.p(Math.abs(_B.z())), T);
+// double pT = htf.pT(Math.abs(_B.z()));
+// Hep3Vector H = VecOp.mult(1. / (Bnorm.magnitude()), Bnorm);
+// Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
+// Hep3Vector J = VecOp.mult(1. / VecOp.cross(T, Z).magnitude(), VecOp.cross(T, Z));
+// Hep3Vector U = VecOp.mult(-1, J);
+// Hep3Vector V = VecOp.cross(T, U);
+// double alpha = VecOp.cross(H, T).magnitude();
+// Hep3Vector N = VecOp.mult(1. / alpha, VecOp.cross(H, T));
+// Hep3Vector K = Z;
+// double Q = -Bnorm.magnitude() * q / p.magnitude();
+// double kappa = -1.0 * q * Bnorm.z() / pT;
+//
+// if (this._debug != 0) {
+// System.out.printf("%s: Bnorm=%s mag(Bnorm)=%f\n", this.getClass().getSimpleName(), Bnorm.toString(), Bnorm.magnitude());
+// System.out.printf("%s: p=%s |p|=%f pT=%f\n", this.getClass().getSimpleName(), p.toString(), p.magnitude(), pT);
+// System.out.printf("%s: q=%f\n", this.getClass().getSimpleName(), q);
+// System.out.printf("%s: q/p=%f\n", this.getClass().getSimpleName(), q / p.magnitude());
+// System.out.printf("%s: T=%s\n", this.getClass().getSimpleName(), T.toString());
+// System.out.printf("%s: H=%s\n", this.getClass().getSimpleName(), H.toString());
+// System.out.printf("%s: kappa=%f\n", this.getClass().getSimpleName(), kappa);
+// System.out.printf("%s: alpha=%f Q=%f \n", this.getClass().getSimpleName(), alpha, Q);
+// System.out.printf("%s: J=%s \n", this.getClass().getSimpleName(), J.toString());
+// System.out.printf("%s: V=%s \n", this.getClass().getSimpleName(), V.toString());
+// System.out.printf("%s: N=%s \n", this.getClass().getSimpleName(), N.toString());
+// System.out.printf("%s: TdotJ=%f \n", this.getClass().getSimpleName(), VecOp.dot(T, J));
+// System.out.printf("%s: VdotN=%f \n", this.getClass().getSimpleName(), VecOp.dot(V, N));
+// System.out.printf("%s: TdotK=%f \n", this.getClass().getSimpleName(), VecOp.dot(T, K));
+// System.out.printf("%s: UdotN=%f \n", this.getClass().getSimpleName(), VecOp.dot(U, N));
+// }
+//
+// j.setElement(0, 0, -1.0 * Math.sin(theta) / Bnorm.z());
+//
+// j.setElement(0, 1, q / (p.magnitude() * Math.tan(theta)));
+//
+// j.setElement(1, 1, -1);
+//
+// j.setElement(1, 3, -alpha * Q * VecOp.dot(T, J) * VecOp.dot(V, N));
+//
+// j.setElement(1, 4, -alpha * Q * VecOp.dot(T, K) * VecOp.dot(V, N));
+//
+// j.setElement(2, 2, 1);
+//
+// j.setElement(2, 3, -alpha * Q * VecOp.dot(T, J) * VecOp.dot(U, N) / Math.cos(lambda));
+//
+// j.setElement(2, 4, -alpha * Q * VecOp.dot(T, K) * VecOp.dot(U, N) / Math.cos(lambda));
+//
+// j.setElement(3, 3, -1);
+//
+// j.setElement(4, 4, VecOp.dot(V, K));
+//
+// if (_debug > 0) {
+// System.out.printf("%s: lambda= J(1,1)=%f * theta + J(1,3)=%f * eps + J(1,4)=%f * z0 \n",
+// this.getClass().getSimpleName(),
+// j.e(1, 1), j.e(1, 3), j.e(1, 4));
+//
+// }
+//
+// return j;
+//
+// }
private void initUnit(BasicMatrix mat) {
- for(int row=0;row!=mat.getNRows();row++) {
- for(int col=0;col!=mat.getNColumns();col++) {
- if(row!=col) mat.setElement(row, col, 0);
- else mat.setElement(row, col, 1);
+ for (int row = 0; row != mat.getNRows(); row++) {
+ for (int col = 0; col != mat.getNColumns(); col++) {
+ if (row != col) {
+ mat.setElement(row, col, 0);
+ } else {
+ mat.setElement(row, col, 1);
+ }
}
}
}
private void initZero(BasicMatrix mat) {
- for(int row=0;row!=mat.getNRows();row++) {
- for(int col=0;col!=mat.getNColumns();col++) {
+ for (int row = 0; row != mat.getNRows(); row++) {
+ for (int col = 0; col != mat.getNColumns(); col++) {
mat.setElement(row, col, 0);
}
}
}
-
-
/**
- * Transform MCParticle into a Helix object.
- * Note that it produces the helix parameters at nominal x=0 and assumes that there is no field at x<0
- *
+ * Transform MCParticle into a Helix object. Note that it produces the helix
+ * parameters at nominal x=0 and assumes that there is no field at x<0
+ *
* @param mcp MC particle to be transformed
* @return helix object based on the MC particle
*/
@@ -771,82 +733,87 @@
// HelicalTrackFit htf = new HelicalTrackFit(par, cov, new double[2], new int[2], null, null);
// return htf;
// }
-
private double truthTrackFitChi2(PerigeeParams perPar, PerigeeParams perParTruth, SymmetricMatrix covariance) {
//re-shuffle the param vector to match the covariance order of parameters
- BasicMatrix p = new BasicMatrix(1,5);
+ BasicMatrix p = new BasicMatrix(1, 5);
p.setElement(0, 0, perPar.getD0());
p.setElement(0, 1, perPar.getPhi());
p.setElement(0, 2, perPar.getKappa());
p.setElement(0, 0, perPar.getZ0());
- p.setElement(0, 4, Math.tan(Math.PI/2.0-perPar.getTheta()));
- BasicMatrix pt = new BasicMatrix(1,5);
+ p.setElement(0, 4, Math.tan(Math.PI / 2.0 - perPar.getTheta()));
+ BasicMatrix pt = new BasicMatrix(1, 5);
pt.setElement(0, 0, perParTruth.getD0());
pt.setElement(0, 1, perParTruth.getPhi());
pt.setElement(0, 2, perParTruth.getKappa());
pt.setElement(0, 0, perParTruth.getZ0());
- pt.setElement(0, 4, Math.tan(Math.PI/2.0-perParTruth.getTheta()));
+ 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);
- BasicMatrix chi2 = (BasicMatrix) MatrixOp.mult(res,MatrixOp.mult(error_matrix, MatrixOp.transposed(res)));
- if(chi2.getNColumns()!=1 ||chi2.getNRows()!=1) {
+ 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!");
}
return chi2.e(0, 0);
}
-
-
private List<MCParticle> getAprimeDecayProducts(List<MCParticle> mcParticles) {
List<MCParticle> pair = new ArrayList<MCParticle>();
- for(MCParticle mcp : mcParticles) {
- if(mcp.getGeneratorStatus()!=MCParticle.FINAL_STATE) continue;
+ for (MCParticle mcp : mcParticles) {
+ if (mcp.getGeneratorStatus() != MCParticle.FINAL_STATE) {
+ continue;
+ }
boolean hasAprimeParent = false;
- for(MCParticle parent : mcp.getParents()) {
- if(Math.abs(parent.getPDGID())==622) hasAprimeParent = true;
- }
- if(hasAprimeParent) pair.add(mcp);
- }
- if(pair.size()!=2) {
- System.out.printf("%s: ERROR this event has %d mcp with 622 as parent!!?? \n",this.getClass().getSimpleName(),pair.size());
+ for (MCParticle parent : mcp.getParents()) {
+ if (Math.abs(parent.getPDGID()) == 622) {
+ hasAprimeParent = true;
+ }
+ }
+ if (hasAprimeParent) {
+ pair.add(mcp);
+ }
+ }
+ if (pair.size() != 2) {
+ System.out.printf("%s: ERROR this event has %d mcp with 622 as parent!!?? \n", this.getClass().getSimpleName(), pair.size());
this.printMCParticles(mcParticles);
System.exit(1);
}
- if( Math.abs(pair.get(0).getPDGID()) != 11 || Math.abs(pair.get(1).getPDGID()) != 11 ) {
- System.out.printf("%s: ERROR decay products are not e+e-? \n",this.getClass().getSimpleName());
+ if (Math.abs(pair.get(0).getPDGID()) != 11 || Math.abs(pair.get(1).getPDGID()) != 11) {
+ System.out.printf("%s: ERROR decay products are not e+e-? \n", this.getClass().getSimpleName());
this.printMCParticles(mcParticles);
System.exit(1);
}
- if(pair.get(0).getPDGID()*pair.get(1).getPDGID() > 0) {
- System.out.printf("%s: ERROR decay products have the same sign? \n",this.getClass().getSimpleName());
+ if (pair.get(0).getPDGID() * pair.get(1).getPDGID() > 0) {
+ System.out.printf("%s: ERROR decay products have the same sign? \n", this.getClass().getSimpleName());
this.printMCParticles(mcParticles);
System.exit(1);
}
return pair;
-
- }
-
- private void printMCParticles(List<MCParticle> mcParticles) {
- System.out.printf("%s: printMCParticles \n",this.getClass().getSimpleName());
- System.out.printf("%s: %d mc particles \n",this.getClass().getSimpleName(),mcParticles.size());
- for(MCParticle mcp : mcParticles) {
- if(mcp.getGeneratorStatus()!=MCParticle.FINAL_STATE) continue;
- System.out.printf("\n%s: (%s) %d p %s org %s %s \n",this.getClass().getSimpleName(),
- mcp.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",mcp.getPDGID(),mcp.getMomentum().toString(),mcp.getOrigin().toString(),
- mcp.getParents().size()>0?"parents:":"");
- for(MCParticle parent : mcp.getParents()) {
- System.out.printf("%s: (%s) %d p %s org %s %s \n",this.getClass().getSimpleName(),
- parent.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",parent.getPDGID(),parent.getMomentum().toString(),parent.getOrigin().toString(),
- parent.getParents().size()>0?"parents:":"");
- for(MCParticle grparent : parent.getParents()) {
- System.out.printf("%s: (%s) %d p %s org %s %s \n",this.getClass().getSimpleName(),
- grparent.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",grparent.getPDGID(),grparent.getMomentum().toString(),grparent.getOrigin().toString(),
- grparent.getParents().size()>0?"parents:":"");
- }
-
- }
- }
- }
+
+ }
+
+ private void printMCParticles(List<MCParticle> mcParticles) {
+ System.out.printf("%s: printMCParticles \n", this.getClass().getSimpleName());
+ System.out.printf("%s: %d mc particles \n", this.getClass().getSimpleName(), mcParticles.size());
+ for (MCParticle mcp : mcParticles) {
+ if (mcp.getGeneratorStatus() != MCParticle.FINAL_STATE) {
+ continue;
+ }
+ System.out.printf("\n%s: (%s) %d p %s org %s %s \n", this.getClass().getSimpleName(),
+ mcp.getGeneratorStatus() == MCParticle.FINAL_STATE ? "F" : "I", mcp.getPDGID(), mcp.getMomentum().toString(), mcp.getOrigin().toString(),
+ mcp.getParents().size() > 0 ? "parents:" : "");
+ for (MCParticle parent : mcp.getParents()) {
+ System.out.printf("%s: (%s) %d p %s org %s %s \n", this.getClass().getSimpleName(),
+ parent.getGeneratorStatus() == MCParticle.FINAL_STATE ? "F" : "I", parent.getPDGID(), parent.getMomentum().toString(), parent.getOrigin().toString(),
+ parent.getParents().size() > 0 ? "parents:" : "");
+ for (MCParticle grparent : parent.getParents()) {
+ System.out.printf("%s: (%s) %d p %s org %s %s \n", this.getClass().getSimpleName(),
+ grparent.getGeneratorStatus() == MCParticle.FINAL_STATE ? "F" : "I", grparent.getPDGID(), grparent.getMomentum().toString(), grparent.getOrigin().toString(),
+ grparent.getParents().size() > 0 ? "parents:" : "");
+ }
+
+ }
+ }
+ }
private double getInvMassTracks(HelicalTrackFit htf1, HelicalTrackFit htf2) {
double p1 = htf1.p(this._B.magnitude());
@@ -854,55 +821,77 @@
Hep3Vector p1vec = VecOp.mult(p1, HelixUtils.Direction(htf1, 0));
Hep3Vector p2vec = VecOp.mult(p2, HelixUtils.Direction(htf2, 0));
double me = 0.000510998910;
- double E1 = Math.sqrt(p1*p1 + me*me);
- double E2 = Math.sqrt(p2*p2 + me*me);
+ double E1 = Math.sqrt(p1 * p1 + me * me);
+ double E2 = Math.sqrt(p2 * p2 + me * me);
//System.out.printf("p1 %f %s E1 %f\n",p1,p1vec.toString(),E1);
//System.out.printf("p2 %f %s E2 %f\n",p2,p2vec.toString(),E2);
- return Math.sqrt( Math.pow(E1+E2,2) - VecOp.add(p1vec, p2vec).magnitudeSquared() );
- }
-
-
+ return Math.sqrt(Math.pow(E1 + E2, 2) - VecOp.add(p1vec, p2vec).magnitudeSquared());
+ }
+
+ private BasicMatrix getPerParVector(HelicalTrackFit htf) {
+ BasicMatrix perPar = new BasicMatrix(1, 5);
+ if (htf != null) {
+ 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);
+ perPar.setElement(0, 1, theta);
+ perPar.setElement(0, 2, htf.phi0());
+ perPar.setElement(0, 3, htf.dca());
+ perPar.setElement(0, 4, htf.z0());
+ }
+ return perPar;
+
+ }
+
public class PerigeeParams {
+
private final BasicMatrix _params;
private PerigeeParams(HelicalTrackFit htf) {
_params = GBLOutput.this.getPerParVector(htf);
}
+
public BasicMatrix getParams() {
return _params;
}
+
public double getKappa() {
return _params.e(0, 0);
}
+
public double getTheta() {
return _params.e(0, 1);
}
+
public double getPhi() {
return _params.e(0, 2);
}
+
public double getD0() {
return _params.e(0, 3);
}
+
public double getZ0() {
return _params.e(0, 4);
}
}
/**
- * Computes the projection matrix from the perigee XY plane variables
- * dca and z0 into the curvilinear xT,yT,zT frame (U,V,T)
+ * Computes the projection matrix from the perigee XY plane variables dca
+ * and z0 into the curvilinear xT,yT,zT frame (U,V,T)
+ *
* @param htf input helix to find the track direction
* @return 3x3 projection matrix
*/
Hep3Matrix getPerToClPrj(HelicalTrackFit htf) {
- Hep3Vector Z = new BasicHep3Vector(0,0,1);
+ Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
Hep3Vector T = HelixUtils.Direction(htf, 0.);
- Hep3Vector J = VecOp.mult(1./VecOp.cross(T,Z).magnitude(), VecOp.cross(T, Z));
+ Hep3Vector J = VecOp.mult(1. / VecOp.cross(T, Z).magnitude(), VecOp.cross(T, Z));
Hep3Vector K = Z;
Hep3Vector U = VecOp.mult(-1, J);
Hep3Vector V = VecOp.cross(T, U);
Hep3Vector I = VecOp.cross(J, K);
-
+
BasicHep3Matrix trans = new BasicHep3Matrix();
trans.setElement(0, 0, VecOp.dot(I, U));
trans.setElement(0, 1, VecOp.dot(J, U));
@@ -914,67 +903,71 @@
trans.setElement(2, 1, VecOp.dot(J, T));
trans.setElement(2, 2, VecOp.dot(K, T));
return trans;
-
- }
-
-
-
-
- public class ClParams {
- private BasicMatrix _params = new BasicMatrix(1,5);
+
+ }
+
+ public class ClParams {
+
+ private BasicMatrix _params = new BasicMatrix(1, 5);
+
private ClParams(HelicalTrackFit htf) {
-
- if (htf == null) return;
+
+ if (htf == null) {
+ return;
+ }
Hep3Matrix perToClPrj = GBLOutput.this.getPerToClPrj(htf);
-
+
double d0 = -1 * htf.dca(); //sign convention for curvilinear frame
double z0 = htf.z0();
- Hep3Vector vecPer = new BasicHep3Vector(0.,d0,z0);
+ Hep3Vector vecPer = new BasicHep3Vector(0., d0, z0);
//System.out.printf("%s: vecPer=%s\n",this.getClass().getSimpleName(),vecPer.toString());
-
+
Hep3Vector vecCl = VecOp.mult(perToClPrj, vecPer);
//System.out.printf("%s: vecCl=%s\n",this.getClass().getSimpleName(),vecCl.toString());
double xT = vecCl.x();
double yT = vecCl.y();
//double zT = vecCl.z();
-
+
Hep3Vector T = HelixUtils.Direction(htf, 0.);
Hep3Vector p = VecOp.mult(htf.p(Math.abs(_B.z())), T);
double lambda = Math.atan(htf.slope());
double q = Math.signum(htf.R());
- double qOverP = q/p.magnitude();
+ double qOverP = q / p.magnitude();
double phi = htf.phi0();
-
+
_params.setElement(0, 0, qOverP);
_params.setElement(0, 1, lambda);
_params.setElement(0, 2, phi);
_params.setElement(0, 3, xT);
_params.setElement(0, 4, yT);
-
+
}
public BasicMatrix getParams() {
return _params;
}
-
+
double getLambda() {
- return _params.e(0,1);
- }
+ return _params.e(0, 1);
+ }
+
double getQoverP() {
- return _params.e(0,0);
- }
+ return _params.e(0, 0);
+ }
+
double getPhi() {
- return _params.e(0,2);
- }
+ return _params.e(0, 2);
+ }
+
double getXt() {
- return _params.e(0,3);
- }
+ return _params.e(0, 3);
+ }
+
double getYt() {
- return _params.e(0,4);
- }
-
- }
-
+ return _params.e(0, 4);
+ }
+
+ }
}
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java Fri Aug 28 16:36:48 2015
@@ -25,25 +25,24 @@
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
import org.lcsim.util.log.LogUtil;
-
/**
- * A class that creates track objects from fitted GBL trajectories and adds them into the event.
- *
+ * A class that creates track objects from fitted GBL trajectories and adds them
+ * into the event.
+ *
* @author Per Hansson Adrian <[log in to unmask]>
*
*/
public class MakeGblTracks {
-
private String _TrkCollectionName = "GblTracks";
private static Logger logger = LogUtil.create(MakeGblTracks.class, new BasicLogFormatter());
-
+
/**
* Creates a new instance of MakeTracks.
*/
public MakeGblTracks() {
//logger = Logger.getLogger(getClass().getName());
- //logger.setUseParentHandlers(false);
+ //logger.setUseParentHandlers(false);
//Handler handler = new StreamHandler(System.out, new SimpleFormatter());
//logger.addHandler(handler);
logger.setLevel(Level.INFO);
@@ -54,25 +53,25 @@
// }
}
-
+
/**
* Process a Gbl track and store it into the event
+ *
* @param event event header
* @param track Gbl trajectory
* @param seed SeedTrack
* @param bfield magnetic field (used to turn curvature into momentum)
*/
public void Process(EventHeader event, List<FittedGblTrajectory> gblTrajectories, double bfield) {
-
+
List<Track> tracks = new ArrayList<Track>();
-
+
logger.info("adding " + gblTrajectories.size() + " of fitted GBL tracks to the event");
-
- for(FittedGblTrajectory fittedTraj : gblTrajectories) {
-
-
+
+ for (FittedGblTrajectory fittedTraj : gblTrajectories) {
+
// Initialize the reference point to the origin
- double[] ref = new double[] {0., 0., 0.};
+ double[] ref = new double[]{0., 0., 0.};
SeedTrack seedTrack = (SeedTrack) fittedTraj.get_seed();
SeedCandidate trackseed = seedTrack.getSeedCandidate();
@@ -86,12 +85,12 @@
// Retrieve the helix and save the relevant bits of helix info
HelicalTrackFit helix = trackseed.getHelix();
- double gblParameters[] = getGblCorrectedHelixParameters(helix,fittedTraj, bfield);
+ double gblParameters[] = getGblCorrectedHelixParameters(helix, fittedTraj, bfield);
trk.setTrackParameters(gblParameters, bfield); // Sets first TrackState.
//TODO Use GBL covariance matrix
//SymmetricMatrix gblCovariance = getGblCorrectedHelixCovariance(helix, fittedTraj, bfield);
trk.setCovarianceMatrix(helix.covariance()); // Modifies first TrackState.
- trk.setChisq( fittedTraj.get_chi2());
+ trk.setChisq(fittedTraj.get_chi2());
trk.setNDF(fittedTraj.get_ndf());
// Flag that the fit was successful and set the reference point
@@ -107,22 +106,21 @@
// Check the track - hook for plugging in external constraint
//if ((_trackCheck != null) && (! _trackCheck.checkTrack(trk))) continue;
-
// Add the track to the list of tracks
tracks.add((Track) trk);
- logger.info(String.format("helix chi2 %f ndf %d gbl chi2 %f ndf %d\n", helix.chisqtot(), helix.ndf()[0]+helix.ndf()[1], trk.getChi2(), trk.getNDF()));
- if(logger.getLevel().intValue()<= Level.INFO.intValue()) {
- for(int i=0;i<5;++i) {
- logger.info(String.format("param %d: %.5f -> %.5f helix-gbl= %f", i, helix.parameters()[i], trk.getTrackParameter(i), helix.parameters()[i]-trk.getTrackParameter(i)));
+ logger.info(String.format("helix chi2 %f ndf %d gbl chi2 %f ndf %d\n", helix.chisqtot(), helix.ndf()[0] + helix.ndf()[1], trk.getChi2(), trk.getNDF()));
+ if (logger.getLevel().intValue() <= Level.INFO.intValue()) {
+ for (int i = 0; i < 5; ++i) {
+ logger.info(String.format("param %d: %.10f -> %.10f helix-gbl= %f", i, helix.parameters()[i], trk.getTrackParameter(i), helix.parameters()[i] - trk.getTrackParameter(i)));
}
}
-
+
}
logger.info("adding " + Integer.toString(tracks.size()) + " Gbl tracks to event with " + event.get(Track.class, "MatchedTracks").size() + " matched tracks");
-
+
// Put the tracks back into the event and exit
- int flag = 1<<LCIOConstants.TRBIT_HITS;
+ int flag = 1 << LCIOConstants.TRBIT_HITS;
event.put(_TrkCollectionName, tracks, Track.class, flag);
return;
@@ -130,6 +128,7 @@
/**
* Compute the track fit covariance matrix
+ *
* @param helix - original seed track
* @param traj - fitted GBL trajectory
* @return covariance matrix.
@@ -142,21 +141,22 @@
/**
* Compute the updated helix parameters.
+ *
* @param helix - original seed track
* @param traj - fitted GBL trajectory
* @return corrected parameters.
*/
private double[] getGblCorrectedHelixParameters(HelicalTrackFit helix, FittedGblTrajectory traj, double bfield) {
-
+
// get seed helix parameters
double d0 = helix.dca();
double z0 = helix.z0();
double phi0 = helix.phi0();
- double slope = helix.slope();
+ double lambda = Math.atan(helix.slope());
double p = helix.p(bfield);
double q = traj.get_seed().getCharge();
- double qOverP = q/p;
-
+ double qOverP = q / p;
+
// get corrections from GBL fit
Vector locPar = new Vector(5);
SymMatrix locCov = new SymMatrix(5);
@@ -166,53 +166,49 @@
double yTCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.YT.getValue());
double xTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
double yTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue());
-
- logger.info((slope>0?"top: ":"bot ") + "qOverPCorr " + qOverPCorr + " xTCorr " + xTCorr + " yTCorr " + yTCorr + " xtPrimeCorr " + xTPrimeCorr + " yTPrimeCorr " + yTPrimeCorr);
-
+
+ logger.info((helix.slope() > 0 ? "top: " : "bot ") + "qOverPCorr " + qOverPCorr + " xTCorr " + xTCorr + " yTCorr " + yTCorr + " xtPrimeCorr " + xTPrimeCorr + " yTPrimeCorr " + yTPrimeCorr);
+
// calculate new d0 and z0
Hep3Matrix perToClPrj = traj.get_track_data().getPrjPerToCl();
Hep3Matrix clToPerPrj = VecOp.inverse(perToClPrj);
Hep3Vector corrPer = VecOp.mult(clToPerPrj, new BasicHep3Vector(xTCorr, yTCorr, 0.0));
//d0
- double d0_corr = -1.0*corrPer.y(); // correct for different sign convention of d0 in curvilinear frame
+ double d0_corr = -1.0 * corrPer.y(); // correct for different sign convention of d0 in curvilinear frame
double d0_gbl = d0 + d0_corr;
-
+
//z0
double z0_corr = corrPer.z();
double z0_gbl = z0 + z0_corr;
-
+
//calculate new phi0
double phi0_gbl = phi0 + xTPrimeCorr;
-
+
//calculate new slope
- double lambda_gbl = Math.atan(slope) + yTPrimeCorr;
- double slope_gbl = Math.tan( lambda_gbl );
+ double lambda_gbl = lambda + yTPrimeCorr;
+ double slope_gbl = Math.tan(lambda_gbl);
// calculate new curvature
-
double qOverP_gbl = qOverP + qOverPCorr;
- double pt_gbl = Math.abs(1.0/qOverP_gbl) * Math.sin((Math.PI/2.0-lambda_gbl));
+ double pt_gbl = Math.abs(1.0 / qOverP_gbl) * Math.sin((Math.PI / 2.0 - lambda_gbl));
double C_gbl = Constants.fieldConversion * bfield / pt_gbl;
//make sure sign is not changed
- C_gbl = Math.signum(helix.curvature())*Math.abs(C_gbl);
-
- logger.info("qOverP="+qOverP+" qOverPCorr="+qOverPCorr+" qOverP_gbl="+qOverP_gbl+" ==> pGbl="+1.0/qOverP_gbl);
-
-
+ C_gbl = Math.signum(helix.curvature()) * Math.abs(C_gbl);
+
+ logger.info("qOverP=" + qOverP + " qOverPCorr=" + qOverPCorr + " qOverP_gbl=" + qOverP_gbl + " ==> pGbl=" + 1.0 / qOverP_gbl);
+
double parameters_gbl[] = new double[5];
parameters_gbl[HelicalTrackFit.dcaIndex] = d0_gbl;
- parameters_gbl[HelicalTrackFit.z0Index] = z0_gbl;
+ parameters_gbl[HelicalTrackFit.z0Index] = z0_gbl;
parameters_gbl[HelicalTrackFit.curvatureIndex] = C_gbl;
parameters_gbl[HelicalTrackFit.slopeIndex] = slope_gbl;
parameters_gbl[HelicalTrackFit.phi0Index] = phi0_gbl;
return parameters_gbl;
}
-
+
public void setTrkCollectionName(String name) {
_TrkCollectionName = name;
}
-
-
}
