Author: [log in to unmask] Date: Wed Sep 30 18:41:59 2015 New Revision: 3746 Log: some cleanup 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/GBLTrackData.java java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.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 Wed Sep 30 18:41:59 2015 @@ -69,8 +69,9 @@ */ GBLOutput(String outputFileName, Hep3Vector bfield) { //System.out.printf("name \"%s\" \n", outputFileName); - if (!outputFileName.equalsIgnoreCase("")) + if (!outputFileName.equalsIgnoreCase("")) { textFile = new GBLFileIO(outputFileName); + } _materialmanager = new MaterialSupervisor(); _scattering = new MultipleScattering(_materialmanager); _B = CoordinateTransformations.transformVectorToTracking(bfield); @@ -91,18 +92,21 @@ } void printNewEvent(int eventNumber, double Bz) { - if (textFile != null) + if (textFile != null) { textFile.printEventInfo(eventNumber, Bz); + } } void printTrackID(int iTrack) { - if (textFile != null) + if (textFile != null) { textFile.printTrackID(iTrack); + } } void close() { - if (textFile != null) + if (textFile != null) { textFile.closeFile(); + } } void setAPrimeEventFlag(boolean flag) { @@ -143,11 +147,13 @@ 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) + } else if (_debug > 0) { System.out.printf("%s: truth particle (pdgif %d ) found in event!\n", this.getClass().getSimpleName(), mcp.getPDGID()); - - if (AprimeEvent) + } + + if (AprimeEvent) { checkAprimeTruth(mcp, mcParticles); + } } // Get track parameters from MC particle @@ -158,13 +164,13 @@ // Get perigee parameters to curvilinear frame PerigeeParams perPar = new PerigeeParams(htf, _B.z()); PerigeeParams perParTruth = new PerigeeParams(htfTruth, _B.z()); - if (textFile != null) { - textFile.printPerTrackParam(perPar); - textFile.printPerTrackParamTruth(perParTruth); - } //GBLDATA gtd.setPerigeeTrackParameters(perPar); + if (textFile != null) { + textFile.printPerTrackParam(gtd.getPerigeeTrackParameters()); + textFile.printPerTrackParamTruth(perParTruth); + } // Get curvilinear parameters if (textFile != null) { @@ -182,37 +188,44 @@ // find the projection from the I,J,K to U,V,T curvilinear coordinates Hep3Matrix perToClPrj = getPerToClPrj(htf); - if (textFile != null) - textFile.printPerToClPrj(perToClPrj); - //GBLDATA - for (int row = 0; row < perToClPrj.getNRows(); ++row) - for (int col = 0; col < perToClPrj.getNColumns(); ++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)); + } + } + if (textFile != null) { + textFile.printPerToClPrj(gtd.getPrjPerToCl()); + } // print chi2 of fit - if (textFile != null) + 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>(); - for (SimTrackerHit sh : simTrackerHits) + for (SimTrackerHit sh : simTrackerHits) { if (sh.getMCParticle() == mcp) { int layer = sh.getIdentifierFieldValue("layer"); - if (!simHitsLayerMap.containsKey(layer) || (sh.getPathLength() < simHitsLayerMap.get(layer).getPathLength())) + if (!simHitsLayerMap.containsKey(layer) || (sh.getPathLength() < simHitsLayerMap.get(layer).getPathLength())) { simHitsLayerMap.put(layer, sh); - } + } + } + } // covariance matrix from the fit - if (textFile != null) + if (textFile != null) { textFile.printPerTrackCov(htf); + } // dummy cov matrix for CL parameters BasicMatrix clCov = GblUtils.unitMatrix(5, 5); clCov = (BasicMatrix) MatrixOp.mult(0.1 * 0.1, clCov); - if (textFile != null) + 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()); @@ -220,8 +233,9 @@ System.out.printf("%s: truth perPar chi2 %f\n", this.getClass().getSimpleName(), chi2_truth); } - if (_debug > 0) + if (_debug > 0) { System.out.printf("%d hits\n", hits.size()); + } int istrip = 0; for (int ihit = 0; ihit != hits.size(); ++ihit) { @@ -236,44 +250,46 @@ // 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 + } else { throw new ClassCastException("Detector element " + de.getName() + " couldn't be casted to " + HpsSiSensor.class.getName()); + } int millepedeId = sensor.getMillepedeId(); - if (_debug > 0) + 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) + } + + 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) + if (textFile != null) { textFile.printOrigin(origin); + } // associated 3D position of the crossing of this and it's stereo partner sensor - if (textFile != null) + 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())); if (trkpos == null) { - if (_debug > 0) + if (_debug > 0) { System.out.println("Can't find track intercept; use sensor origin"); + } trkpos = strip.origin(); } Hep3Vector trkposTruth = htfTruth != null ? TrackUtils.getHelixPlaneIntercept(htfTruth, strip, Math.abs(_B.z())) : new BasicHep3Vector(-999999.9, -999999.9, -999999.9); - if (textFile != null) + 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"); @@ -287,8 +303,9 @@ 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) + 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 @@ -298,42 +315,39 @@ 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()) + 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) + 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 (_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()); } - } + } + } + + //GBLDATA + GBLStripClusterData stripData = new GBLStripClusterData(millepedeId); + //Add to output list + stripClusterDataList.add(stripData); //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 (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()); - } //GBLDATA stripData.setU(strip.u()); @@ -344,58 +358,59 @@ Hep3Vector tDir = HelixUtils.Direction(htf, s); 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); - } //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())) { - SiTrackerHitStrip1D local = stripHit.getTransformedHit(TrackerHitType.CoordinateSystem.SENSOR); - double d = Math.abs(strip.umeas() - local.getPosition()[0]); - if (d < stripIsoMin && d > 0) - stripIsoMin = d; + if (textFile != null) { + textFile.printStripPathLen(stripData.getPath()); + textFile.printStripPathLen3D(stripData.getPath3D()); + textFile.printMeasDir(stripData.getU()); + textFile.printNonMeasDir(stripData.getV()); + textFile.printNormalDir(stripData.getW()); + textFile.printStripTrackDir(Math.sin(stripData.getTrackPhi()), Math.sin(stripData.getTrackLambda())); + textFile.printStripTrackDirFull(stripData.getTrackDirection()); + } + + if (_debug > 0 || textFile != null) { + // calculate isolation to other strip clusters + double stripIsoMin = 9999.9; + for (SiTrackerHitStrip1D stripHit : stripHits) { + 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) { + stripIsoMin = d; + } + } } - if (_debug > 0) - System.out.printf("%s: stripIsoMin = %f \n", this.getClass().getSimpleName(), stripIsoMin); - - // Add isolation to text file output - if (textFile != null) - textFile.printStripIso(stripIsoMin); + if (_debug > 0) { + System.out.printf("%s: stripIsoMin = %f \n", this.getClass().getSimpleName(), stripIsoMin); + } + + // Add isolation to text file output + 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; // 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; - - // 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()); - - //if(textFile != null) { - // textFile.printStripTrackPosMeasFrame(trkpos_meas); - //} + //GBLDATA stripData.setMeas(strip.umeas()); stripData.setTrackPos(trkpos_meas); + stripData.setMeasErr(strip.du()); if (_debug > 1) { System.out.printf("%s: rotation matrix to meas frame\n%s\n", getClass().getSimpleName(), VecOp.toString(trkToStripRot)); @@ -404,57 +419,65 @@ 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()); - } - - //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()); + textFile.printStripMeas(stripData.getMeas()); + textFile.printStripMeasRes(stripData.getMeas() - stripData.getTrackPos().x(), stripData.getMeasErr()); + } + + if (textFile != null) { + Hep3Vector vdiffTrkTruth = htfTruth != null ? VecOp.sub(trkposTruth, origin) : null; + Hep3Vector trkposTruth_meas = vdiffTrkTruth != null ? VecOp.mult(trkToStripRot, vdiffTrkTruth) : null; + // residual in measurement frame + Double resTruth_meas_x = trkposTruth_meas != null ? strip.umeas() - trkposTruth_meas.x() : null; + textFile.printStripMeasResTruth(resTruth_meas_x != null ? resTruth_meas_x : -9999999.9, strip.du()); + } + + if (_debug > 0) { + System.out.printf("layer %d millePedeId %d uRes %.10f\n", strip.layer(), millepedeId, stripData.getMeas() - stripData.getTrackPos().x()); + } // sim hit residual - if (simHit != null) { + if (isMC && simHit != null) { Hep3Vector simHitPos = CoordinateTransformations.transformVectorToTracking(simHit.getPositionVec()); - if (_debug > 0) + 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()); - } else if (textFile != null) + if (textFile != null) { + textFile.printStripMeasResSimHit(simHitPos_meas.x(), stripData.getMeasErr()); + } + } else if (textFile != null) { textFile.printStripMeasResSimHit(-999999.9, -999999.9); + } // find scattering angle ScatterPoint scatter = scatters.getScatterPoint(((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement()); double scatAngle; - if (scatter != null) + if (scatter != null) { scatAngle = scatter.getScatterAngle().Angle(); - else { - if (_debug > 0) + } 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()); + } scatAngle = GblUtils.estimateScatter(strip, htf, _scattering, _B.magnitude()); } - - //print scatterer to file - if (textFile != null) - textFile.printStripScat(scatAngle); //GBLDATA stripData.setScatterAngle(scatAngle); + //print scatterer to file + if (textFile != null) { + textFile.printStripScat(stripData.getScatterAngle()); + } ++istrip; } } - if (_addBeamspot) + if (_addBeamspot) { addBeamspotToHitList(htf, stripClusterDataList, istrip, htfTruth); + } } @@ -499,7 +522,7 @@ private void addBeamspotToHitList(HelicalTrackFit htf, List<GBLStripClusterData> stripClusterDataList, int istrip, HelicalTrackFit htfTruth) { Hep3Vector targetNormal = new BasicHep3Vector(Math.cos(beamTilt), Math.sin(beamTilt), 0); - int millepedeId = 665; + int millepedeId = 665; List<HelicalTrackStrip> beamspotStrips = makeHelicalTrackStripsFromBeamSpot(); for (HelicalTrackStrip stripOld : beamspotStrips) { HelicalTrackStripGbl strip = new HelicalTrackStripGbl(stripOld, false); @@ -512,13 +535,14 @@ // else // throw new ClassCastException("Detector element " + de.getName() + " couldn't be casted to " + HpsSiSensor.class.getName()); // int millepedeId = sensor.getMillepedeId(); - - millepedeId ++; - if (_debug > 0) + millepedeId++; + if (_debug > 0) { System.out.printf("%s: layer %d millepede %d (DE=\"%s\", origin %s) \n", this.getClass().getSimpleName(), strip.layer(), millepedeId, "BeamSpot", strip.origin().toString()); - - if (textFile != null) + } + + if (textFile != null) { textFile.printStrip(istrip, millepedeId, "BeamSpot"); + } //GBLDATA GBLStripClusterData stripData = new GBLStripClusterData(millepedeId); @@ -528,24 +552,28 @@ //Center of the sensor Hep3Vector origin = strip.origin(); - if (textFile != null) + if (textFile != null) { textFile.printOrigin(origin); + } // associated 3D position of beamspot on target - if (textFile != null) + if (textFile != null) { textFile.printHitPos3D(new BasicHep3Vector(0, 0, 0)); + } //Find intercept point with sensor in tracking frame // Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(htf, strip, Math.abs(_B.z())); Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(htf, strip.w(), strip.origin(), Math.abs(_B.z())); if (trkpos == null) { - if (_debug > 0) + if (_debug > 0) { System.out.println("Can't find track intercept; use sensor origin"); + } trkpos = strip.origin(); } Hep3Vector trkposTruth = htfTruth != null ? TrackUtils.getHelixPlaneIntercept(htfTruth, strip.w(), strip.origin(), Math.abs(_B.z())) : new BasicHep3Vector(-999999.9, -999999.9, -999999.9); - if (textFile != null) + 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"); @@ -559,8 +587,9 @@ 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) + 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 @@ -635,12 +664,14 @@ // stripIsoMin = d; // } - if (_debug > 0) + 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 @@ -649,9 +680,9 @@ // then find the rotation from tracking to measurement frame // Hep3Matrix trkToStripRot = _trackerHitUtils.getTrackToStripRotation(strip.getStrip()); - Hep3Matrix trkToBeamSpot = new BasicHep3Matrix(Math.cos(beamTilt),Math.sin(beamTilt),0, - -Math.sin(beamTilt),Math.cos(beamTilt),0, - 0,0,1); + Hep3Matrix trkToBeamSpot = new BasicHep3Matrix(Math.cos(beamTilt), Math.sin(beamTilt), 0, + -Math.sin(beamTilt), Math.cos(beamTilt), 0, + 0, 0, 1); // then rotate that vector into the measurement frame to get the predicted measurement position Hep3Vector trkpos_meas = VecOp.mult(trkToBeamSpot, vdiffTrk); Hep3Vector trkposTruth_meas = vdiffTrkTruth != null ? VecOp.mult(trkToBeamSpot, vdiffTrkTruth) : null; @@ -660,8 +691,9 @@ 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) + if (textFile != null) { textFile.printStripMeas(m_meas.x()); + } //if(textFile != null) { // textFile.printStripTrackPosMeasFrame(trkpos_meas); @@ -688,15 +720,17 @@ //GBLDATA stripData.setMeasErr(res_err_meas.x()); - if (_debug > 0) + if (_debug > 0) { System.out.printf("layer %d millePedeId %d uRes %.10f\n", strip.layer(), millepedeId, res_meas.x()); + } // double scatAngle = 0.0001; //print scatterer to file - if (textFile != null) + if (textFile != null) { textFile.printStripScat(scatAngle); + } //GBLDATA stripData.setScatterAngle(scatAngle); @@ -718,16 +752,19 @@ // 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) + 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) + } else if (_debug > 0) { System.out.printf("%s: the matched MC particle is the recoil electron\n", this.getClass().getSimpleName()); + } } } @@ -736,19 +773,22 @@ Map<MCParticle, Integer> particlesOnTrack = new HashMap<MCParticle, Integer>(); - if (debug) + 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) + 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) + } 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)) + if (!particlesOnTrack.containsKey(mcp)) { particlesOnTrack.put(mcp, 0); + } int c = particlesOnTrack.get(mcp); particlesOnTrack.put(mcp, c + 1); } @@ -757,23 +797,26 @@ 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()) + 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) + 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; - //} + } // if(entry.getValue().compareTo(maxEntry.getValue()) < 0) continue; + } //} //maxEntry = entry; - if (debug) - if (maxEntry != null) + 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 + } else { System.out.printf("No truth particle found on this track\n"); + } + } return maxEntry == null ? null : maxEntry.getKey(); } @@ -906,22 +949,27 @@ 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) + 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) + if (mcp.getGeneratorStatus() != MCParticle.FINAL_STATE) { continue; + } boolean hasAprimeParent = false; - for (MCParticle parent : mcp.getParents()) - if (Math.abs(parent.getPDGID()) == 622) + for (MCParticle parent : mcp.getParents()) { + if (Math.abs(parent.getPDGID()) == 622) { hasAprimeParent = true; - if (hasAprimeParent) + } + } + 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()); @@ -946,8 +994,9 @@ 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) + 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:" : ""); @@ -955,10 +1004,11 @@ 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()) + 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:" : ""); + } } } @@ -977,19 +1027,23 @@ return Math.sqrt(Math.pow(E1 + E2, 2) - VecOp.add(p1vec, p2vec).magnitudeSquared()); } + private static BasicMatrix getPerParVector(double kappa, double theta, double phi, double d0, double z0) { + BasicMatrix perPar = new BasicMatrix(1, 5); + perPar.setElement(0, 0, kappa); + perPar.setElement(0, 1, theta); + perPar.setElement(0, 2, phi); + perPar.setElement(0, 3, d0); + perPar.setElement(0, 4, z0); + return perPar; + } + private static BasicMatrix getPerParVector(HelicalTrackFit htf, double B) { - BasicMatrix perPar = new BasicMatrix(1, 5); if (htf != null) { double kappa = -1.0 * Math.signum(B) / htf.R(); 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; - + return getPerParVector(kappa, theta, htf.phi0(), htf.dca(), htf.z0()); + } + return new BasicMatrix(1, 5); } public static class PerigeeParams { @@ -998,6 +1052,10 @@ public PerigeeParams(HelicalTrackFit htf, double B) { _params = getPerParVector(htf, B); + } + + public PerigeeParams(double kappa, double theta, double phi, double d0, double z0) { + this._params = getPerParVector(kappa, theta, phi, d0, z0); } public BasicMatrix getParams() { @@ -1082,8 +1140,9 @@ public ClParams(HelicalTrackFit htf, double B) { - if (htf == null) + if (htf == null) { return; + } Hep3Matrix perToClPrj = getPerToClPrj(htf); Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java ============================================================================= --- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java (original) +++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java Wed Sep 30 18:41:59 2015 @@ -7,143 +7,146 @@ import org.lcsim.event.GenericObject; /** - * - * + * + * * @version $Id: */ public class GBLTrackData implements GenericObject { - - /* - * - * Interface enumerator to access the correct data - * - */ - private static class GBLINT { - public static final int ID = 0; - public static final int BANK_INT_SIZE = 1; - } - public static class GBLDOUBLE { - public static final int PERKAPPA =0; - public static final int PERTHETA = 1; - public static final int PERPHI = 2; - public static final int PERD0 = 3; - public static final int PERZ0 = 4; - // 9 entries from projection matrix from perigee to curvilinear frame - public static final int BANK_DOUBLE_SIZE = 5+9; - } - // array holding the integer data - private int bank_int[] = new int[GBLINT.BANK_INT_SIZE]; - // array holding the double data - private double bank_double[] = new double[GBLDOUBLE.BANK_DOUBLE_SIZE]; - - /** - * Default constructor - */ - public GBLTrackData(int id) { - setTrackId(id); - } - - /* - * Constructor from GenericObject - * TODO add size checks for backwards compatability - */ - public GBLTrackData(GenericObject o) - { - for(int i=0; i<GBLINT.BANK_INT_SIZE; ++i) - { - bank_int[i] = o.getIntVal(i); - } - for(int i=0; i<GBLDOUBLE.BANK_DOUBLE_SIZE; ++i) - { - bank_double[i] = o.getDoubleVal(i); - } - - } - - /** - * @param set track id to val - */ - public void setTrackId(int val) { - bank_int[GBLINT.ID] = val; - } - - /** - * @return track id for this object - */ - public int getTrackId() { - return this.getIntVal(GBLINT.ID); - } - - /** - * @param perPar is the perigee parameters that is added to object - */ - public void setPerigeeTrackParameters(PerigeeParams perPar) { - this.bank_double[GBLDOUBLE.PERKAPPA] = perPar.getKappa(); - this.bank_double[GBLDOUBLE.PERTHETA] = perPar.getTheta(); - this.bank_double[GBLDOUBLE.PERPHI] = perPar.getPhi(); - this.bank_double[GBLDOUBLE.PERD0] = perPar.getD0(); - this.bank_double[GBLDOUBLE.PERZ0] = perPar.getZ0(); - } - - - - public void setPrjPerToCl(int row, int col, double val) { - int idx = col + row*3; - if(idx>8) { - System.out.printf("%s: ERROR to large matrix\n", this.getClass().getSimpleName()); - System.exit(1); - } - this.bank_double[idx+5] = val; - } - public Hep3Matrix getPrjPerToCl() { - BasicHep3Matrix matrix = new BasicHep3Matrix(); - for(int row=0; row<3; ++row) { - for(int col=0; col<3; ++ col) { - matrix.setElement(row, col, getPrjPerToClVal(row, col)); - } - } - return matrix; - } + /* + * + * Interface enumerator to access the correct data + * + */ + private static class GBLINT { - private double getPrjPerToClVal(int row, int col) { - int idx = col + row*3; - return this.bank_double[idx+5]; + public static final int ID = 0; + public static final int BANK_INT_SIZE = 1; } - /* - * The functions below are all overide from - * @see org.lcsim.event.GenericObject#getNInt() - */ - - public int getNInt() { - return GBLINT.BANK_INT_SIZE; - } + public static class GBLDOUBLE { - public int getNFloat() { - return 0; - } + public static final int PERKAPPA = 0; + public static final int PERTHETA = 1; + public static final int PERPHI = 2; + public static final int PERD0 = 3; + public static final int PERZ0 = 4; + // 9 entries from projection matrix from perigee to curvilinear frame + public static final int BANK_DOUBLE_SIZE = 5 + 9; + } + // array holding the integer data + private int bank_int[] = new int[GBLINT.BANK_INT_SIZE]; + // array holding the double data + private double bank_double[] = new double[GBLDOUBLE.BANK_DOUBLE_SIZE]; - public int getNDouble() { - return GBLDOUBLE.BANK_DOUBLE_SIZE; - } + /** + * Default constructor + */ + public GBLTrackData(int id) { + setTrackId(id); + } - public int getIntVal(int index) { - return bank_int[index]; - } + /* + * Constructor from GenericObject + * TODO add size checks for backwards compatability + */ + public GBLTrackData(GenericObject o) { + for (int i = 0; i < GBLINT.BANK_INT_SIZE; ++i) { + bank_int[i] = o.getIntVal(i); + } + for (int i = 0; i < GBLDOUBLE.BANK_DOUBLE_SIZE; ++i) { + bank_double[i] = o.getDoubleVal(i); + } - public float getFloatVal(int index) { - return 0; - } + } - public double getDoubleVal(int index) { - return bank_double[index]; - } + /** + * @param set track id to val + */ + public void setTrackId(int val) { + bank_int[GBLINT.ID] = val; + } - public boolean isFixedSize() { - return false; - } + /** + * @return track id for this object + */ + public int getTrackId() { + return this.getIntVal(GBLINT.ID); + } - + /** + * @param perPar is the perigee parameters that is added to object + */ + public void setPerigeeTrackParameters(PerigeeParams perPar) { + this.bank_double[GBLDOUBLE.PERKAPPA] = perPar.getKappa(); + this.bank_double[GBLDOUBLE.PERTHETA] = perPar.getTheta(); + this.bank_double[GBLDOUBLE.PERPHI] = perPar.getPhi(); + this.bank_double[GBLDOUBLE.PERD0] = perPar.getD0(); + this.bank_double[GBLDOUBLE.PERZ0] = perPar.getZ0(); + } -} + public PerigeeParams getPerigeeTrackParameters() { + return new PerigeeParams(this.bank_double[GBLDOUBLE.PERKAPPA], + this.bank_double[GBLDOUBLE.PERTHETA], + this.bank_double[GBLDOUBLE.PERPHI], + this.bank_double[GBLDOUBLE.PERD0], + this.bank_double[GBLDOUBLE.PERZ0]); + } + + public void setPrjPerToCl(int row, int col, double val) { + int idx = col + row * 3; + if (idx > 8) { + System.out.printf("%s: ERROR to large matrix\n", this.getClass().getSimpleName()); + System.exit(1); + } + this.bank_double[idx + 5] = val; + } + + public Hep3Matrix getPrjPerToCl() { + BasicHep3Matrix matrix = new BasicHep3Matrix(); + for (int row = 0; row < 3; ++row) { + for (int col = 0; col < 3; ++col) { + matrix.setElement(row, col, getPrjPerToClVal(row, col)); + } + } + return matrix; + } + + private double getPrjPerToClVal(int row, int col) { + int idx = col + row * 3; + return this.bank_double[idx + 5]; + } + + /* + * The functions below are all overide from + * @see org.lcsim.event.GenericObject#getNInt() + */ + public int getNInt() { + return GBLINT.BANK_INT_SIZE; + } + + public int getNFloat() { + return 0; + } + + public int getNDouble() { + return GBLDOUBLE.BANK_DOUBLE_SIZE; + } + + public int getIntVal(int index) { + return bank_int[index]; + } + + public float getFloatVal(int index) { + return 0; + } + + public double getDoubleVal(int index) { + return bank_double[index]; + } + + public boolean isFixedSize() { + return false; + } + +} Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java ============================================================================= --- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java (original) +++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java Wed Sep 30 18:41:59 2015 @@ -158,12 +158,16 @@ int trackNum = 0; logger.info("Trying to fit " + stripsGblMap.size() + " tracks"); for (GBLTrackData t : stripsGblMap.keySet()) { - FittedGblTrajectory traj = fit(stripsGblMap.get(t), bfac); + FittedGblTrajectory traj = fit(stripsGblMap.get(t), bfac, _debug); ++trackNum; if (traj != null) { logger.info("GBL fit successful"); if (_debug) { System.out.printf("%s: GBL fit successful.\n", getClass().getSimpleName()); + } + // write to MP binary file + if (writeMilleBinary) { + traj.get_traj().milleOut(mille); } traj.set_seed(gblToSeedMap.get(t)); traj.set_track_data(t); @@ -190,7 +194,7 @@ } - private FittedGblTrajectory fit(List<GBLStripClusterData> hits, double bfac) { + private static FittedGblTrajectory fit(List<GBLStripClusterData> hits, double bfac, boolean debug) { // path length along trajectory double s = 0.; @@ -213,14 +217,15 @@ for (int istrip = 0; istrip != n_strips; ++istrip) { GBLStripClusterData strip = hits.get(istrip); //MG--9/18/2015--beamspot has Id=666/667...don't include it in the GBL fit - if(strip.getId()>99) + if (strip.getId() > 99) { continue; - if (_debug) { + } + if (debug) { System.out.println("HpsGblFitter: Processing strip " + istrip + " with id/layer " + strip.getId()); } // Path length step for this cluster double step = strip.getPath3D() - s; - if (_debug) { + if (debug) { System.out.println("HpsGblFitter: " + "Path length step " + step + " from " + s + " to " + strip.getPath3D()); } @@ -237,13 +242,13 @@ mDir.set(1, 0, v.x()); mDir.set(1, 1, v.y()); mDir.set(1, 2, v.z()); - if (_debug) { + if (debug) { System.out.println("HpsGblFitter: " + "mDir"); mDir.print(4, 6); } Matrix mDirT = mDir.copy().transpose(); - if (_debug) { + if (debug) { System.out.println("HpsGblFitter: " + "mDirT"); mDirT.print(4, 6); } @@ -254,7 +259,7 @@ double sinPhi = sin(strip.getTrackPhi());//->GetPhi()); double cosPhi = sqrt(1.0 - sinPhi * sinPhi); - if (_debug) { + if (debug) { System.out.println("HpsGblFitter: " + "Track direction sinLambda=" + sinLambda + " sinPhi=" + sinPhi); } @@ -268,7 +273,7 @@ uvDir.set(1, 1, -sinLambda * sinPhi); uvDir.set(1, 2, cosLambda); - if (_debug) { + if (debug) { System.out.println("HpsGblFitter: " + "uvDir"); uvDir.print(6, 4); } @@ -281,7 +286,7 @@ proL2m = proL2m.inverse(); proL2m_list.put(strip.getId(), proL2m.copy()); // is a copy needed or is that just a C++/root thing? - if (_debug) { + if (debug) { System.out.println("HpsGblFitter: " + "proM2l:"); proM2l.print(4, 6); System.out.println("HpsGblFitter: " + "proL2m:"); @@ -305,7 +310,7 @@ measPrec.set(0, 1.0 / (measErr.get(0) * measErr.get(0))); measPrec.set(1, 0.); - if (_debug) { + if (debug) { System.out.println("HpsGblFitter: " + "meas: "); meas.print(4, 6); System.out.println("HpsGblFitter: " + "measErr:"); @@ -317,7 +322,7 @@ //Find the Jacobian to be able to propagate the covariance matrix to this strip position jacPointToPoint = gblSimpleJacobianLambdaPhi(step, cosLambda, abs(bfac)); - if (_debug) { + if (debug) { System.out.println("HpsGblFitter: " + "jacPointToPoint to extrapolate to this point:"); jacPointToPoint.print(4, 6); } @@ -341,7 +346,7 @@ Matrix measMsCov = proL2m.times((msCov.getMatrix(3, 4, 3, 4)).times(proL2mTransposed)); - if (_debug) { + if (debug) { System.out.println("HpsGblFitter: " + " msCov at this point:"); msCov.print(4, 6); System.out.println("HpsGblFitter: " + "measMsCov at this point:"); @@ -367,7 +372,7 @@ // add scatterer if not using the uncorrelated MS covariances for testing if (!useUncorrMS) { point.addScatterer(scat, scatPrec); - if (_debug) { + if (debug) { System.out.println("HpsGblFitter: " + "adding scatError to this point:"); scatErr.print(4, 6); } @@ -440,7 +445,7 @@ } // print the trajectory - if (_debug) { + if (debug) { System.out.println("%%%% Gbl Trajectory "); traj.printTrajectory(1); traj.printData(); @@ -454,7 +459,7 @@ Vector aCorrection = new Vector(5); SymMatrix aCovariance = new SymMatrix(5); traj.getResults(1, aCorrection, aCovariance); - if (_debug) { + if (debug) { System.out.println(" cor "); aCorrection.print(6, 4); System.out.println(" cov "); @@ -463,10 +468,6 @@ logger.fine("locPar " + aCorrection.toString()); -// // write to MP binary file - if (writeMilleBinary) { - traj.milleOut(mille); - } // return new FittedGblTrajectory(traj, dVals[0], iVals[0], dVals[1]); } @@ -475,7 +476,7 @@ protected void detectorChanged(Detector detector) { } - private Matrix gblSimpleJacobianLambdaPhi(double ds, double cosl, double bfac) { + private static Matrix gblSimpleJacobianLambdaPhi(double ds, double cosl, double bfac) { /** * Simple jacobian: quadratic in arc length difference. using lambda phi * as directions