Author: [log in to unmask]
Date: Tue Jan 26 18:37:32 2016
New Revision: 4143
Log:
Merge trunk changes.
Added:
java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v3-5-3-fieldmap/
- copied from r4139, java/trunk/detector-data/detectors/HPS-EngRun2015-Nominal-v3-5-3-fieldmap/
java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/SamplingFractions/
- copied from r4139, java/trunk/detector-data/detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/SamplingFractions/
java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v5-fieldmap/
- copied from r4139, java/trunk/detector-data/detectors/HPS-EngRun2015-Nominal-v5-fieldmap/
java/branches/jeremy-dev/recon/src/main/java/org/hps/recon/filtering/SvtRawHitMultiplicityFilter.java
- copied unchanged from r4139, java/trunk/recon/src/main/java/org/hps/recon/filtering/SvtRawHitMultiplicityFilter.java
java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/users/phansson/EventInfo.lcsim
- copied unchanged from r4139, java/trunk/steering-files/src/main/resources/org/hps/steering/users/phansson/EventInfo.lcsim
java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/users/phansson/HPSReconNoReadout.lcsim
- copied unchanged from r4139, java/trunk/steering-files/src/main/resources/org/hps/steering/users/phansson/HPSReconNoReadout.lcsim
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/HpsHelicalTrackFit.java
- copied unchanged from r4139, java/trunk/tracking/src/main/java/org/hps/recon/tracking/HpsHelicalTrackFit.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/byale/
- copied from r4139, java/trunk/users/src/main/java/org/hps/users/byale/
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/GblResidualDriver.java
- copied unchanged from r4139, java/trunk/users/src/main/java/org/hps/users/phansson/GblResidualDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/PrintEventInfoDriver.java
- copied unchanged from r4139, java/trunk/users/src/main/java/org/hps/users/phansson/PrintEventInfoDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/alignment/
- copied from r4139, java/trunk/users/src/main/java/org/hps/users/phansson/alignment/
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/apps/
- copied from r4139, java/trunk/users/src/main/java/org/hps/users/phansson/apps/
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/daq/
- copied from r4139, java/trunk/users/src/main/java/org/hps/users/phansson/daq/
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/gbl/
- copied from r4139, java/trunk/users/src/main/java/org/hps/users/phansson/gbl/
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/testrun/
- copied from r4139, java/trunk/users/src/main/java/org/hps/users/phansson/testrun/
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/tools/
- copied from r4139, java/trunk/users/src/main/java/org/hps/users/phansson/tools/
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/ExtractFormFactors.java
- copied unchanged from r4139, java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ExtractFormFactors.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/FEESpectrumGenerator.java
- copied unchanged from r4139, java/trunk/users/src/main/java/org/hps/users/spaul/feecc/FEESpectrumGenerator.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/FormFactor.java
- copied unchanged from r4139, java/trunk/users/src/main/java/org/hps/users/spaul/feecc/FormFactor.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/MottIntegral.java
- copied unchanged from r4139, java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MottIntegral.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/MultipleScattering.java
- copied unchanged from r4139, java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MultipleScattering.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/RemoveDuplicateParticles.java
- copied unchanged from r4139, java/trunk/users/src/main/java/org/hps/users/spaul/feecc/RemoveDuplicateParticles.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinningXY.java
- copied unchanged from r4139, java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinningXY.java
Removed:
java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/Ecal.properties
java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/production/DataQualityRecon_Pass2.lcsim
java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/recon/EngineeringRun2015FullReconMC_Pass2.lcsim
java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/recon/EngineeringRun2015FullRecon_Pass2.lcsim
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/HPSTrack.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLDriver.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblFitter.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/AlignmentUtils.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/BeamCurrentData.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/Count.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/DAQDeadTimeData.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/DumpAIDATextFiles.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/ECalExtrapolationDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/EcalHitMapPlots.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/GlobalParameter.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/GlobalParameters.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/MPAlignmentInputCalculator.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/ModuleMPAlignmentInput.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/PolarCount.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/ROOTFlatTupleDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/ReadSurveyRotations.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/RegExpMatcherTester.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/ResLimit.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/RunMPAlignment.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/SimpleHPSConditions.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/StripMPAlignmentInput.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/SvtHeaderAnalysisDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/SvtHeaderMetaDataReaderDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/SvtOldHeaderAnalysisDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/SvtOldHeaderDataInfo.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TestSort.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TrackingGeometryChecker.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TrigRateAna.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TrigRateDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/WTrack.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/dataMCPlots.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/ecalPlots.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/mergeSimpleAIDA.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/trigRate.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/analysis/
Modified:
java/branches/jeremy-dev/ (props changed)
java/branches/jeremy-dev/analysis/src/main/java/org/hps/analysis/dataquality/TridentMonitoring.java
java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/HPS-EngRun2015-Nominal-v4-4-fieldmap.lcdd
java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/compact.xml
java/branches/jeremy-dev/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterEnergyCorrection.java
java/branches/jeremy-dev/evio/src/main/java/org/hps/evio/RfFitterDriver.java
java/branches/jeremy-dev/monitoring-app/ (props changed)
java/branches/jeremy-dev/monitoring-app/src/main/java/org/hps/monitoring/application/EventProcessing.java
java/branches/jeremy-dev/monitoring-drivers/src/main/java/org/hps/monitoring/drivers/trackrecon/TrackingReconPlots.java
java/branches/jeremy-dev/recon/src/main/java/org/hps/recon/vertexing/TwoTrackFringeVertexer.java
java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/production/FilterHitsPerSensor.lcsim
java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/recon/EngineeringRun2014EcalReconMC.lcsim
java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/recon/EngineeringRun2015FullReconMC.lcsim
java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/users/phansson/EngineeringRun2015FullRecon_Pass2_Gbl.lcsim
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/WTrack.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLFileIO.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GblTrajectory.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java
java/branches/jeremy-dev/tracking/src/main/java/org/hps/svt/alignment/BuildMillepedeCompact.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/meeg/SvtChargeIntegrator.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/ExamplePlotter.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/HelicalTrackHitResidualsDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/TrackExtrapolationAnalysis.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/TwoTrackAnalysis.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/CmpGenToFittedTracksDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/DataTrackerFakeHitDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/ParticleHelixProducer.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TrackingReconstructionPlots.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TriggerTurnOnAnalysis.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TruthMomentumResolutionDriver.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/StyleUtil.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/SumEverything.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/CustomBinning.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/DisplayHistograms.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/EcalUtil.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/MakeHistograms.java
java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java
Modified: java/branches/jeremy-dev/analysis/src/main/java/org/hps/analysis/dataquality/TridentMonitoring.java
=============================================================================
--- java/branches/jeremy-dev/analysis/src/main/java/org/hps/analysis/dataquality/TridentMonitoring.java (original)
+++ java/branches/jeremy-dev/analysis/src/main/java/org/hps/analysis/dataquality/TridentMonitoring.java Tue Jan 26 18:37:32 2016
@@ -9,17 +9,18 @@
import hep.physics.vec.BasicHep3Matrix;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
-
import java.util.ArrayList;
import java.util.EnumSet;
import java.util.List;
import java.util.Map.Entry;
import java.util.logging.Logger;
-
import org.hps.recon.ecal.cluster.ClusterUtilities;
import org.hps.recon.particle.ReconParticleDriver;
import org.hps.recon.tracking.TrackType;
import org.hps.recon.tracking.TrackUtils;
+import org.hps.recon.vertexing.BilliorTrack;
+import org.hps.recon.vertexing.BilliorVertex;
+import org.hps.recon.vertexing.BilliorVertexer;
import org.lcsim.event.EventHeader;
import org.lcsim.event.ReconstructedParticle;
import org.lcsim.event.RelationalTable;
@@ -155,10 +156,10 @@
private IHistogram1D vertMass;
private IHistogram2D vertZVsMass;
// private IHistogram1D vertX;
-// private IHistogram1D vertY;
+ private IHistogram1D vertY;
// private IHistogram1D vertZ;
-// private IHistogram2D vertZY;
-// private IHistogram2D vertXY;
+ private IHistogram2D vertZY;
+ private IHistogram2D vertXY;
// private IHistogram1D vertPx;
// private IHistogram1D vertPy;
// private IHistogram1D vertPz;
@@ -189,6 +190,8 @@
private IHistogram1D vertRadU;
private IHistogram1D vertRadV;
+ private IHistogram2D vertRadUnconBsconChi2;
+
private IHistogram1D nTriCand;
private IHistogram1D nVtxCand;
// IHistogram1D vertexW;
@@ -224,15 +227,25 @@
private final double maxChi2GBLTrack = 15.0;
private final double maxVertChi2 = 7.0;
- //v0 cuts
+ //v0 plot ranges
private final double v0PzMax = 1.25 * ebeam;//GeV
private final double v0PzMin = 0.1;// GeV
private final double v0PyMax = 0.04;//GeV absolute value
private final double v0PxMax = 0.04;//GeV absolute value
private final double v0VzMax = 50.0;// mm from target...someday make mass dependent
- private final double v0VyMax = 1.0;// mm from target...someday make mass dependent
+ private final double v0VyMax = 2.0;// mm from target...someday make mass dependent
private final double v0VxMax = 2.0;// mm from target...someday make mass dependent
- // track quality cuts
+
+ //v0 cuts
+ private final double v0PzMaxCut = 1.25 * ebeam;//GeV
+ private final double v0PzMinCut = 0.1;// GeV
+ private final double v0PyCut = 0.04;//GeV absolute value
+ private final double v0PxCut = 0.04;//GeV absolute value
+ private final double v0VzCut = 50.0;// mm from target...someday make mass dependent
+ private final double v0VyCut = 2.0;// mm from target...someday make mass dependent
+ private final double v0VxCut = 2.0;// mm from target...someday make mass dependent
+
+// track quality cuts
private final double beamPCut = 0.85;
private final double minPCut = 0.05;
// private double trkPyMax = 0.2;
@@ -242,12 +255,14 @@
private final double clusterTimeDiffCut = 2.5;
private final double l1IsoMin = 1.0;
+
+ private double[] beamSize = {0.001, 0.130, 0.050}; //rough estimate from harp scans during engineering run production running
+
//cluster matching
// private boolean reqCluster = false;
// private int nClustMax = 3;
// private double eneLossFactor = 0.7; //average E/p roughly
// private double eneOverPCut = 0.3; //|(E/p)_meas - (E/p)_mean|<eneOverPCut
-
//counters
private float nEvents = 0;
private float nRecoV0 = 0;
@@ -346,10 +361,10 @@
vertMass = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex mass", 100, 0, 0.11);
vertZVsMass = aida.histogram2D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex Z vs. mass", 100, 0, 0.11, 100, -v0VzMax, v0VzMax);
// vertX = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex X", 100, -v0VxMax, v0VxMax);
-// vertY = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex Y", 100, -v0VyMax, v0VyMax);
+ vertY = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex Y", 100, -v0VyMax, v0VyMax);
// vertZ = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex Z", 100, -v0VzMax, v0VzMax);
-// vertXY = aida.histogram2D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex Y vs. X", 100, -v0VxMax, v0VxMax, 100, -v0VyMax, v0VyMax);
-// vertZY = aida.histogram2D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex Z vs. Y", 100, -v0VyMax, v0VyMax, 100, -v0VzMax, v0VzMax);
+ vertXY = aida.histogram2D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex Y vs. X", 100, -v0VxMax, v0VxMax, 100, -v0VyMax, v0VyMax);
+ vertZY = aida.histogram2D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex Z vs. Y", 100, -v0VyMax, v0VyMax, 100, -v0VzMax, v0VzMax);
// vertPx = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex Px", 100, -v0PxMax, v0PxMax);
// vertPy = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex Py", 100, -v0PyMax, v0PyMax);
// vertPz = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Vertex: Vertex Pz", 100, v0PzMin, v0PzMax);
@@ -382,6 +397,8 @@
vertRadU = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Radiative vertex: Vertex Px over Ptot", 100, -0.1, 0.1);
vertRadV = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Radiative vertex: Vertex Py over Ptot", 100, -0.1, 0.1);
+ vertRadUnconBsconChi2 = aida.histogram2D(plotDir + trkType + triggerType + "/" + "Radiative vertex: beamspot chi2 vs. uncon chi2", 100, 0, 25.0, 100, 0, 25.0);
+
nVtxCand = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Number of Vertexing Candidates", 5, 0, 4);
maxTrkChi2 = aida.histogram1D(plotDir + trkType + triggerType + "/" + "Cut: Trk Chi2", 50, 0.0, 50.0);
@@ -527,8 +544,8 @@
bits.add(Cut.VTX_QUALITY);
}
- boolean vertexMomentumCut = v0MomRot.z() < v0PzMax && v0MomRot.z() > v0PzMin && Math.abs(v0MomRot.x()) < v0PxMax && Math.abs(v0MomRot.y()) < v0PyMax;
- boolean vertexPositionCut = Math.abs(v0Vtx.x()) < v0VxMax && Math.abs(v0Vtx.y()) < v0VyMax && Math.abs(v0Vtx.z()) < v0VzMax;
+ boolean vertexMomentumCut = v0MomRot.z() < v0PzMaxCut && v0MomRot.z() > v0PzMinCut && Math.abs(v0MomRot.x()) < v0PxCut && Math.abs(v0MomRot.y()) < v0PyCut;
+ boolean vertexPositionCut = Math.abs(v0Vtx.x()) < v0VxCut && Math.abs(v0Vtx.y()) < v0VyCut && Math.abs(v0Vtx.z()) < v0VzCut;
if (vertexMomentumCut && vertexPositionCut) {
bits.add(Cut.VERTEX_CUTS);
}
@@ -715,6 +732,7 @@
if (!vertCandidateList.isEmpty()) {
// pick the best candidate...for now just pick a random one.
ReconstructedParticle bestCandidate = vertCandidateList.get((int) (Math.random() * vertCandidateList.size()));
+ Vertex unconVertex = bestCandidate.getStartVertex();
//fill some stuff:
ReconstructedParticle electron = bestCandidate.getParticles().get(ReconParticleDriver.ELECTRON);
@@ -728,7 +746,7 @@
Hep3Vector pBestV0Rot = VecOp.mult(beamAxisRotation, bestCandidate.getMomentum());
Hep3Vector pEleRot = VecOp.mult(beamAxisRotation, electron.getMomentum());
Hep3Vector pPosRot = VecOp.mult(beamAxisRotation, positron.getMomentum());
- Hep3Vector v0Vtx = VecOp.mult(beamAxisRotation, bestCandidate.getStartVertex().getPosition());
+ Hep3Vector v0Vtx = VecOp.mult(beamAxisRotation, unconVertex.getPosition());
// vertTrackTime2D.fill(tEle, tPos);
// vertTrackTimeDiff.fill(tEle - tPos);
@@ -744,16 +762,26 @@
vertMass.fill(bestCandidate.getMass());
vertZVsMass.fill(bestCandidate.getMass(), v0Vtx.z());
// vertX.fill(v0Vtx.x());
-// vertY.fill(v0Vtx.y());
+ vertY.fill(v0Vtx.y());
// vertZ.fill(v0Vtx.z());
// vertPx.fill(pBestV0Rot.x());
// vertPy.fill(pBestV0Rot.y());
// vertPz.fill(pBestV0Rot.z());
// vertU.fill(pBestV0Rot.x() / pBestV0Rot.magnitude());
// vertV.fill(pBestV0Rot.y() / pBestV0Rot.magnitude());
-// vertXY.fill(v0Vtx.x(), v0Vtx.y());
-// vertZY.fill(v0Vtx.y(), v0Vtx.z());
+ vertXY.fill(v0Vtx.x(), v0Vtx.y());
+ vertZY.fill(v0Vtx.y(), v0Vtx.z());
if (bestCandidate.getMomentum().magnitude() > radCut) {
+
+ BilliorVertexer vtxFitter = new BilliorVertexer(TrackUtils.getBField(event.getDetector()).y());
+ vtxFitter.setBeamSize(beamSize);
+ List<BilliorTrack> billiorTracks = new ArrayList<BilliorTrack>();
+ billiorTracks.add(new BilliorTrack(electron.getTracks().get(0)));
+ billiorTracks.add(new BilliorTrack(positron.getTracks().get(0)));
+ vtxFitter.doBeamSpotConstraint(true);
+ BilliorVertex bsconVertex = vtxFitter.fitVertex(billiorTracks);
+ vertRadUnconBsconChi2.fill(unconVertex.getChi2(), bsconVertex.getChi2());
+
vertRadTrackTime2D.fill(tEle, tPos);
vertRadTrackTimeDiff.fill(tEle - tPos);
vertRadZVsMomentum.fill(bestCandidate.getMomentum().magnitude(), v0Vtx.z());
Modified: java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/HPS-EngRun2015-Nominal-v4-4-fieldmap.lcdd
=============================================================================
--- java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/HPS-EngRun2015-Nominal-v4-4-fieldmap.lcdd (original)
+++ java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/HPS-EngRun2015-Nominal-v4-4-fieldmap.lcdd Tue Jan 26 18:37:32 2016
@@ -4,7 +4,7 @@
<detector name="HPS-EngRun2015-Nominal-v4-4-fieldmap" />
<generator name="GeomConverter" version="1.0" file="detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/compact.xml" checksum="546271450" />
<author name="NONE" />
- <comment>HPS detector for 2015 Engineering Run with tracker fully open this detector uses the corrected fieldmap scaled to -0.24T for 1.056 GeV running including SVT mechanical survey. The ECal top is shifted up 2.7 mm from nominal and the bottom is shifted down 0.9 mm from nominal. SVT alignment: 1) start from v3-2-fieldmap 2) beamspot BS_D: 3)BS_D-20k.gbl-5mradScatAngle-L01_L12_L23_L01_L12_L34_L45_L01_L12_tu</comment>
+ <comment>HPS detector for 2015 Engineering Run with tracker fully open this detector uses the corrected fieldmap scaled to -0.24T for 1.056 GeV running including SVT mechanical survey. The ECal top is shifted up 2.3 mm from nominal and the bottom is shifted down 2.4 mm from nominal, based only on tracking. SVT alignment: 1) start from v3-2-fieldmap 2) beamspot BS_D: 3)BS_D-20k.gbl-5mradScatAngle-L01_L12_L23_L01_L12_L34_L45_L01_L12_tu</comment>
</header>
<iddict>
<idspec name="TrackerHitsFieldDef" length="46">
Modified: java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/compact.xml
=============================================================================
--- java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/compact.xml (original)
+++ java/branches/jeremy-dev/detector-data/detectors/HPS-EngRun2015-Nominal-v4-4-fieldmap/compact.xml Tue Jan 26 18:37:32 2016
@@ -4,8 +4,8 @@
<info name="HPS-EngRun2015-Nominal-v4-4-fieldmap">
<comment>HPS detector for 2015 Engineering Run with tracker fully open
this detector uses the corrected fieldmap scaled to -0.24T for 1.056 GeV
- running including SVT mechanical survey. The ECal top is shifted up 2.7 mm
- from nominal and the bottom is shifted down 0.9 mm from nominal.
+ running including SVT mechanical survey. The ECal top is shifted up 2.3 mm
+ from nominal and the bottom is shifted down 2.4 mm from nominal, based only on tracking.
SVT alignment:
1) start from v3-2-fieldmap
2) beamspot BS_D:
Modified: java/branches/jeremy-dev/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterEnergyCorrection.java
=============================================================================
--- java/branches/jeremy-dev/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterEnergyCorrection.java (original)
+++ java/branches/jeremy-dev/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterEnergyCorrection.java Tue Jan 26 18:37:32 2016
@@ -12,9 +12,8 @@
/**
* This is the cluster energy correction requiring the particle id
- * uncorrected cluster energy. This is not accurate for edge crystals
- * and should be used after cluster-track matching and after position
- * corrections.
+ * uncorrected cluster energy. This is now updated to include edge
+ * corrections and sampling fractions derived from data.
*
* @author Holly Vance <[log in to unmask]>
* @author Jeremy McCormick <[log in to unmask]>
@@ -22,19 +21,19 @@
public final class ClusterEnergyCorrection {
// Variables for electron energy corrections.
- static final double par0_em = 0.02203;
- static final double par1_em[] = {35,-0.1795,-0.0005613,16.42,0.3431,-2.021,74.85,-0.3626};
- static final double par2_em[] = {35, 1.011, 0.003234, 18.06, 0.2592, 8.586, 75.08, -0.3771};
+ static final double par0_em = -0.017;
+ static final double par1_em[] = {35,-0.06738,-0.0005613,16.42,0.3431,-2.021,74.85,-0.3626};
+ static final double par2_em[] = {35, 0.933, 0.003234, 18.06, 0.24, 8.586, 75.08, -0.39};
// Variables for positron energy corrections.
- static final double par0_ep = 0.01307;
- static final double par1_ep[] = {35,-0.1415,-0.0008183,17.88,0.2886,-1.192,73.12,-0.3747};
- static final double par2_ep[] = {35, 0.9733, 0.003713, 18.19, 0.2557, 8.342, 72.44, -0.3834};
+ static final double par0_ep = -0.0131;
+ static final double par1_ep[] = {35,-0.076,-0.0008183,17.88,0.2886,-1.192,73.12,-0.3747};
+ static final double par2_ep[] = {35, 0.94, 0.003713, 18.19, 0.24, 8.342, 72.44, -0.39};
// Variables for photon energy corrections.
- static final double par0_p = 0.01604;
- static final double par1_p[] = {35,-0.1268,-0.0008572,16.76,0.2784,-0.07232,72.88,-0.1685};
- static final double par2_p[] = {35, 0.965, 0.004, 18.05, 0.24, 3.027, 74.93, -0.3221};
+ static final double par0_p = -0.0113;
+ static final double par1_p[] = {35,-0.0585,-0.0008572,16.76,0.2784,-0.07232,72.88,-0.1685};
+ static final double par2_p[] = {35, 0.9307, 0.004, 18.05, 0.23, 3.027, 74.93, -0.34};
/**
* Calculate the corrected energy for the cluster.
Modified: java/branches/jeremy-dev/evio/src/main/java/org/hps/evio/RfFitterDriver.java
=============================================================================
--- java/branches/jeremy-dev/evio/src/main/java/org/hps/evio/RfFitterDriver.java (original)
+++ java/branches/jeremy-dev/evio/src/main/java/org/hps/evio/RfFitterDriver.java Tue Jan 26 18:37:32 2016
@@ -43,38 +43,37 @@
IFunction fitFunction=new RfFitFunction();
/*
- * Check the event for an RF pulse, and, if found, fit it to get
- * RF time and then dump it in the lcsim event.
+ * Check the event for an RF pulse, and, if found, fit it to get RF time.
*/
public void process(EventHeader event) {
- if (!event.hasCollection(GenericObject.class,"FADCGenericHits")) return;
+
+ List <RfHit> rfHits=new ArrayList<RfHit>();
boolean foundRf=false;
- double times[]={-9999,-9999};
-
- for (GenericObject gob : event.get(GenericObject.class,"FADCGenericHits")) {
- FADCGenericHit hit=(FADCGenericHit)gob;
-
- // ignore hits not from proper RF signals based on crate/slot/channel:
- if (hit.getCrate()!=CRATE || hit.getSlot()!=SLOT) continue;
- for (int ii=0; ii<CHANNELS.length; ii++) {
- if (hit.getChannel()==CHANNELS[ii]) {
-
- // we found a RF readout, fit it:
- foundRf=true;
- times[ii] = fitPulse(hit);
-
- break;
+ double times[]={-9999,-9999};
+
+ if (event.hasCollection(GenericObject.class,"FADCGenericHits")) {
+
+ for (GenericObject gob : event.get(GenericObject.class,"FADCGenericHits")) {
+ FADCGenericHit hit=(FADCGenericHit)gob;
+
+ // ignore hits not from proper RF signals based on crate/slot/channel:
+ if (hit.getCrate()!=CRATE || hit.getSlot()!=SLOT) continue;
+
+ for (int ii=0; ii<CHANNELS.length; ii++) {
+ if (hit.getChannel()==CHANNELS[ii]) {
+
+ // we found a RF readout, fit it:
+ foundRf=true;
+ times[ii] = fitPulse(hit);
+
+ break;
+ }
}
}
}
-
- // if we found an RF readout, dump the fit result in the event:
- if (foundRf) {
- List <RfHit> rfHits=new ArrayList<RfHit>();
- rfHits.add(new RfHit(times));
- event.put("RFHits", rfHits, RfHit.class, 1);
- }
+ if (foundRf) rfHits.add(new RfHit(times));
+ event.put("RFHits", rfHits, RfHit.class, 1);
}
/*
Modified: java/branches/jeremy-dev/monitoring-app/src/main/java/org/hps/monitoring/application/EventProcessing.java
=============================================================================
--- java/branches/jeremy-dev/monitoring-app/src/main/java/org/hps/monitoring/application/EventProcessing.java (original)
+++ java/branches/jeremy-dev/monitoring-app/src/main/java/org/hps/monitoring/application/EventProcessing.java Tue Jan 26 18:37:32 2016
@@ -464,9 +464,12 @@
steering = configurationModel.getSteeringFile();
} else {
steering = configurationModel.getSteeringResource();
- }
-
- this.logger.config("set steering " + steering + " with type "
+ if (!steering.startsWith("/")) {
+ steering = "/" + steering;
+ }
+ }
+
+ this.logger.config("set steering " + steering + " with type "
+ (steeringType == SteeringType.RESOURCE ? "RESOURCE" : "FILE"));
try {
Modified: java/branches/jeremy-dev/monitoring-drivers/src/main/java/org/hps/monitoring/drivers/trackrecon/TrackingReconPlots.java
=============================================================================
--- java/branches/jeremy-dev/monitoring-drivers/src/main/java/org/hps/monitoring/drivers/trackrecon/TrackingReconPlots.java (original)
+++ java/branches/jeremy-dev/monitoring-drivers/src/main/java/org/hps/monitoring/drivers/trackrecon/TrackingReconPlots.java Tue Jan 26 18:37:32 2016
@@ -12,14 +12,17 @@
import java.util.List;
import java.util.logging.Level;
import java.util.logging.Logger;
+
import static org.hps.monitoring.drivers.trackrecon.PlotAndFitUtilities.plot;
-import org.hps.recon.tracking.HPSTrack;
+
+import org.hps.recon.tracking.HpsHelicalTrackFit;
+import org.hps.recon.tracking.TrackUtils;
import org.lcsim.event.Cluster;
-
import org.lcsim.event.EventHeader;
import org.lcsim.event.LCIOParameters.ParameterName;
import org.lcsim.event.RawTrackerHit;
import org.lcsim.event.Track;
+import org.lcsim.event.TrackState;
import org.lcsim.event.TrackerHit;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.geometry.Detector;
@@ -267,10 +270,12 @@
SeedTrack stEle = (SeedTrack) trk;
SeedCandidate seedEle = stEle.getSeedCandidate();
HelicalTrackFit ht = seedEle.getHelix();
- HPSTrack hpstrk = new HPSTrack(ht);
+ HpsHelicalTrackFit hpstrk = new HpsHelicalTrackFit(ht);
double svt_l12 = 900.00;//mm ~approximately...this doesn't matter much
double ecal_face = 1393.00;//mm ~approximately ... this matters! Should use typical shower depth...or, once have cluster match, use that value of Z
- Hep3Vector posAtEcal = hpstrk.getPositionAtZMap(svt_l12, ecal_face, 5.0, event.getDetector().getFieldMap())[0];
+ TrackState stateAtEcal = TrackUtils.extrapolateTrackUsingFieldMap(trk, svt_l12, ecal_face, 5.0, event.getDetector().getFieldMap());
+ Hep3Vector posAtEcal = new BasicHep3Vector(stateAtEcal.getReferencePoint());
+ //Hep3Vector posAtEcal = hpstrk.getPositionAtZMap(svt_l12, ecal_face, 5.0, event.getDetector().getFieldMap())[0];
List<Cluster> clusters = event.get(Cluster.class, ecalCollectionName);
if (clusters != null) {
if (debug)
@@ -282,7 +287,8 @@
Hep3Vector clusterPos = new BasicHep3Vector(clust.getPosition());
double zCluster = clusterPos.z();
//improve the extrapolation...use the reconstructed cluster z-position
- posAtEcal = hpstrk.getPositionAtZMap(svt_l12, zCluster, 5.0, event.getDetector().getFieldMap())[0];
+ stateAtEcal = TrackUtils.extrapolateTrackUsingFieldMap(trk, svt_l12, zCluster, 5.0, event.getDetector().getFieldMap());
+ posAtEcal = new BasicHep3Vector(stateAtEcal.getReferencePoint());
double eOverP = clust.getEnergy() / pmag;
double dx = posAtEcal.x() - clusterPos.x();
double dy = posAtEcal.y() - clusterPos.y();
Modified: java/branches/jeremy-dev/recon/src/main/java/org/hps/recon/vertexing/TwoTrackFringeVertexer.java
=============================================================================
--- java/branches/jeremy-dev/recon/src/main/java/org/hps/recon/vertexing/TwoTrackFringeVertexer.java (original)
+++ java/branches/jeremy-dev/recon/src/main/java/org/hps/recon/vertexing/TwoTrackFringeVertexer.java Tue Jan 26 18:37:32 2016
@@ -4,14 +4,14 @@
import hep.physics.vec.Hep3Vector;
import org.hps.recon.tracking.BeamlineConstants;
-import org.hps.recon.tracking.HPSTrack;
+import org.hps.recon.tracking.HpsHelicalTrackFit;
import org.hps.recon.tracking.HelixConverter;
-import org.hps.recon.tracking.StraightLineTrack;
+import org.hps.recon.tracking.TrackUtils;
import org.lcsim.event.Track;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.fit.helicaltrack.HelixUtils;
+import org.lcsim.geometry.FieldMap;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
-import org.lcsim.util.swim.Helix;
/**
*
@@ -20,38 +20,42 @@
* @author phansson
*
*/
+
public class TwoTrackFringeVertexer extends TwoTrackVertexer {
protected HelixConverter converter = new HelixConverter(0.);
- public void setTracks(Track track1, Track track2) {
+ public void setTracks(Track track1, Track track2, FieldMap fieldMap) {
SeedTrack s1 = (SeedTrack) track1;
HelicalTrackFit htf1 = s1.getSeedCandidate().getHelix();
- HPSTrack hpstrk1 = new HPSTrack(htf1);
+ HpsHelicalTrackFit hpstrk1 = new HpsHelicalTrackFit(htf1);
SeedTrack s2 = (SeedTrack) track2;
HelicalTrackFit htf2 = s2.getSeedCandidate().getHelix();
- HPSTrack hpstrk2 = new HPSTrack(htf2);
+ HpsHelicalTrackFit hpstrk2 = new HpsHelicalTrackFit(htf2);
boolean debug = false;
+
+ Hep3Vector posAtConv1 = new BasicHep3Vector( TrackUtils.extrapolateTrackUsingFieldMap(track1, 100.0, BeamlineConstants.HARP_POSITION_TESTRUN, 5.0, fieldMap).getReferencePoint());
+ Hep3Vector posAtConv2 = new BasicHep3Vector( TrackUtils.extrapolateTrackUsingFieldMap(track2, 100.0, BeamlineConstants.HARP_POSITION_TESTRUN, 5.0, fieldMap).getReferencePoint());
+
+ //FIXME the straight line objects are not working
- Hep3Vector posAtConv1 = hpstrk1.getPositionAtZMap(100.0, BeamlineConstants.HARP_POSITION_TESTRUN, 5.0)[0];
- Hep3Vector posAtConv2 = hpstrk2.getPositionAtZMap(100.0, BeamlineConstants.HARP_POSITION_TESTRUN, 5.0)[0];
+ //StraightLineTrack slt1_conv = converter.Convert((Helix)hpstrk1.getTrajectory());
+ //StraightLineTrack slt2_conv = converter.Convert((Helix)hpstrk2.getTrajectory());
+
+ //A1 = new BasicHep3Vector(slt1_conv.x0(),slt1_conv.y0(),slt1_conv.z0());
+ //B1 = new BasicHep3Vector(slt2_conv.x0(),slt2_conv.y0(),slt2_conv.z0());
- StraightLineTrack slt1_conv = converter.Convert((Helix)hpstrk1.getTrajectory());
- StraightLineTrack slt2_conv = converter.Convert((Helix)hpstrk2.getTrajectory());
+ //double YZAtConv1[] = slt1_conv.getYZAtX(BeamlineConstants.HARP_POSITION_TESTRUN);
+ //double YZAtConv2[] = slt2_conv.getYZAtX(BeamlineConstants.HARP_POSITION_TESTRUN);
- A1 = new BasicHep3Vector(slt1_conv.x0(),slt1_conv.y0(),slt1_conv.z0());
- B1 = new BasicHep3Vector(slt2_conv.x0(),slt2_conv.y0(),slt2_conv.z0());
-
- double YZAtConv1[] = slt1_conv.getYZAtX(BeamlineConstants.HARP_POSITION_TESTRUN);
- double YZAtConv2[] = slt2_conv.getYZAtX(BeamlineConstants.HARP_POSITION_TESTRUN);
+ //A2 = new BasicHep3Vector(BeamlineConstants.HARP_POSITION_TESTRUN,YZAtConv1[0],YZAtConv1[1]);
+ //B2 = new BasicHep3Vector(BeamlineConstants.HARP_POSITION_TESTRUN,YZAtConv2[0],YZAtConv2[1]);
- A2 = new BasicHep3Vector(BeamlineConstants.HARP_POSITION_TESTRUN,YZAtConv1[0],YZAtConv1[1]);
- B2 = new BasicHep3Vector(BeamlineConstants.HARP_POSITION_TESTRUN,YZAtConv2[0],YZAtConv2[1]);
if(debug) {
System.out.printf("%s: original track1 direction at x=0 %s \n",this.getClass().getSimpleName(),HelixUtils.Direction(hpstrk1,0.).toString());
System.out.printf("%s: original track2 direction at x=0 %s \n",this.getClass().getSimpleName(),HelixUtils.Direction(hpstrk2,0.).toString());
- System.out.printf("%s: track1 direction at conv %s \n",this.getClass().getSimpleName(),hpstrk1.getTrajectory().getUnitTangentAtLength(0.).toString());
- System.out.printf("%s: track2 direction at conv %s \n",this.getClass().getSimpleName(),hpstrk2.getTrajectory().getUnitTangentAtLength(0.).toString());
+ //System.out.printf("%s: track1 direction at conv %s \n",this.getClass().getSimpleName(),hpstrk1.getTrajectory().getUnitTangentAtLength(0.).toString());
+ //System.out.printf("%s: track2 direction at conv %s \n",this.getClass().getSimpleName(),hpstrk2.getTrajectory().getUnitTangentAtLength(0.).toString());
System.out.printf("%s: pos at converter track1 %s \n",this.getClass().getSimpleName(),posAtConv1.toString());
Modified: java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/production/FilterHitsPerSensor.lcsim
=============================================================================
--- java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/production/FilterHitsPerSensor.lcsim (original)
+++ java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/production/FilterHitsPerSensor.lcsim Tue Jan 26 18:37:32 2016
@@ -7,14 +7,14 @@
<execute>
<driver name="CleanupDriver"/>
<driver name="EventMarkerDriver"/>
+ <driver name="RawHitsPerSensorFilter"/>
<!-- SVT reconstruction drivers -->
+ <!--
<driver name="RawTrackerHitSensorSetup"/>
<driver name="RawTrackerHitFitterDriver" />
<driver name="TrackerHitDriver"/>
-
- <!-- filter driver -->
- <driver name="HitsPerSensorFilter"/>
-
+ <driver name="ClustersPerSensorFilter"/>
+ -->
<!-- write and cleanup -->
<driver name="LCIOWriter"/>
</execute>
@@ -39,7 +39,11 @@
<driver name="TrackerHitDriver" type="org.hps.recon.tracking.DataTrackerHitDriver">
<neighborDeltaT>8.0</neighborDeltaT>
</driver>
- <driver name="HitsPerSensorFilter" type="org.hps.recon.filtering.SvtHitMultiplicityFilter"/>
+ <driver name="RawHitsPerSensorFilter" type="org.hps.recon.filtering.SvtRawHitMultiplicityFilter">
+ <minHitsPerSensor>1</minHitsPerSensor>
+ <maxHitsPerSensor>2</maxHitsPerSensor>
+ </driver>
+ <driver name="ClustersPerSensorFilter" type="org.hps.recon.filtering.SvtHitMultiplicityFilter"/>
<driver name="LCIOWriter" type="org.lcsim.util.loop.LCIODriver">
<outputFilePath>${outputFile}.slcio</outputFilePath>
</driver>
Modified: java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/recon/EngineeringRun2014EcalReconMC.lcsim
=============================================================================
--- java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/recon/EngineeringRun2014EcalReconMC.lcsim (original)
+++ java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/recon/EngineeringRun2014EcalReconMC.lcsim Tue Jan 26 18:37:32 2016
@@ -18,6 +18,7 @@
<execute>
<!--<driver name="EventMarkerDriver" />-->
<driver name="EcalRunningPedestal"/>
+ <driver name="RfFitter"/>
<driver name="EcalRawConverter" />
<driver name="ReconClusterer" />
<driver name="GTPOnlineClusterer" />
@@ -28,6 +29,8 @@
<driver name="EventMarkerDriver" type="org.lcsim.job.EventMarkerDriver">
<eventInterval>100</eventInterval>
</driver>
+ <driver name="RfFitter" type="org.hps.evio.RfFitterDriver"/>
+
<driver name="ConditionsDriver" type="org.hps.conditions.ConditionsDriver">
<detectorName>HPS-ECalCommissioning-v3-fieldmap</detectorName>
<runNumber>3422</runNumber>
Modified: java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/recon/EngineeringRun2015FullReconMC.lcsim
=============================================================================
--- java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/recon/EngineeringRun2015FullReconMC.lcsim (original)
+++ java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/recon/EngineeringRun2015FullReconMC.lcsim Tue Jan 26 18:37:32 2016
@@ -11,6 +11,7 @@
<driver name="ConditionsDriver"/>
<driver name="EventMarkerDriver"/>
+ <driver name="RfFitter"/>
<!-- Ecal reconstruction drivers -->
<driver name="EcalRawConverter" />
<driver name="ReconClusterer" />
@@ -68,7 +69,9 @@
<driver name="EventMarkerDriver" type="org.lcsim.job.EventMarkerDriver">
<eventInterval>1000</eventInterval>
- </driver>
+ </driver>
+ <driver name="RfFitter" type="org.hps.evio.RfFitterDriver"/>
+
<driver name="RawTrackerHitSensorSetup" type="org.lcsim.recon.tracking.digitization.sisim.config.RawTrackerHitSensorSetup">
<readoutCollections>SVTRawTrackerHits</readoutCollections>
</driver>
Modified: java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/users/phansson/EngineeringRun2015FullRecon_Pass2_Gbl.lcsim
=============================================================================
--- java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/users/phansson/EngineeringRun2015FullRecon_Pass2_Gbl.lcsim (original)
+++ java/branches/jeremy-dev/steering-files/src/main/resources/org/hps/steering/users/phansson/EngineeringRun2015FullRecon_Pass2_Gbl.lcsim Tue Jan 26 18:37:32 2016
@@ -56,7 +56,7 @@
<driver name="GBLOutputDriver" />
<driver name="GBLRefitterDriver" />
<driver name="TrackDataDriver" />
- <!--<driver name="GblResidualEcalDriver" />-->
+ <driver name="GblResidualEcalDriver" />
<!--<driver name="TrackExtrapolationTestDriver"/>-->
<driver name="ReconParticleDriver" />
<!--<driver name="TrackingReconstructionPlots" />-->
@@ -162,13 +162,13 @@
<outputFileName>${outputFile}.root</outputFileName>
</driver>
<driver name="AlignmentFilterDriver" type="org.hps.recon.filtering.SvtAlignmentFilter"/>
- <driver name="GblResidualEcalDriver" type="org.hps.users.phansson.ECalExtrapolationDriver"/>
+ <driver name="GblResidualEcalDriver" type="org.hps.users.phansson.GblResidualDriver"/>
<driver name="TrackExtrapolationTestDriver" type="org.hps.users.phansson.TrackExtrapolationTestDriver"/>
<driver name="TrackingReconstructionPlots" type="org.hps.users.phansson.TrackingReconstructionPlots">
<showPlots>False</showPlots>
</driver>
<driver name="TimerDriver1" type="org.hps.util.TimerDriver"/>
- <driver name="GeomChecker" type="org.hps.users.phansson.TrackingGeometryChecker"/>
+ <driver name="GeomChecker" type="org.hps.users.phansson.tools.TrackingGeometryChecker"/>
</drivers>
</lcsim>
Modified: java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
=============================================================================
--- java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java (original)
+++ java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java Tue Jan 26 18:37:32 2016
@@ -6,6 +6,7 @@
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.SpacePoint;
import hep.physics.vec.VecOp;
+
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
@@ -15,10 +16,13 @@
import java.util.List;
import java.util.Map;
import java.util.Set;
+
import org.apache.commons.math3.util.Pair;
import org.hps.recon.tracking.EventQuality.Quality;
import org.hps.recon.tracking.gbl.HelicalTrackStripGbl;
+
import static org.lcsim.constants.Constants.fieldConversion;
+
import org.lcsim.detector.ITransform3D;
import org.lcsim.detector.solids.Box;
import org.lcsim.detector.solids.Point3D;
@@ -86,6 +90,67 @@
return extrapolateHelixToXPlane(getHTF(track), x);
}
+ /**
+ * Change reference point of helix (following L3 Internal Note 1666.)
+ * @param newRefPoint - The new reference point in XY
+ */
+ public static double[] getParametersAtNewRefPoint(double[] newRefPoint, HpsHelicalTrackFit helicalTrackFit) {
+ return getParametersAtNewRefPoint(newRefPoint, helicalTrackFit.getRefPoint(),helicalTrackFit.parameters());
+ }
+
+ /**
+ * Change reference point of helix (following L3 Internal Note 1666.)
+ * @param newRefPoint - The new reference point in XY
+ */
+ public static double[] getParametersAtNewRefPoint(double[] newRefPoint, double[] __refPoint,
+ double[] parameters) {
+
+ double phi0 = parameters[HelicalTrackFit.phi0Index];
+ double curvature =parameters[HelicalTrackFit.curvatureIndex];
+ double dca = parameters[HelicalTrackFit.dcaIndex];
+ double slope = parameters[HelicalTrackFit.slopeIndex];
+ double z0 = parameters[HelicalTrackFit.z0Index];
+
+ //take care of phi0 range if needed (this matters for dphi below I think)
+ // L3 defines it in the range [-pi,pi]
+ if(phi0 > Math.PI)
+ phi0 -= Math.PI*2;
+
+ double dx = newRefPoint[0] - __refPoint[0];
+ double dy = newRefPoint[1] - __refPoint[1];
+ double sinphi = Math.sin(phi0);
+ double cosphi = Math.cos(phi0);
+ double R = 1.0/curvature;
+
+ // calculate new phi
+ double phinew = Math.atan2( sinphi - dx/(R-dca) , cosphi + dy/(R-dca) );
+
+ // difference in phi
+ // watch out for ambiguity
+ double dphi = phinew - phi0;
+ if (Math.abs( dphi ) > Math.PI)
+ throw new RuntimeException("dphi is large " + dphi + " from phi0 " + phi0
+ + " and phinew " + phinew + " take care of the ambiguity!!??");
+
+ // calculate new dca
+ double dcanew = dca + dx*sinphi - dy*cosphi + (dx*cosphi + dy*sinphi)*Math.tan( dphi/2. );
+
+ // path length from old to new point
+ double s = -1.0*dphi/curvature;
+
+ // new z0
+ double z0new = z0 + s*slope;
+
+ // new array
+ double[] params = new double[5];
+ params[HelicalTrackFit.phi0Index] = phinew;
+ params[HelicalTrackFit.curvatureIndex] = curvature;
+ params[HelicalTrackFit.dcaIndex] = dcanew;
+ params[HelicalTrackFit.slopeIndex] = slope;
+ params[HelicalTrackFit.z0Index] = z0new;
+ return params;
+ }
+
/**
* Extrapolate helix to a position along the x-axis. Re-use HelixUtils.
@@ -801,11 +866,11 @@
/**
- * 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
- * @param org origin to be used for the track
+ * Transform MCParticle into a {@link HelicalTrackFit} object. Note that it produces the {@link HelicalTrackFit}
+ * parameters at nominal reference point at origin and assumes that there is no field at x<0
+ *
+ * @param mcp - MC particle to be transformed
+ * @param origin - origin to be used for the track
* @return {@link HelicalTrackFit} object based on the MC particle
*/
public static HelicalTrackFit getHTF(MCParticle mcp, Hep3Vector origin, double Bz) {
@@ -871,27 +936,7 @@
return htf;
}
- public static StraightLineTrack findSLTAtZ(Track trk1, double zVal, boolean useFringe) {
- SeedTrack s1 = (SeedTrack) trk1;
- HelicalTrackFit htf1 = s1.getSeedCandidate().getHelix();
- HPSTrack hpstrk1 = new HPSTrack(htf1);
- Hep3Vector pos1;
- if (useFringe) {
- // broken because you need ot provide the Field Map to get this...
-// pos1 = hpstrk1.getPositionAtZMap(100.0, zVal, 5.0)[0];
- } else
- pos1 = TrackUtils.extrapolateTrack(trk1, zVal);
- // System.out.printf("%s: Position1 at edge of fringe %s\n",this.getClass().getSimpleName(),pos1.toString());
- Helix traj = (Helix) hpstrk1.getTrajectory();
- if (traj == null) {
- SpacePoint r0 = new SpacePoint(HelixUtils.PointOnHelix(htf1, 0));
- traj = new Helix(r0, htf1.R(), htf1.phi0(), Math.atan(htf1.slope()));
- }
- HelixConverter converter = new HelixConverter(0.);
- StraightLineTrack slt1 = converter.Convert(traj);
- // System.out.printf("%s: straight line track: x0=%f,y0=%f,z0=%f dz/dx=%f dydx=%f targetY=%f targetZ=%f \n",this.getClass().getSimpleName(),slt1.x0(),slt1.y0(),slt1.z0(),slt1.dzdx(),slt1.dydx(),slt1.TargetYZ()[0],slt1.TargetYZ()[1]);
- return slt1;
- }
+
public static MCParticle getMatchedTruthParticle(Track track) {
boolean debug = false;
Modified: java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/WTrack.java
=============================================================================
--- java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/WTrack.java (original)
+++ java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/WTrack.java Tue Jan 26 18:37:32 2016
@@ -12,9 +12,13 @@
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
/**
- * Track parameterization representation.
- *
- * @author phansson
+ * Track representation based on paper
+ * Paul Avery, CBX 98-39, June 9, 1998
+ *
+ * Used primarily for the algorithm to intersect a helix with a generic plane in space.
+ *
+ * @author phansson <[log in to unmask]>
+
*/
public class WTrack {
Modified: java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java
=============================================================================
--- java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java (original)
+++ java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java Tue Jan 26 18:37:32 2016
@@ -1,6 +1,20 @@
package org.hps.recon.tracking.gbl;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.Map;
+import java.util.logging.Logger;
+
+import org.apache.commons.math3.util.Pair;
+import org.hps.recon.tracking.HpsHelicalTrackFit;
+import org.hps.recon.tracking.TrackUtils;
+import org.hps.recon.tracking.gbl.matrix.Matrix;
+import org.hps.recon.tracking.gbl.matrix.SymMatrix;
+import org.hps.recon.tracking.gbl.matrix.Vector;
import org.lcsim.event.Track;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelixUtils;
/**
* A class that collects information about a fitted GBL trajectory.
@@ -9,6 +23,9 @@
*
*/
public class FittedGblTrajectory {
+
+ public static Logger LOGGER = Logger.getLogger(FittedGblTrajectory.class.getName());
+
public enum GBLPOINT {
IP(0), LAST(1), VERTEX(2);
private int numVal;
@@ -34,6 +51,7 @@
}
}
+
public static enum GBLPARIDX {
QOVERP(0),YTPRIME(1),XTPRIME(2),XT(3),YT(4);
private int _value;
@@ -44,18 +62,101 @@
return _value;
}
};
+
private GblTrajectory _traj;
private double _chi2;
private double _lost;
private int _ndf;
private Track _seed = null;
- private GBLTrackData _t = null;
+ private Map<Integer, Double> pathLengthMap = null;
+
+ /**
+ * Default constructor.
+ *
+ * @param traj
+ * @param chi2
+ * @param ndf
+ * @param lost
+ */
public FittedGblTrajectory(GblTrajectory traj, double chi2, int ndf, double lost) {
_traj = traj;
_chi2 = chi2;
_ndf = ndf;
_lost = lost;
}
+
+ /**
+ * Find the index (or label) of the GBL point on the trajectory from the {@link GBLPOINT}.
+ * @param point
+ * @return
+ */
+ public int getPointIndex(GBLPOINT point) {
+ int gblPointIndex;
+ if (point.compareTo(GBLPOINT.IP) == 0)
+ gblPointIndex = 1;
+ else if (point.compareTo(GBLPOINT.LAST) == 0)
+ gblPointIndex = _traj.getNumPoints();
+ else
+ throw new RuntimeException("This GBL point " + point.toString() + "( " + point.name() + ") is not valid");
+ return gblPointIndex;
+ }
+
+
+ /**
+ * Find the corrections and covariance matrix for a particular {@link GBLPOINT}
+ * @param point
+ * @param locPar
+ * @param locCov
+ */
+ public void getResults(GBLPOINT point, Vector locPar, SymMatrix locCov) {
+
+ // find the GBL point index
+ int gblPointIndex = getPointIndex(point);
+
+ // get the results
+ getResults(gblPointIndex, locPar, locCov);
+
+ }
+
+
+ /**
+ * Find the corrections and covariance matrix for a particular point on the GBL trajectory
+ * @param iLabel
+ * @param locPar
+ * @param locCov
+ */
+ public void getResults(int iLabel, Vector locPar, SymMatrix locCov) {
+
+ // Get the result from the trajectory
+ int ok = _traj.getResults(iLabel, locPar, locCov);
+
+ // check that the fit was ok
+ if( ok != 0)
+ throw new RuntimeException("Trying to extract GBL corrections for fit that failed!?");
+ }
+
+
+ /**
+ * Find the path length to this point.
+ * @param point - {@link GBLPOINT} point
+ * @return path length
+ */
+ public double getPathLength(GBLPOINT point) {
+ int gblPointIndex = getPointIndex(point);
+ return getPathLength(gblPointIndex);
+ }
+
+ /**
+ * Find the path length to this point.
+ * @param iLabel - GBL point index
+ * @return path length
+ */
+ public double getPathLength(int iLabel) {
+ if( !this.pathLengthMap.containsKey(iLabel) )
+ throw new RuntimeException("This iLabel " + iLabel + " doesn't exists in the path length map.");
+ return this.pathLengthMap.get(iLabel);
+ }
+
public void set_seed(Track seed) {
_seed = seed;
}
@@ -74,5 +175,149 @@
public int get_ndf() {
return _ndf;
}
+
+ public void setPathLengthMap(Map<Integer, Double> pathLengthMap) {
+ this.pathLengthMap = pathLengthMap;
+ }
+
+ public Map<Integer, Double> getPathLengthMap() {
+ if (this.pathLengthMap == null)
+ throw new RuntimeException("No path length map has been set on this trajectory!");
+ return this.pathLengthMap;
+ }
+
+
+
+ /**
+ * Get the corrected perigee parameters and covariance matrix for a point on the {@link GblTrajectory}.
+ *
+ * FIXME the covariance matrix is not properly propagated along the trajectory right now!
+ *
+ * @param htf - helix to be corrected
+ * @param point - {@link GBLPOINT} on the trajectory
+ * @param bfield - magnitude of B-field.
+ * @return
+ */
+ public Pair<double[], SymmetricMatrix> getCorrectedPerigeeParameters(HelicalTrackFit htf, GBLPOINT point, double bfield) {
+
+ // find the point on the trajectory from the GBLPOINT
+ int iLabel = getPointIndex(point);
+
+ return getCorrectedPerigeeParameters(htf, iLabel, bfield);
+
+ }
+
+
+ /**
+ * Get the corrected perigee parameters and covariance matrix for a point on the {@link GblTrajectory}.
+ *
+ * FIXME the covariance matrix is not properly propagated along the trajectory right now!
+ *
+ * @param htf - helix to be corrected
+ * @param iLabel - label of the point on the {@link GblTrajectory}
+ * @param bfield - magnitude of B-field.
+ * @return the corrected perigee parameters
+ */
+ public Pair<double[], SymmetricMatrix> getCorrectedPerigeeParameters(HelicalTrackFit htf, int iLabel, double bfield) {
+
+ // Get corrections from GBL fit
+ Vector locPar = new Vector(5);
+ SymMatrix locCov = new SymMatrix(5);
+
+ // Extract the corrections to the track parameters and the covariance matrix from the GBL trajectory
+ getResults(iLabel, locPar, locCov);
+
+ // Use the super class to keep track of reference point of the helix
+ HpsHelicalTrackFit helicalTrackFit = new HpsHelicalTrackFit(htf);
+ double[] refIP = helicalTrackFit.getRefPoint();
+
+ // Calculate new reference point for this point
+ // This is the intersection of the helix with the plane
+ // The trajectory has this information already in the form of a map between GBL point and path length
+ double pathLength = getPathLength(iLabel);
+ Hep3Vector refPointVec = HelixUtils.PointOnHelix(helicalTrackFit, pathLength);
+ double[] refPoint = new double[]{refPointVec.x(), refPointVec.y()};
+
+ LOGGER.finest("pathLength " + pathLength + " -> refPointVec " + refPointVec.toString());
+
+ // Propagate the helix to new reference point
+ double[] helixParametersAtPoint = TrackUtils.getParametersAtNewRefPoint(refPoint, helicalTrackFit);
+
+ // Create a new helix with the new parameters and the new reference point
+ HpsHelicalTrackFit helicalTrackFitAtPoint = new HpsHelicalTrackFit(helixParametersAtPoint, helicalTrackFit.covariance(),
+ helicalTrackFit.chisq(), helicalTrackFit.ndf(), helicalTrackFit.PathMap(),
+ helicalTrackFit.ScatterMap(), refPoint);
+
+ // find the corrected perigee track parameters at this point
+ double[] helixParametersAtPointCorrected = GblUtils.getCorrectedPerigeeParameters(locPar, helicalTrackFitAtPoint, bfield);
+
+ // create a new helix
+ HpsHelicalTrackFit helicalTrackFitAtPointCorrected = new HpsHelicalTrackFit(helixParametersAtPointCorrected, helicalTrackFit.covariance(),
+ helicalTrackFit.chisq(), helicalTrackFit.ndf(), helicalTrackFit.PathMap(),
+ helicalTrackFit.ScatterMap(), refPoint);
+
+ // change reference point back to the original one
+ double[] helixParametersAtIPCorrected = TrackUtils.getParametersAtNewRefPoint(refIP, helicalTrackFitAtPointCorrected);
+
+ // create a new helix for the new parameters at the IP reference point
+ HpsHelicalTrackFit helicalTrackFitAtIPCorrected = new HpsHelicalTrackFit(helixParametersAtIPCorrected, helicalTrackFit.covariance(),
+ helicalTrackFit.chisq(), helicalTrackFit.ndf(), helicalTrackFit.PathMap(),
+ helicalTrackFit.ScatterMap(), refIP);
+
+
+ // Calculate the updated covariance
+ Matrix jacobian = GblUtils.getCLToPerigeeJacobian(helicalTrackFit, helicalTrackFitAtIPCorrected, bfield);
+ Matrix helixCovariance = jacobian.times(locCov.times(jacobian.transpose()));
+ SymmetricMatrix cov = new SymmetricMatrix(5);
+ for (int i = 0; i < 5; i++) {
+ for (int j = 0; j < 5; j++) {
+ if (i >= j) {
+ cov.setElement(i, j, helixCovariance.get(i, j));
+ }
+ }
+ }
+ LOGGER.finest("corrected helix covariance:\n" + cov);
+
+ double parameters_gbl[] = helicalTrackFitAtIPCorrected.parameters();
+
+ return new Pair<double[], SymmetricMatrix>(parameters_gbl,cov);
+ }
+
+
+
+ /**
+ * Extract kinks across the trajectory.
+ * @return kinks in a {@link GBLKinkData} object.
+ */
+ public GBLKinkData getKinks() {
+ GblTrajectory traj = this._traj;
+ // get corrections from GBL fit
+ Vector locPar = new Vector(5);
+ SymMatrix locCov = new SymMatrix(5);
+ float[] lambdaKinks = new float[traj.getNumPoints() - 1];
+ double[] phiKinks = new double[traj.getNumPoints() - 1];
+
+ double oldPhi = 0, oldLambda = 0;
+ for (int i = 0; i < traj.getNumPoints(); i++) {
+ traj.getResults(i + 1, locPar, locCov); // vertex point
+ double newPhi = locPar.get(GBLPARIDX.XTPRIME.getValue());
+ double newLambda = locPar.get(GBLPARIDX.YTPRIME.getValue());
+ if (i > 0) {
+ lambdaKinks[i - 1] = (float) (newLambda - oldLambda);
+ phiKinks[i - 1] = newPhi - oldPhi;
+ // System.out.println("phikink: " + (newPhi - oldPhi));
+ }
+ oldPhi = newPhi;
+ oldLambda = newLambda;
+ }
+
+ return new GBLKinkData(lambdaKinks, phiKinks);
+ }
+
+
+
+
+
+
}
Modified: java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLFileIO.java
=============================================================================
--- java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLFileIO.java (original)
+++ java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLFileIO.java Tue Jan 26 18:37:32 2016
@@ -11,8 +11,8 @@
import java.util.logging.Level;
import java.util.logging.Logger;
-import org.hps.recon.tracking.gbl.GBLOutput.ClParams;
-import org.hps.recon.tracking.gbl.GBLOutput.PerigeeParams;
+import org.hps.recon.tracking.gbl.GblUtils.ClParams;
+import org.hps.recon.tracking.gbl.GblUtils.PerigeeParams;
import org.hps.svt.alignment.RunAlignment;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
@@ -67,34 +67,34 @@
addLine(String.format("Track perPar (R phi0 slope d0 z0) %.12f %.12f %.12f %.12f %.12f",htf.R(),htf.phi0(),htf.slope(),htf.dca(),htf.z0()));
}
- String getPerTrackParamStr(PerigeeParams perPar) {
+ String getPerTrackParamStr(GblUtils.PerigeeParams perPar) {
return String.format("Track perPar (R theta phi d0 z0) %.12f %.12f %.12f %.12f %.12f",1.0/perPar.getKappa(),perPar.getTheta(),perPar.getPhi(),perPar.getD0(),perPar.getZ0());
}
- void printPerTrackParam(PerigeeParams perPar) {
+ void printPerTrackParam(GblUtils.PerigeeParams perPar) {
addLine(this.getPerTrackParamStr(perPar));
}
- String getPerTrackParamTruthStr(PerigeeParams perPar) {
+ String getPerTrackParamTruthStr(GblUtils.PerigeeParams perPar) {
return String.format("Truth perPar (kappa theta phi d0 z0) %.12f %.12f %.12f %.12f %.12f",perPar.getKappa(),perPar.getTheta(),perPar.getPhi(),perPar.getD0(),perPar.getZ0());
}
- void printPerTrackParamTruth(PerigeeParams perPar) {
+ void printPerTrackParamTruth(GblUtils.PerigeeParams perPar) {
addLine(this.getPerTrackParamTruthStr(perPar));
}
- String getClTrackParamTruthStr(ClParams perPar) {
+ String getClTrackParamTruthStr(GblUtils.ClParams perPar) {
return String.format("Truth clPar (q/p lambda phi xT yT) %.12f %.12f %.12f %.12f %.12f",perPar.getQoverP(),perPar.getLambda(),perPar.getPhi(),perPar.getXt(),perPar.getYt());
}
- void printClTrackParamTruth(ClParams perPar) {
+ void printClTrackParamTruth(GblUtils.ClParams perPar) {
addLine(this.getClTrackParamTruthStr(perPar));
}
- String getClTrackParamStr(ClParams perPar) {
+ String getClTrackParamStr(GblUtils.ClParams perPar) {
return String.format("Track clPar (q/p lambda phi xT yT) %.12f %.12f %.12f %.12f %.12f",perPar.getQoverP(),perPar.getLambda(),perPar.getPhi(),perPar.getXt(),perPar.getYt());
}
- void printClTrackParam(ClParams perPar) {
+ void printClTrackParam(GblUtils.ClParams perPar) {
addLine(String.format("%s",this.getClTrackParamStr(perPar)));
}
Modified: java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java
=============================================================================
--- java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java (original)
+++ java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java Tue Jan 26 18:37:32 2016
@@ -169,7 +169,7 @@
if (isMC) {
// find the truth particle for this track
- mcp = getMatchedTruthParticle(trk);
+ mcp = TrackUtils.getMatchedTruthParticle(trk);
// check if this is an A' event
for(MCParticle part : mcParticles) {
@@ -222,8 +222,8 @@
// Use the truth helix as the initial track for GBL?
//htf = htfTruth;
// Get perigee parameters to curvilinear frame
- PerigeeParams perPar = new PerigeeParams(htf, bFieldVector.z());
- PerigeeParams perParTruth = new PerigeeParams(htfTruth, bFieldVector.z());
+ GblUtils.PerigeeParams perPar = new GblUtils.PerigeeParams(htf, bFieldVector.z());
+ GblUtils.PerigeeParams perParTruth = new GblUtils.PerigeeParams(htfTruth, bFieldVector.z());
//GBLDATA
gtd.setPerigeeTrackParameters(perPar);
@@ -234,8 +234,8 @@
// Get curvilinear parameters
if (textFile != null) {
- ClParams clPar = new ClParams(htf, bFieldVector.z());
- ClParams clParTruth = new ClParams(htfTruth, bFieldVector.z());
+ GblUtils.ClParams clPar = new GblUtils.ClParams(htf, bFieldVector.z());
+ GblUtils.ClParams clParTruth = new GblUtils.ClParams(htfTruth, bFieldVector.z());
textFile.printClTrackParam(clPar);
textFile.printClTrackParamTruth(clParTruth);
@@ -246,7 +246,7 @@
}
// find the projection from the I,J,K to U,V,T curvilinear coordinates
- Hep3Matrix perToClPrj = getPerToClPrj(htf);
+ Hep3Matrix perToClPrj = GblUtils.getPerToClPrj(htf);
//GBLDATA
for (int row = 0; row < perToClPrj.getNRows(); ++row) {
@@ -780,171 +780,10 @@
}
}
- 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)) {
- 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; //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");
- }
- }
- return maxEntry == null ? null : maxEntry.getKey();
- }
-
-// 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;
-//
-// }
- /**
- * 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
- */
-// private HelicalTrackFit getHTF(MCParticle mcp) {
-// Hep3Vector org = this._hpstrans.transformVectorToTracking(mcp.getOrigin());
-// Hep3Vector p = this._hpstrans.transformVectorToTracking(mcp.getMomentum());
-// // Move to x=0 if needed
-// if(org.x() < 0.) {
-// double dydx = p.y()/p.x();
-// double dzdx = p.z()/p.x();
-// double delta_x = -1. * org.x();
-// double y = delta_x * dydx;
-// double z = delta_x * dzdx;
-// double x = org.x() + delta_x;
-// if( Math.abs(x) > 1e-8) throw new RuntimeException("Error: origin is not zero!");
-// Hep3Vector old = org;
-// org = new BasicHep3Vector(x,y,z);
-// System.out.printf("org %s p %s -> org %s\n", old.toString(),p.toString(),org.toString());
-// } else {
-// org = this._hpstrans.transformVectorToTracking(mcp.getOrigin());
-// }
-//
-//
-//
-// HelixParamCalculator helixParamCalculator = new HelixParamCalculator(p, org, -1*((int)mcp.getCharge()), -1.0*this._B.z());
-// double par[] = new double[5];
-// par[HelicalTrackFit.dcaIndex] = helixParamCalculator.getDCA();
-// par[HelicalTrackFit.slopeIndex] = helixParamCalculator.getSlopeSZPlane();
-// par[HelicalTrackFit.phi0Index] = helixParamCalculator.getPhi0();
-// par[HelicalTrackFit.curvatureIndex] = 1.0/helixParamCalculator.getRadius();
-// par[HelicalTrackFit.z0Index] = helixParamCalculator.getZ0();
-// SymmetricMatrix cov = new SymmetricMatrix(5);
-// for(int i=0;i<cov.getNRows();++i) cov.setElement(i, i, 1.);
-// 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) {
+
+
+
+ private double truthTrackFitChi2(GblUtils.PerigeeParams perPar, GblUtils.PerigeeParams perParTruth, SymmetricMatrix covariance) {
//re-shuffle the param vector to match the covariance order of parameters
BasicMatrix p = new BasicMatrix(1, 5);
p.setElement(0, 0, perPar.getD0());
@@ -1005,190 +844,6 @@
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) {
- if (htf != null) {
- double kappa = -1.0 * Math.signum(B) / htf.R();
- double theta = Math.PI / 2.0 - Math.atan(htf.slope());
- return getPerParVector(kappa, theta, htf.phi0(), htf.dca(), htf.z0());
- }
- return new BasicMatrix(1, 5);
- }
-
- /**
- *
- * Store perigee track parameters.
- *
- * @author Per Hansson Adrian <[log in to unmask]>
- *
- */
- public static class PerigeeParams {
-
- private final BasicMatrix _params;
-
- 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() {
- 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)
- *
- * @param htf input helix to find the track direction
- * @return 3x3 projection matrix
- */
- static Hep3Matrix getPerToClPrj(HelicalTrackFit htf) {
- 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 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));
- trans.setElement(0, 2, VecOp.dot(K, U));
- trans.setElement(1, 0, VecOp.dot(I, V));
- trans.setElement(1, 1, VecOp.dot(J, V));
- trans.setElement(1, 2, VecOp.dot(K, V));
- trans.setElement(2, 0, VecOp.dot(I, T));
- trans.setElement(2, 1, VecOp.dot(J, T));
- trans.setElement(2, 2, VecOp.dot(K, T));
- return trans;
-
- /*
- Hep3Vector B = new BasicHep3Vector(0, 0, 1); // TODO sign convention?
- Hep3Vector H = VecOp.mult(1 / bfield, B);
- Hep3Vector T = HelixUtils.Direction(helix, 0.);
- Hep3Vector HcrossT = VecOp.cross(H, T);
- double alpha = HcrossT.magnitude(); // this should be Bvec cross TrackDir/|B|
- double Q = Math.abs(bfield) * q / p;
- Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
- 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);
- Hep3Vector N = VecOp.mult(1 / alpha, VecOp.cross(H, T)); //-cross(T,H)/alpha = -cross(T,Z) = -J
- double UdotI = VecOp.dot(U, I); // 0,0
- double NdotV = VecOp.dot(N, V); // 1,1?
- double NdotU = VecOp.dot(N, U); // 0,1?
- double TdotI = VecOp.dot(T, I); // 2,0
- double VdotI = VecOp.dot(V, I); // 1,0
- double VdotK = VecOp.dot(V, K); // 1,2
- */
- }
-
-
- /**
- *
- * Store curvilinear track parameters.
- *
- * @author Per Hansson Adrian <[log in to unmask]>
- *
- */
- public static class ClParams {
-
- private BasicMatrix _params = new BasicMatrix(1, 5);
-
- public ClParams(HelicalTrackFit htf, double B) {
-
- if (htf == null) {
- return;
- }
-
- Hep3Matrix perToClPrj = getPerToClPrj(htf);
-
- double d0 = -1 * htf.dca(); //sign convention for curvilinear frame
- double z0 = htf.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();
-
- double lambda = Math.atan(htf.slope());
- double q = Math.signum(htf.R());
- double qOverP = q / htf.p(Math.abs(B));
- double phi = htf.phi0();
-
- _params.setElement(0, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), qOverP);
- _params.setElement(0, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), lambda);
- _params.setElement(0, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), phi);
- _params.setElement(0, FittedGblTrajectory.GBLPARIDX.XT.getValue(), xT);
- _params.setElement(0, FittedGblTrajectory.GBLPARIDX.YT.getValue(), yT);
- }
-
- public BasicMatrix getParams() {
- return _params;
- }
-
- double getQoverP() {
- return _params.e(0, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue());
- }
-
- double getLambda() {
- return _params.e(0, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue());
- }
-
- double getPhi() {
- return _params.e(0, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
- }
-
- double getXt() {
- return _params.e(0, FittedGblTrajectory.GBLPARIDX.XT.getValue());
- }
-
- double getYt() {
- return _params.e(0, FittedGblTrajectory.GBLPARIDX.YT.getValue());
- }
-
- }
-
-
/**
*
* {@link HelicalTrackStripGbl} that explicitly uses the given unit vectors when accessed.
Modified: java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java
=============================================================================
--- java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java (original)
+++ java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java Tue Jan 26 18:37:32 2016
@@ -3,7 +3,7 @@
import hep.physics.vec.BasicHep3Matrix;
import hep.physics.vec.Hep3Matrix;
-import org.hps.recon.tracking.gbl.GBLOutput.PerigeeParams;
+import org.hps.recon.tracking.gbl.GblUtils.PerigeeParams;
import org.lcsim.event.GenericObject;
/**
@@ -77,7 +77,7 @@
/**
* @param perPar is the perigee parameters that is added to object
*/
- public void setPerigeeTrackParameters(PerigeeParams perPar) {
+ public void setPerigeeTrackParameters(GblUtils.PerigeeParams perPar) {
this.bank_double[GBLDOUBLE.PERKAPPA] = perPar.getKappa();
this.bank_double[GBLDOUBLE.PERTHETA] = perPar.getTheta();
this.bank_double[GBLDOUBLE.PERPHI] = perPar.getPhi();
@@ -85,8 +85,8 @@
this.bank_double[GBLDOUBLE.PERZ0] = perPar.getZ0();
}
- public PerigeeParams getPerigeeTrackParameters() {
- return new PerigeeParams(this.bank_double[GBLDOUBLE.PERKAPPA],
+ public GblUtils.PerigeeParams getPerigeeTrackParameters() {
+ return new GblUtils.PerigeeParams(this.bank_double[GBLDOUBLE.PERKAPPA],
this.bank_double[GBLDOUBLE.PERTHETA],
this.bank_double[GBLDOUBLE.PERPHI],
this.bank_double[GBLDOUBLE.PERD0],
Modified: java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GblTrajectory.java
=============================================================================
--- java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GblTrajectory.java (original)
+++ java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GblTrajectory.java Tue Jan 26 18:37:32 2016
@@ -214,7 +214,7 @@
}
/// Retrieve validity of trajectory
- boolean isValid()
+ public boolean isValid()
{
return constructOK;
}
@@ -841,7 +841,7 @@
/**
* \param [in] level print level (0: minimum, >0: more)
*/
- void printPoints(int level)
+ public void printPoints(int level)
{
System.out.println("GblPoints ");
Modified: java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java
=============================================================================
--- java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java (original)
+++ java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java Tue Jan 26 18:37:32 2016
@@ -1,23 +1,234 @@
package org.hps.recon.tracking.gbl;
+import java.util.logging.Logger;
+
import hep.physics.matrix.BasicMatrix;
+import hep.physics.vec.BasicHep3Matrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Matrix;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import org.hps.recon.tracking.HpsHelicalTrackFit;
import org.hps.recon.tracking.MaterialSupervisor;
import org.hps.recon.tracking.MultipleScattering;
+import org.hps.recon.tracking.gbl.matrix.Matrix;
+import org.hps.recon.tracking.gbl.matrix.Vector;
+import org.lcsim.constants.Constants;
import org.lcsim.detector.IDetectorElement;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.fit.helicaltrack.HelixUtils;
import org.lcsim.recon.tracking.seedtracker.ScatterAngle;
/**
- * A class providing various utilities related to GBL
+ * A class with only static utilities related to GBL
*
* @author Per Hansson Adrian <[log in to unmask]>
*
*/
public class GblUtils {
-
+
+ public static Logger LOGGER = Logger.getLogger(GblUtils.class.getName());
+
+
+
+
+ /**
+ * Private constructor to avoid instantiation.
+ */
private GblUtils() {
}
+
+
+ /**
+ *
+ * Store local curvilinear track parameters.
+ *
+ * @author Per Hansson Adrian <[log in to unmask]>
+ *
+ */
+ public static class ClParams {
+
+ private BasicMatrix _params = new BasicMatrix(1, 5);
+
+ public ClParams(HelicalTrackFit htf, double B) {
+
+ if (htf == null) {
+ return;
+ }
+
+ Hep3Matrix perToClPrj = getPerToClPrj(htf);
+
+ double d0 = -1 * htf.dca(); //sign convention for curvilinear frame
+ double z0 = htf.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();
+
+ double lambda = Math.atan(htf.slope());
+ double q = Math.signum(htf.R());
+ double qOverP = q / htf.p(Math.abs(B));
+ double phi = htf.phi0();
+
+ _params.setElement(0, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), qOverP);
+ _params.setElement(0, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), lambda);
+ _params.setElement(0, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), phi);
+ _params.setElement(0, FittedGblTrajectory.GBLPARIDX.XT.getValue(), xT);
+ _params.setElement(0, FittedGblTrajectory.GBLPARIDX.YT.getValue(), yT);
+ }
+
+ public BasicMatrix getParams() {
+ return _params;
+ }
+
+ double getQoverP() {
+ return _params.e(0, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue());
+ }
+
+ double getLambda() {
+ return _params.e(0, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue());
+ }
+
+ double getPhi() {
+ return _params.e(0, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
+ }
+
+ double getXt() {
+ return _params.e(0, FittedGblTrajectory.GBLPARIDX.XT.getValue());
+ }
+
+ double getYt() {
+ return _params.e(0, FittedGblTrajectory.GBLPARIDX.YT.getValue());
+ }
+
+ }
+
+
+
+
+
+ /**
+ *
+ * Store perigee track parameters.
+ *
+ * @author Per Hansson Adrian <[log in to unmask]>
+ *
+ */
+ public static class PerigeeParams {
+
+ private final BasicMatrix _params;
+
+ public PerigeeParams(HelicalTrackFit htf, double B) {
+ _params = GblUtils.getPerParVector(htf, B);
+ }
+
+ public PerigeeParams(double kappa, double theta, double phi, double d0, double z0) {
+ this._params = GblUtils.getPerParVector(kappa, theta, phi, d0, z0);
+ }
+
+ 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);
+ }
+ }
+
+
+
+
+
+
+ /**
+ * Get corrected perigee parameters.
+ * @param locPar - GBL local curvilinear corrections
+ * @param helicalTrackFit - helix
+ * @param bfield - B-field strength
+ * @return corrected parameters
+ */
+ public static double[] getCorrectedPerigeeParameters(Vector locPar, HelicalTrackFit helicalTrackFit, double bfield) {
+
+
+ // Explicitly assign corrections to local variables
+ double qOverPCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.QOVERP.getValue());
+ double xTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
+ double yTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue());
+ double xTCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.XT.getValue());
+ double yTCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.YT.getValue());
+
+
+ // Get helix parameters
+ double qOverP = helicalTrackFit.curvature() / (Constants.fieldConversion * Math.abs(bfield) * Math.sqrt(1 + Math.pow(helicalTrackFit.slope(), 2)));
+ double d0 = -1.0 * helicalTrackFit.dca(); // correct for different sign convention of d0 in perigee frame
+ double z0 = helicalTrackFit.z0();
+ double phi0 = helicalTrackFit.phi0();
+ double lambda = Math.atan(helicalTrackFit.slope());
+
+ // calculate new d0 and z0
+ Hep3Matrix perToClPrj = GblUtils.getPerToClPrj(helicalTrackFit);
+
+ Hep3Matrix clToPerPrj = VecOp.inverse(perToClPrj);
+ Hep3Vector corrPer = VecOp.mult(clToPerPrj, new BasicHep3Vector(xTCorr, yTCorr, 0.0));
+
+ //d0
+ double d0_corr = corrPer.y();
+ double dca_gbl = -1.0 * (d0 + d0_corr);
+
+ //z0
+ double z0_corr = corrPer.z();
+ double z0_gbl = z0 + z0_corr;
+
+ //calculate new slope
+ double lambda_gbl = lambda + yTPrimeCorr;
+ double slope_gbl = Math.tan(lambda_gbl);
+
+ // calculate new curvature
+ double qOverP_gbl = qOverP + qOverPCorr;
+ double C_gbl = Constants.fieldConversion * Math.abs(bfield) * qOverP_gbl / Math.cos(lambda_gbl);
+
+ //calculate new phi0
+ double phi0_gbl = phi0 + xTPrimeCorr - corrPer.x() * C_gbl;
+
+ LOGGER.info("qOverP=" + qOverP + " qOverPCorr=" + qOverPCorr + " qOverP_gbl=" + qOverP_gbl + " ==> pGbl=" + 1.0 / qOverP_gbl + " C_gbl=" + C_gbl);
+
+ LOGGER.info(String.format("corrected helix: d0=%f, z0=%f, omega=%f, tanlambda=%f, phi0=%f, p=%f", dca_gbl, z0_gbl, C_gbl, slope_gbl, phi0_gbl, Math.abs(1 / qOverP_gbl)));
+
+ double parameters_gbl[] = new double[5];
+ parameters_gbl[HelicalTrackFit.dcaIndex] = dca_gbl;
+ parameters_gbl[HelicalTrackFit.phi0Index] = phi0_gbl;
+ parameters_gbl[HelicalTrackFit.curvatureIndex] = C_gbl;
+ parameters_gbl[HelicalTrackFit.z0Index] = z0_gbl;
+ parameters_gbl[HelicalTrackFit.slopeIndex] = slope_gbl;
+
+ return parameters_gbl;
+
+ }
+
+
+
+
public static BasicMatrix gblSimpleJacobianLambdaPhi(double ds, double cosl, double bfac) {
/*
@@ -105,4 +316,172 @@
throw new UnsupportedOperationException("Should not happen. This problem is only solved with the MaterialSupervisor.");
}
}
+
+
+
+
+
+ /**
+ * Calculate the Jacobian from Curvilinear to Perigee frame.
+ * @param helicalTrackFit - original helix
+ * @param helicalTrackFitAtIPCorrected - corrected helix at this point
+ * @param bfield - magnitude of B-field
+ * @return
+ */
+ public static Matrix getCLToPerigeeJacobian(HelicalTrackFit helicalTrackFit, HpsHelicalTrackFit helicalTrackFitAtIPCorrected, double bfield) {
+
+ /*
+ * This part is taken from:
+ // Strandlie, Wittek, NIMA 566, 2006
+ Matrix covariance_gbl = new Matrix(5, 5);
+ //helpers
+ double Bz = -Constants.fieldConversion * Math.abs(bfield); // TODO sign convention and should it be it scaled from Telsa?
+ double p = Math.abs(1 / qOverP_gbl);
+ double q = Math.signum(qOverP_gbl);
+ double tanLambda = Math.tan(lambda_gbl);
+ double cosLambda = Math.cos(lambda_gbl);
+ // Hep3Vector B = new BasicHep3Vector(0, 0, Bz); // TODO sign convention?
+ Hep3Vector H = new BasicHep3Vector(0, 0, 1);
+ Hep3Vector T = HelixUtils.Direction(helix, 0.);
+ Hep3Vector HcrossT = VecOp.cross(H, T);
+ double alpha = HcrossT.magnitude(); // this should be Bvec cross TrackDir/|B|
+ double Q = Bz * q / p;
+ Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
+ 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);
+ Hep3Vector N = VecOp.mult(1 / alpha, VecOp.cross(H, T));
+ double UdotI = VecOp.dot(U, I);
+ double NdotV = VecOp.dot(N, V);
+ double NdotU = VecOp.dot(N, U);
+ double TdotI = VecOp.dot(T, I);
+ double VdotI = VecOp.dot(V, I);
+ double VdotK = VecOp.dot(V, K);
+ covariance_gbl.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), VdotK / TdotI);
+ covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1);
+ covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), -alpha * Q * UdotI * NdotU / (cosLambda * TdotI));
+ covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -alpha * Q * VdotI * NdotU / (cosLambda * TdotI));
+ covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), -1 * Bz / cosLambda);
+ // covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 0);
+ covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), -1 * q * Bz * tanLambda / (p * cosLambda));
+ covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), q * Bz * alpha * Q * tanLambda * UdotI * NdotV / (p * cosLambda * TdotI));
+ covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), q * Bz * alpha * Q * tanLambda * VdotI * NdotV / (p * cosLambda * TdotI));
+ covariance_gbl.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -1 / TdotI);
+ covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), -1);
+ covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), alpha * Q * UdotI * NdotV / TdotI);
+ covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), alpha * Q * VdotI * NdotV / TdotI);
+
+ covariance_gbl.print(15, 13);
+ */
+
+ // Sho's magic below
+
+ // Use projection matrix
+ //TODO should this not be the corrected helix?
+ Hep3Matrix perToClPrj = getPerToClPrj(helicalTrackFit);
+ Hep3Matrix clToPerPrj = VecOp.inverse(perToClPrj);
+ double C_gbl = helicalTrackFitAtIPCorrected.curvature();
+ double lambda_gbl = Math.atan(helicalTrackFitAtIPCorrected.slope());
+ double qOverP_gbl = helicalTrackFitAtIPCorrected.curvature() / (Constants.fieldConversion * Math.abs(bfield) * Math.sqrt(1 + Math.pow(helicalTrackFitAtIPCorrected.slope(), 2)));
+
+ Matrix jacobian = new Matrix(5, 5);
+ jacobian.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), -clToPerPrj.e(1, 0));
+ jacobian.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -clToPerPrj.e(1, 1));
+ jacobian.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1.0);
+ jacobian.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), clToPerPrj.e(0, 1) * C_gbl);
+ jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), Constants.fieldConversion * Math.abs(bfield) / Math.cos(lambda_gbl));
+ jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Constants.fieldConversion * Math.abs(bfield) * qOverP_gbl * Math.tan(lambda_gbl) / Math.cos(lambda_gbl));
+ jacobian.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), clToPerPrj.e(2, 0));
+ jacobian.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), clToPerPrj.e(2, 1));
+ jacobian.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Math.pow(Math.cos(lambda_gbl), -2.0));
+
+ return jacobian;
+ }
+
+
+
+
+
+ /**
+ * Computes the projection matrix from the perigee XY plane variables dca
+ * and z0 into the curvilinear xT,yT,zT frame (U,V,T) with reference point (0,0,0)
+ * for the perigee frame.
+ *
+ * @param htf input helix to find the track direction
+ * @return 3x3 projection matrix
+ */
+ static Hep3Matrix getPerToClPrj(HelicalTrackFit htf) {
+ 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 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));
+ trans.setElement(0, 2, VecOp.dot(K, U));
+ trans.setElement(1, 0, VecOp.dot(I, V));
+ trans.setElement(1, 1, VecOp.dot(J, V));
+ trans.setElement(1, 2, VecOp.dot(K, V));
+ trans.setElement(2, 0, VecOp.dot(I, T));
+ trans.setElement(2, 1, VecOp.dot(J, T));
+ trans.setElement(2, 2, VecOp.dot(K, T));
+ return trans;
+
+ /*
+ Hep3Vector B = new BasicHep3Vector(0, 0, 1); // TODO sign convention?
+ Hep3Vector H = VecOp.mult(1 / bfield, B);
+ Hep3Vector T = HelixUtils.Direction(helix, 0.);
+ Hep3Vector HcrossT = VecOp.cross(H, T);
+ double alpha = HcrossT.magnitude(); // this should be Bvec cross TrackDir/|B|
+ double Q = Math.abs(bfield) * q / p;
+ Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
+ 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);
+ Hep3Vector N = VecOp.mult(1 / alpha, VecOp.cross(H, T)); //-cross(T,H)/alpha = -cross(T,Z) = -J
+ double UdotI = VecOp.dot(U, I); // 0,0
+ double NdotV = VecOp.dot(N, V); // 1,1?
+ double NdotU = VecOp.dot(N, U); // 0,1?
+ double TdotI = VecOp.dot(T, I); // 2,0
+ double VdotI = VecOp.dot(V, I); // 1,0
+ double VdotK = VecOp.dot(V, K); // 1,2
+ */
+ }
+
+
+
+
+
+ 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) {
+ if (htf != null) {
+ double kappa = -1.0 * Math.signum(B) / htf.R();
+ double theta = Math.PI / 2.0 - Math.atan(htf.slope());
+ return getPerParVector(kappa, theta, htf.phi0(), htf.dca(), htf.z0());
+ }
+ return new BasicMatrix(1, 5);
+ }
+
+
}
Modified: java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java
=============================================================================
--- java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java (original)
+++ java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java Tue Jan 26 18:37:32 2016
@@ -14,10 +14,12 @@
import java.util.logging.Formatter;
import java.util.logging.Level;
import java.util.logging.Logger;
+
import org.apache.commons.math3.util.Pair;
-
import org.hps.recon.tracking.TrackUtils;
+
import static org.hps.recon.tracking.gbl.MakeGblTracks.makeCorrectedTrack;
+
import org.hps.recon.tracking.gbl.matrix.Matrix;
import org.hps.recon.tracking.gbl.matrix.SymMatrix;
import org.hps.recon.tracking.gbl.matrix.Vector;
@@ -52,6 +54,8 @@
private final String track2GblTrackRelationName = "TrackToGBLTrack";
private final String gblTrack2StripRelationName = "GBLTrackToStripData";
private final String outputTrackCollectionName = "GBLTracks";
+ private final String trackRelationCollectionName = "MatchedToGBLTrackRelations";
+
private MilleBinary mille;
private String milleBinaryFileName = MilleBinary.DEFAULT_OUTPUT_FILE_NAME;
@@ -184,6 +188,8 @@
LOGGER.info(trackFits.size() + " fitted GBL tracks before adding to event");
List<Track> newTracks = new ArrayList<Track>();
+
+ List<LCRelation> trackRelations = new ArrayList<LCRelation>();
List<GBLKinkData> kinkDataCollection = new ArrayList<GBLKinkData>();
@@ -201,6 +207,10 @@
// Add the track to the list of tracks
newTracks.add(trk.getFirst());
+
+ // Create relation from seed to GBL track
+ trackRelations.add(new BaseLCRelation(fittedTraj.get_seed(), trk.getFirst()));
+
kinkDataCollection.add(trk.getSecond());
kinkDataRelations.add(new BaseLCRelation(trk.getSecond(), trk.getFirst()));
}
@@ -210,6 +220,7 @@
// Put the tracks back into the event and exit
int flag = 1 << LCIOConstants.TRBIT_HITS;
event.put(outputTrackCollectionName, newTracks, Track.class, flag);
+ event.put(trackRelationCollectionName, trackRelations, LCRelation.class, 0);
event.put(GBLKinkData.DATA_COLLECTION, kinkDataCollection, GBLKinkData.class, 0);
event.put(GBLKinkData.DATA_RELATION_COLLECTION, kinkDataRelations, LCRelation.class, 0);
@@ -222,12 +233,15 @@
public static FittedGblTrajectory fit(List<GBLStripClusterData> hits, double bfac, boolean debug) {
// path length along trajectory
double s = 0.;
+ int iLabel;
+
// jacobian to transport errors between points along the path
Matrix jacPointToPoint = new Matrix(5, 5);
jacPointToPoint.UnitMatrix();
// Vector of the strip clusters used for the GBL fit
List<GblPoint> listOfPoints = new ArrayList<GblPoint>();
+ Map<Integer, Double> pathLengthMap = new HashMap<Integer, Double>();
// Store the projection from local to measurement frame for each strip cluster
Map< Integer, Matrix> proL2m_list = new HashMap<Integer, Matrix>();
@@ -236,6 +250,10 @@
//start trajectory at refence point (s=0) - this point has no measurement
GblPoint ref_point = new GblPoint(jacPointToPoint);
listOfPoints.add(ref_point);
+
+ // save path length to each point
+ iLabel = listOfPoints.size();
+ pathLengthMap.put(iLabel, s);
// Loop over strips
int n_strips = hits.size();
@@ -405,7 +423,10 @@
// Add this GBL point to list that will be used in fit
listOfPoints.add(point);
- int iLabel = listOfPoints.size();
+ iLabel = listOfPoints.size();
+
+ // save path length to each point
+ pathLengthMap.put(iLabel, s);
// Update MS covariance matrix
msCov.set(1, 1, msCov.get(1, 1) + scatErr.get(0) * scatErr.get(0));
@@ -493,8 +514,10 @@
LOGGER.fine("locPar " + aCorrection.toString());
-//
- return new FittedGblTrajectory(traj, dVals[0], iVals[0], dVals[1]);
+ FittedGblTrajectory fittedTraj = new FittedGblTrajectory(traj, dVals[0], iVals[0], dVals[1]);
+ fittedTraj.setPathLengthMap(pathLengthMap);
+
+ return fittedTraj;
}
@Override
Modified: java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java
=============================================================================
--- java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java (original)
+++ java/branches/jeremy-dev/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java Tue Jan 26 18:37:32 2016
@@ -5,22 +5,18 @@
import hep.physics.vec.Hep3Matrix;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
+
import java.util.ArrayList;
import java.util.Collection;
-import java.util.Collections;
-import java.util.Comparator;
import java.util.List;
import java.util.logging.Level;
import java.util.logging.Logger;
+
import org.apache.commons.math3.util.Pair;
import org.hps.recon.tracking.CoordinateTransformations;
import org.hps.recon.tracking.MultipleScattering;
import org.hps.recon.tracking.TrackType;
import org.hps.recon.tracking.TrackUtils;
-import static org.hps.recon.tracking.gbl.GBLOutput.getPerToClPrj;
-import org.hps.recon.tracking.gbl.matrix.Matrix;
-import org.hps.recon.tracking.gbl.matrix.SymMatrix;
-import org.hps.recon.tracking.gbl.matrix.Vector;
import org.lcsim.constants.Constants;
import org.lcsim.detector.ITransform3D;
import org.lcsim.detector.tracker.silicon.ChargeCarrier;
@@ -48,7 +44,11 @@
public class MakeGblTracks {
private final static Logger LOGGER = Logger.getLogger(MakeGblTracks.class.getPackage().getName());
-
+ static {
+ LOGGER.setLevel(Level.WARNING);
+ }
+
+
private MakeGblTracks() {
}
@@ -60,7 +60,16 @@
}
}
- public static Pair<Track, GBLKinkData> makeCorrectedTrack(FittedGblTrajectory fittedTraj, HelicalTrackFit helix, List<TrackerHit> trackHits, int trackType, double bfield) {
+ /**
+ * Create a new {@link BaseTrack} from a {@link FittedGblTrajectory}.
+ * @param fittedGblTrajectory
+ * @param helicalTrackFit
+ * @param hitsOnTrack
+ * @param trackType
+ * @param bfield
+ * @return the new {@link BaseTrack} and the kinks along the {@link GblTrajectory} as a {@link Pair}.
+ */
+ public static Pair<Track, GBLKinkData> makeCorrectedTrack(FittedGblTrajectory fittedGblTrajectory, HelicalTrackFit helicalTrackFit, List<TrackerHit> hitsOnTrack, int trackType, double bfield) {
// Initialize the reference point to the origin
double[] ref = new double[]{0., 0., 0.};
@@ -68,221 +77,41 @@
BaseTrack trk = new BaseTrack();
// Add the hits to the track
- for (TrackerHit hit : trackHits) {
+ for (TrackerHit hit : hitsOnTrack) {
trk.addHit(hit);
}
- // Set state at vertex
- Pair<double[], SymmetricMatrix> correctedHelixParams = getGblCorrectedHelixParameters(helix, fittedTraj.get_traj(), bfield, FittedGblTrajectory.GBLPOINT.IP);
+ // Set base track parameters
+ Pair<double[], SymmetricMatrix> correctedHelixParams = fittedGblTrajectory.getCorrectedPerigeeParameters(helicalTrackFit, FittedGblTrajectory.GBLPOINT.IP, bfield);
trk.setTrackParameters(correctedHelixParams.getFirst(), bfield);// hack to set the track charge
trk.getTrackStates().clear();
- TrackState stateVertex = new BaseTrackState(correctedHelixParams.getFirst(), ref, correctedHelixParams.getSecond().asPackedArray(true), TrackState.AtIP, bfield);
- trk.getTrackStates().add(stateVertex);
-
- // Set state at last point
- Pair<double[], SymmetricMatrix> correctedHelixParamsLast = getGblCorrectedHelixParameters(helix, fittedTraj.get_traj(), bfield, FittedGblTrajectory.GBLPOINT.LAST);
+ // Set state at IP
+ TrackState stateIP = new BaseTrackState(correctedHelixParams.getFirst(), ref, correctedHelixParams.getSecond().asPackedArray(true), TrackState.AtIP, bfield);
+ trk.getTrackStates().add(stateIP);
+
+ // Set state at last point on trajectory
+ Pair<double[], SymmetricMatrix> correctedHelixParamsLast = fittedGblTrajectory.getCorrectedPerigeeParameters(helicalTrackFit, FittedGblTrajectory.GBLPOINT.LAST, bfield);
TrackState stateLast = new BaseTrackState(correctedHelixParamsLast.getFirst(), ref, correctedHelixParamsLast.getSecond().asPackedArray(true), TrackState.AtLastHit, bfield);
trk.getTrackStates().add(stateLast);
- GBLKinkData kinkData = getKinks(fittedTraj.get_traj());
+ // Extract kinks from trajectory
+ GBLKinkData kinkData = fittedGblTrajectory.getKinks();
// Set other info needed
- trk.setChisq(fittedTraj.get_chi2());
- trk.setNDF(fittedTraj.get_ndf());
+ trk.setChisq(fittedGblTrajectory.get_chi2());
+ trk.setNDF(fittedGblTrajectory.get_ndf());
trk.setFitSuccess(true);
trk.setRefPointIsDCA(true);
trk.setTrackType(TrackType.setGBL(trackType, true));
// Add the track to the list of tracks
// tracks.add(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: %.10f -> %.10f helix-gbl= %f", i, helix.parameters()[i], trk.getTrackParameter(i), helix.parameters()[i] - trk.getTrackParameter(i)));
- }
- }
+ LOGGER.fine(String.format("helix chi2 %f ndf %d gbl chi2 %f ndf %d\n", helicalTrackFit.chisqtot(), helicalTrackFit.ndf()[0] + helicalTrackFit.ndf()[1], trk.getChi2(), trk.getNDF()));
return new Pair<Track, GBLKinkData>(trk, kinkData);
}
- /**
- * Compute the updated helix parameters and covariance matrix at a given
- * point along the trajectory.
- *
- * @param helix - original seed track
- * @param traj - fitted GBL trajectory
- * @param point - the point along the track where the result is computed.
- * @return corrected parameters.
- */
- public static Pair<double[], SymmetricMatrix> getGblCorrectedHelixParameters(HelicalTrackFit helix, GblTrajectory traj, double bfield, FittedGblTrajectory.GBLPOINT point) {
-
- // get seed helix parameters
- double qOverP = helix.curvature() / (Constants.fieldConversion * Math.abs(bfield) * Math.sqrt(1 + Math.pow(helix.slope(), 2)));
- double d0 = -1.0 * helix.dca(); // correct for different sign convention of d0 in perigee frame
- double z0 = helix.z0();
- double phi0 = helix.phi0();
- double lambda = Math.atan(helix.slope());
-
- LOGGER.info("GblPoint: " + point.toString() + "( " + point.name() + ")");
- LOGGER.info(String.format("original helix: d0=%f, z0=%f, omega=%f, tanlambda=%f, phi0=%f, p=%f", helix.dca(), helix.z0(), helix.curvature(), helix.slope(), helix.phi0(), helix.p(Math.abs(bfield))));
- LOGGER.info("original helix covariance:\n" + helix.covariance());
-
- // get corrections from GBL fit
- Vector locPar = new Vector(5);
- SymMatrix locCov = new SymMatrix(5);
- int pointIndex;
- if (point.compareTo(FittedGblTrajectory.GBLPOINT.IP) == 0) {
- pointIndex = 1;
- } else if (point.compareTo(FittedGblTrajectory.GBLPOINT.LAST) == 0) {
- pointIndex = traj.getNumPoints();
- } else {
- throw new RuntimeException("This GBLPOINT " + point.toString() + "( " + point.name() + ") is not valid");
- }
-
- traj.getResults(pointIndex, locPar, locCov); // vertex point
-// locCov.print(10, 8);
- double qOverPCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.QOVERP.getValue());
- double xTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
- double yTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue());
- double xTCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.XT.getValue());
- double yTCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.YT.getValue());
-// System.out.println((helix.slope() > 0 ? "top: " : "bot ") + "qOverPCorr " + qOverPCorr + " xtPrimeCorr " + xTPrimeCorr + " yTPrimeCorr " + yTPrimeCorr + " xTCorr " + xTCorr + " yTCorr " + yTCorr);
- LOGGER.info((helix.slope() > 0 ? "top: " : "bot ") + "qOverPCorr " + qOverPCorr + " xtPrimeCorr " + xTPrimeCorr + " yTPrimeCorr " + yTPrimeCorr + " xTCorr " + xTCorr + " yTCorr " + yTCorr);
-
- // calculate new d0 and z0
-// Hep3Matrix perToClPrj = traj.get_track_data().getPrjPerToCl();
- Hep3Matrix perToClPrj = getPerToClPrj(helix);
-
- Hep3Matrix clToPerPrj = VecOp.inverse(perToClPrj);
- Hep3Vector corrPer = VecOp.mult(clToPerPrj, new BasicHep3Vector(xTCorr, yTCorr, 0.0));
-
- //d0
- double d0_corr = corrPer.y();
- double dca_gbl = -1.0 * (d0 + d0_corr);
-
- //z0
- double z0_corr = corrPer.z();
- double z0_gbl = z0 + z0_corr;
-
- //calculate new slope
- double lambda_gbl = lambda + yTPrimeCorr;
- double slope_gbl = Math.tan(lambda_gbl);
-
- // calculate new curvature
- double qOverP_gbl = qOverP + qOverPCorr;
-// double pt_gbl = (1.0 / qOverP_gbl) * Math.cos(lambda_gbl);
-// double C_gbl = Constants.fieldConversion * Math.abs(bfield) / pt_gbl;
- double C_gbl = Constants.fieldConversion * Math.abs(bfield) * qOverP_gbl / Math.cos(lambda_gbl);
-
- //calculate new phi0
- double phi0_gbl = phi0 + xTPrimeCorr - corrPer.x() * C_gbl;
-
- LOGGER.info("qOverP=" + qOverP + " qOverPCorr=" + qOverPCorr + " qOverP_gbl=" + qOverP_gbl + " ==> pGbl=" + 1.0 / qOverP_gbl + " C_gbl=" + C_gbl);
-
- LOGGER.info(String.format("corrected helix: d0=%f, z0=%f, omega=%f, tanlambda=%f, phi0=%f, p=%f", dca_gbl, z0_gbl, C_gbl, slope_gbl, phi0_gbl, Math.abs(1 / qOverP_gbl)));
-
- /*
- // Strandlie, Wittek, NIMA 566, 2006
- Matrix covariance_gbl = new Matrix(5, 5);
- //helpers
- double Bz = -Constants.fieldConversion * Math.abs(bfield); // TODO sign convention and should it be it scaled from Telsa?
- double p = Math.abs(1 / qOverP_gbl);
- double q = Math.signum(qOverP_gbl);
- double tanLambda = Math.tan(lambda_gbl);
- double cosLambda = Math.cos(lambda_gbl);
- // Hep3Vector B = new BasicHep3Vector(0, 0, Bz); // TODO sign convention?
- Hep3Vector H = new BasicHep3Vector(0, 0, 1);
- Hep3Vector T = HelixUtils.Direction(helix, 0.);
- Hep3Vector HcrossT = VecOp.cross(H, T);
- double alpha = HcrossT.magnitude(); // this should be Bvec cross TrackDir/|B|
- double Q = Bz * q / p;
- Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
- 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);
- Hep3Vector N = VecOp.mult(1 / alpha, VecOp.cross(H, T));
- double UdotI = VecOp.dot(U, I);
- double NdotV = VecOp.dot(N, V);
- double NdotU = VecOp.dot(N, U);
- double TdotI = VecOp.dot(T, I);
- double VdotI = VecOp.dot(V, I);
- double VdotK = VecOp.dot(V, K);
- covariance_gbl.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), VdotK / TdotI);
- covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1);
- covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), -alpha * Q * UdotI * NdotU / (cosLambda * TdotI));
- covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -alpha * Q * VdotI * NdotU / (cosLambda * TdotI));
- covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), -1 * Bz / cosLambda);
- // covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 0);
- covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), -1 * q * Bz * tanLambda / (p * cosLambda));
- covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), q * Bz * alpha * Q * tanLambda * UdotI * NdotV / (p * cosLambda * TdotI));
- covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), q * Bz * alpha * Q * tanLambda * VdotI * NdotV / (p * cosLambda * TdotI));
- covariance_gbl.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -1 / TdotI);
- covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), -1);
- covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), alpha * Q * UdotI * NdotV / TdotI);
- covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), alpha * Q * VdotI * NdotV / TdotI);
-
- covariance_gbl.print(15, 13);
- */
- // Sho's magic
- Matrix jacobian = new Matrix(5, 5);
- jacobian.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), -clToPerPrj.e(1, 0));
- jacobian.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -clToPerPrj.e(1, 1));
- jacobian.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1.0);
- jacobian.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), clToPerPrj.e(0, 1) * C_gbl);
- jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), Constants.fieldConversion * Math.abs(bfield) / Math.cos(lambda_gbl));
- jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Constants.fieldConversion * Math.abs(bfield) * qOverP_gbl * Math.tan(lambda_gbl) / Math.cos(lambda_gbl));
- jacobian.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), clToPerPrj.e(2, 0));
- jacobian.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), clToPerPrj.e(2, 1));
- jacobian.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Math.pow(Math.cos(lambda_gbl), -2.0));
-
-// jacobian.print(15, 13);
- Matrix helixCovariance = jacobian.times(locCov.times(jacobian.transpose()));
- SymmetricMatrix cov = new SymmetricMatrix(5);
- for (int i = 0; i < 5; i++) {
- for (int j = 0; j < 5; j++) {
- if (i >= j) {
- cov.setElement(i, j, helixCovariance.get(i, j));
- }
- }
- }
- LOGGER.info("corrected helix covariance:\n" + cov);
-
- double parameters_gbl[] = new double[5];
- parameters_gbl[HelicalTrackFit.dcaIndex] = dca_gbl;
- parameters_gbl[HelicalTrackFit.phi0Index] = phi0_gbl;
- parameters_gbl[HelicalTrackFit.curvatureIndex] = C_gbl;
- parameters_gbl[HelicalTrackFit.z0Index] = z0_gbl;
- parameters_gbl[HelicalTrackFit.slopeIndex] = slope_gbl;
-
- return new Pair<double[], SymmetricMatrix>(parameters_gbl, cov);
- }
-
- public static GBLKinkData getKinks(GblTrajectory traj) {
-
- // get corrections from GBL fit
- Vector locPar = new Vector(5);
- SymMatrix locCov = new SymMatrix(5);
- float[] lambdaKinks = new float[traj.getNumPoints() - 1];
- double[] phiKinks = new double[traj.getNumPoints() - 1];
-
- double oldPhi = 0, oldLambda = 0;
- for (int i = 0; i < traj.getNumPoints(); i++) {
- traj.getResults(i + 1, locPar, locCov); // vertex point
- double newPhi = locPar.get(FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
- double newLambda = locPar.get(FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue());
- if (i > 0) {
- lambdaKinks[i - 1] = (float) (newLambda - oldLambda);
- phiKinks[i - 1] = newPhi - oldPhi;
-// System.out.println("phikink: " + (newPhi - oldPhi));
- }
- oldPhi = newPhi;
- oldLambda = newLambda;
- }
-
- return new GBLKinkData(lambdaKinks, phiKinks);
- }
+
/**
* Do a GBL fit to an arbitrary set of strip hits, with a starting value of
@@ -312,6 +141,15 @@
return mergedTrack;
}
+ /**
+ * Do a GBL fit to a list of {@link TrackerHit}.
+ * @param htf - seed fit
+ * @param stripHits - list of {@link TrackerHit}.
+ * @param _scattering - estimation of the multiple scattering {@link MultipleScattering}.
+ * @param bfield - magnitude of B-field.
+ * @param debug - debug flag.
+ * @return
+ */
public static FittedGblTrajectory doGBLFit(HelicalTrackFit htf, List<TrackerHit> stripHits, MultipleScattering _scattering, double bfield, int debug) {
List<GBLStripClusterData> stripData = makeStripData(htf, stripHits, _scattering, bfield, debug);
double bfac = Constants.fieldConversion * bfield;
@@ -320,6 +158,15 @@
return fit;
}
+ /**
+ * Create a list of {@link GBLStripClusterData} objects that can be used as input to the GBL fitter.
+ * @param htf
+ * @param stripHits
+ * @param _scattering
+ * @param _B
+ * @param _debug
+ * @return
+ */
public static List<GBLStripClusterData> makeStripData(HelicalTrackFit htf, List<TrackerHit> stripHits, MultipleScattering _scattering, double _B, int _debug) {
List<GBLStripClusterData> stripClusterDataList = new ArrayList<GBLStripClusterData>();
Modified: java/branches/jeremy-dev/tracking/src/main/java/org/hps/svt/alignment/BuildMillepedeCompact.java
=============================================================================
--- java/branches/jeremy-dev/tracking/src/main/java/org/hps/svt/alignment/BuildMillepedeCompact.java (original)
+++ java/branches/jeremy-dev/tracking/src/main/java/org/hps/svt/alignment/BuildMillepedeCompact.java Tue Jan 26 18:37:32 2016
@@ -56,6 +56,7 @@
private static boolean replaceConstant = false;
private static boolean calcNewValue = true;
private static boolean ignoreBeamspot = true;
+ private static boolean flipRotationCorrection = false;
@@ -65,6 +66,7 @@
options.addOption(new Option("o",true,"The name of the new compact xml file."));
options.addOption(new Option("r", false, "Replace correction instead of adding to it."));
options.addOption(new Option("t", false, "Add a text string as a new value instead of adding to it."));
+ options.addOption(new Option("f",false, "Flip sign of rotation corrections."));
return options;
}
@@ -113,6 +115,23 @@
calcNewValue = false;
}
+ if(cl.hasOption("f")) {
+ flipRotationCorrection = true;
+ }
+
+ if (calcNewValue)
+ System.out.println("DO CALCULATE NEW VALUE");
+ else
+ System.out.println("DO NOT CALCULATE NEW VALUE");
+
+
+
+ if (flipRotationCorrection)
+ System.out.println("DO FLIP ROTATIONS");
+ else
+ System.out.println("DO NOT FLIP ROTATIONS");
+
+
File compactFile = new File(compactFilename);
@@ -197,10 +216,12 @@
if(replaceConstant) {
newValue = correction;
} else {
- if (p.getType() == MilleParameter.Type.ROTATION.getType()) {
+ if (p.getType() == MilleParameter.Type.ROTATION.getType() && !flipRotationCorrection) {
newValue = oldValue - correction;
+ System.out.println("NOFLIP");
} else {
newValue = oldValue + correction;
+ System.out.println("FLIP");
}
}
System.out.println("Update " + p.getId() + ": " + oldValue + " (corr. " + correction + ") -> " + newValue );
@@ -214,7 +235,16 @@
throw new RuntimeException("Doesn't make sense to try and replace with the string option?");
if( correction != 0.0) {
- String newValue = oldValue + " + " + String.format("%.6f",correction);
+ String newValue;
+
+ if (p.getType() == MilleParameter.Type.ROTATION.getType() && !flipRotationCorrection) {
+ newValue = oldValue + " - " + String.format("%.6f",correction);
+ System.out.println("NOFLIP");
+ } else {
+ newValue = oldValue + " + " + String.format("%.6f",correction);
+ System.out.println("FLIP");
+ }
+
System.out.println("Update " + p.getId() + ": " + oldValue + " (corr. " + correction + ") -> " + newValue );
node.setAttribute("value", newValue);
}
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/meeg/SvtChargeIntegrator.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/meeg/SvtChargeIntegrator.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/meeg/SvtChargeIntegrator.java Tue Jan 26 18:37:32 2016
@@ -25,6 +25,7 @@
import org.hps.conditions.svt.SvtBiasConstant.SvtBiasConstantCollection;
import org.hps.conditions.svt.SvtMotorPosition;
import org.hps.conditions.svt.SvtMotorPosition.SvtMotorPositionCollection;
+import org.hps.conditions.svt.SvtTimingConstants;
import org.hps.run.database.RunManager;
/**
@@ -39,15 +40,9 @@
private static Logger LOGGER = Logger.getLogger(SvtChargeIntegrator.class.getPackage().getName());
private static final double angleTolerance = 1e-4;
+ private static final double burstModeNoiseEfficiency = 0.965;
/**
- * Default constructor
- */
- public SvtChargeIntegrator() {
- }
-
- /**
- * Load SVT HV bias constants into the conditions database.
*
* @param args the command line arguments (requires a CVS run log file and a
* MYA dump file.)
@@ -115,14 +110,18 @@
BufferedReader br = new BufferedReader(new FileReader(cl.getArgs()[1]));
String line;
System.err.println("myaData header: " + br.readLine()); //discard the first line
-// System.out.println("run\ttotalQ\ttotalQBias\tfracBias\ttotalQNom\tfracNom\ttotalQ1pt5\tfrac1pt5\ttotalGatedQ\ttotalGatedQBias\tfracGatedBias\ttotalGatedQNom\tfracGatedNom\ttotalGatedQ1pt5\tfracGated1pt5");
- System.out.println("run_num\tfile_num\tnominal_position\tnEvents\ttotalQ\ttotalQ_withbias\ttotalQ_atnom\tgatedQ\tgatedQ_withbias\tgatedQ_atnom");
+ if (perRun) {
+ System.out.println("run_num\tnominal_position\tnEvents\ttotalQ\ttotalQ_withbias\ttotalQ_atnom\tgatedQ\tgatedQ_withbias\tgatedQ_atnom\tgoodQ\tgoodQ_withbias\tgoodQ_atnom");
+ } else {
+ System.out.println("run_num\tfile_num\tnominal_position\tnEvents\ttotalQ\ttotalQ_withbias\ttotalQ_atnom\tgatedQ\tgatedQ_withbias\tgatedQ_atnom\tgoodQ\tgoodQ_withbias\tgoodQ_atnom");
+ }
int currentRun = 0;
double nominalAngleTop = -999;
double nominalAngleBottom = -999;
String nominalPosition = null;
long tiTimeOffset = 0;
+ double efficiency = 0;
SvtBiasConstantCollection svtBiasConstants = null;
SvtMotorPositionCollection svtPositionConstants = null;
SvtAlignmentConstant.SvtAlignmentConstantCollection alignmentConstants = null;
@@ -170,6 +169,20 @@
} else {
nominalPosition = "unknown";
}
+ efficiency = burstModeNoiseEfficiency;
+ SvtTimingConstants svtTimingConstants;
+ try {
+ svtTimingConstants = DatabaseConditionsManager.getInstance().getCachedConditions(SvtTimingConstants.SvtTimingConstantsCollection.class, "svt_timing_constants").getCachedData().get(0);
+ } catch (Exception ex) {
+ svtTimingConstants = null;
+ }
+ if (svtTimingConstants != null) {
+ if (svtTimingConstants.getOffsetTime() > 27) {
+ efficiency *= 2.0 / 3.0; // bad latency: drop 2 out of 6 trigger phases
+ }// otherwise, we have good latency
+ } else {
+ efficiency = 0;
+ }//no latency info in conditions: give up
currentRun = runNum;
}
@@ -182,6 +195,9 @@
double totalGatedCharge = 0;
double totalGatedChargeWithBias = 0;
double totalGatedChargeWithBiasAtNominal = 0;
+ double totalGoodCharge = 0;
+ double totalGoodChargeWithBias = 0;
+ double totalGoodChargeWithBiasAtNominal = 0;
br.mark(1000);
while ((line = br.readLine()) != null) {
@@ -227,17 +243,18 @@
if (lastDate != null) {
double dt = (Math.min(date.getTime(), endDate.getTime()) - Math.max(startDate.getTime(), lastDate.getTime())) / 1000.0;
- double dq = dt * current; // nC
- double dqGated = dt * current * livetime; // nC
// System.out.format("start %d end %d date %d lastDate %d current %f dt %f\n", startDate.getTime(), endDate.getTime(), date.getTime(), lastDate.getTime(), current, dt);
- totalCharge += dq;
- totalGatedCharge += dqGated;
+ totalCharge += dt * current; // nC
+ totalGatedCharge += dt * current * livetime;
+ totalGoodCharge += dt * current * livetime * efficiency;
if (biasGood) {
totalChargeWithBias += dq;
totalGatedChargeWithBias += dqGated;
+ totalGoodChargeWithBias += biasDt * current * livetime * efficiency;
if (positionGood) {
totalChargeWithBiasAtNominal += dq;
totalGatedChargeWithBiasAtNominal += dqGated;
+ totalGoodChargeWithBiasAtNominal += positionDt * current * livetime * efficiency;
}
}
}
@@ -248,8 +265,14 @@
}
br.mark(1000);
}
- int nEvents = Integer.parseInt(record.get(2));
- System.out.format("%d\t%d\t%s\t%d\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\n", runNum, fileNum, nominalPosition, nEvents, totalCharge, totalChargeWithBias, totalChargeWithBiasAtNominal, totalGatedCharge, totalGatedChargeWithBias, totalGatedChargeWithBiasAtNominal);
+ if (perRun) {
+ int nEvents = Integer.parseInt(record.get(9));
+ System.out.format("%d\t%s\t%d\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\n", runNum, nominalPosition, nEvents, totalCharge, totalChargeWithBias, totalChargeWithBiasAtNominal, totalGatedCharge, totalGatedChargeWithBias, totalGatedChargeWithBiasAtNominal, totalGoodCharge, totalGoodChargeWithBias, totalGoodChargeWithBiasAtNominal);
+ } else {
+ int fileNum = Integer.parseInt(record.get(1));
+ int nEvents = Integer.parseInt(record.get(2));
+ System.out.format("%d\t%d\t%s\t%d\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\n", runNum, fileNum, nominalPosition, nEvents, totalCharge, totalChargeWithBias, totalChargeWithBiasAtNominal, totalGatedCharge, totalGatedChargeWithBias, totalGatedChargeWithBiasAtNominal, totalGoodCharge, totalGoodChargeWithBias, totalGoodChargeWithBiasAtNominal);
+ }
}
} catch (Exception ex) {
Logger.getLogger(SvtChargeIntegrator.class.getName()).log(Level.SEVERE, null, ex);
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/ExamplePlotter.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/ExamplePlotter.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/ExamplePlotter.java Tue Jan 26 18:37:32 2016
@@ -4,6 +4,7 @@
import hep.aida.IHistogram1D;
import hep.aida.IPlotter;
import hep.aida.IPlotterStyle;
+import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import java.io.IOException;
@@ -12,9 +13,10 @@
import java.util.logging.Logger;
import org.hps.recon.tracking.BeamlineConstants;
-import org.hps.recon.tracking.HPSTrack;
+import org.hps.recon.tracking.HpsHelicalTrackFit;
import org.hps.recon.tracking.HelixConverter;
import org.hps.recon.tracking.StraightLineTrack;
+import org.hps.recon.tracking.TrackUtils;
import org.lcsim.event.EventHeader;
import org.lcsim.event.Track;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
@@ -76,15 +78,9 @@
aida.histogram1D("Track Momentum (Pz)").fill(trk.getTrackStates().get(0).getMomentum()[0]);
aida.histogram1D("Track Chi2").fill(trk.getChi2());
- SeedTrack stEle = (SeedTrack) trk;
- SeedCandidate seedEle = stEle.getSeedCandidate();
- HelicalTrackFit ht = seedEle.getHelix();
- HelixConverter converter = new HelixConverter(0);
- StraightLineTrack slt = converter.Convert(ht);
- HPSTrack hpstrack = new HPSTrack(ht);
- Hep3Vector[] trkatconver = hpstrack.getPositionAtZMap(100, BeamlineConstants.HARP_POSITION_TESTRUN, 1);
- aida.histogram1D("X (mm) @ Converter").fill(trkatconver[0].x()); // y tracker frame?
- aida.histogram1D("Y (mm) @ Converter").fill(trkatconver[0].y()); // z tracker frame?
+ Hep3Vector trkatconver = new BasicHep3Vector(TrackUtils.extrapolateTrackUsingFieldMap(trk, 100.0, BeamlineConstants.HARP_POSITION_TESTRUN, 5.0, event.getDetector().getFieldMap()).getReferencePoint());
+ aida.histogram1D("X (mm) @ Converter").fill(trkatconver.x()); // y tracker frame?
+ aida.histogram1D("Y (mm) @ Converter").fill(trkatconver.y()); // z tracker frame?
}
}
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/HelicalTrackHitResidualsDriver.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/HelicalTrackHitResidualsDriver.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/HelicalTrackHitResidualsDriver.java Tue Jan 26 18:37:32 2016
@@ -16,10 +16,11 @@
import java.util.logging.Logger;
+
//===> import org.hps.conditions.deprecated.SvtUtils;
import org.hps.recon.tracking.EventQuality;
import org.hps.recon.tracking.TrackUtils;
-import org.hps.users.phansson.TrigRateDriver;
+import org.hps.users.phansson.testrun.TrigRateDriver;
import org.lcsim.detector.tracker.silicon.HpsSiSensor;
import org.lcsim.event.EventHeader;
import org.lcsim.event.RawTrackerHit;
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/TrackExtrapolationAnalysis.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/TrackExtrapolationAnalysis.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/TrackExtrapolationAnalysis.java Tue Jan 26 18:37:32 2016
@@ -5,13 +5,14 @@
import hep.aida.IHistogram2D;
import hep.aida.IPlotter;
import hep.aida.IPlotterStyle;
+import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import java.util.List;
import java.util.Map;
import org.hps.recon.tracking.BeamlineConstants;
-import org.hps.recon.tracking.HPSTrack;
+import org.hps.recon.tracking.HpsHelicalTrackFit;
import org.hps.recon.tracking.HelixConverter;
import org.hps.recon.tracking.StraightLineTrack;
import org.hps.recon.tracking.TrackUtils;
@@ -99,10 +100,9 @@
charge = 0;//make plot look pretty
// System.out.println("Charge = " + charge + "; isTop = " + isTop);
- HPSTrack hpstrk=null;
- hpstrk = new HPSTrack(ht);
-// Hep3Vector posAtConv = hpstrk.getPositionAtZ(zAtConverter, -101, -100, 0.1);
- Hep3Vector posAtConv = hpstrk.getPositionAtZMap(100,BeamlineConstants.HARP_POSITION_TESTRUN , 5.0)[0];
+ HpsHelicalTrackFit hpstrk=null;
+ hpstrk = new HpsHelicalTrackFit(ht);
+ Hep3Vector posAtConv = new BasicHep3Vector(TrackUtils.extrapolateTrackUsingFieldMap(trk, 100, BeamlineConstants.HARP_POSITION_TESTRUN, 5.0, event.getDetector().getFieldMap()).getReferencePoint());
double useThisx=posAtConv.x();
double useThisy=posAtConv.y();
@@ -123,7 +123,7 @@
aida.histogram1D("Negative Y (mm) @ Converter").fill(useThisy);
}
// Hep3Vector posAtConvShort = hpstrk.getPositionAtZ(zAtConverter, -0.1, 0, 0.01);
- Hep3Vector posAtConvShort = hpstrk.getPositionAtZMap(0,BeamlineConstants.HARP_POSITION_TESTRUN, 5.0)[0];
+ Hep3Vector posAtConvShort = new BasicHep3Vector(TrackUtils.extrapolateTrackUsingFieldMap(trk, 0, BeamlineConstants.HARP_POSITION_TESTRUN, 5.0, event.getDetector().getFieldMap()).getReferencePoint());
aida.histogram2D("Extrapolated X: short vs long fringe").fill(posAtConvShort.x(), posAtConv.x());
aida.histogram2D("Extrapolated Y: short vs long fringe").fill(posAtConvShort.y(), posAtConv.y());
@@ -144,7 +144,8 @@
// Hep3Vector posAtEcalHPS = hpstrk.getPositionAtZMap(750,zCluster, 5.0);
double zCluster=clust.getPosition()[2];
// double zCluster=1450.0;
- Hep3Vector posAtEcalHPS = hpstrk.getPositionAtZMap(750,zCluster, 5.0)[0];
+ Hep3Vector posAtEcalHPS = new BasicHep3Vector(TrackUtils.extrapolateTrackUsingFieldMap(trk, 750, zCluster, 5.0, event.getDetector().getFieldMap()).getReferencePoint());
+
Hep3Vector posAtEcalExtend= TrackUtils.extrapolateTrack(trk,zCluster);
aida.histogram2D("ECal Extrapolation X : HPS vs Extend").fill( posAtEcalExtend.y(),posAtEcalHPS.x()-posAtEcalExtend.y());
aida.histogram2D("ECal Extrapolation Y : HPS vs Extend").fill( posAtEcalExtend.z(),posAtEcalHPS.y()-posAtEcalExtend.z());
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/TwoTrackAnalysis.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/TwoTrackAnalysis.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/mgraham/TwoTrackAnalysis.java Tue Jan 26 18:37:32 2016
@@ -17,10 +17,12 @@
import java.util.logging.Level;
import java.util.logging.Logger;
+import org.apache.commons.lang3.NotImplementedException;
import org.hps.recon.tracking.BeamlineConstants;
-import org.hps.recon.tracking.HPSTrack;
+import org.hps.recon.tracking.HpsHelicalTrackFit;
import org.hps.recon.tracking.HelixConverter;
import org.hps.recon.tracking.StraightLineTrack;
+import org.hps.recon.tracking.TrackUtils;
import org.lcsim.event.EventHeader;
import org.lcsim.event.Track;
import org.lcsim.event.TrackerHit;
@@ -273,10 +275,11 @@
HelicalTrackFit ht = seedEle.getHelix();
HelixConverter converter = new HelixConverter(0);
StraightLineTrack slt = converter.Convert(ht);
- HPSTrack hpstrack = new HPSTrack(ht);
- Hep3Vector[] trkatconver = hpstrack.getPositionAtZMap(100, BeamlineConstants.HARP_POSITION_TESTRUN, 1);
- aida.histogram1D("X (mm) @ Converter using Map").fill(trkatconver[0].x()); // y tracker frame?
- aida.histogram1D("Y (mm) @ Converter using Map").fill(trkatconver[0].y()); // z tracker frame?
+ HpsHelicalTrackFit hpstrack = new HpsHelicalTrackFit(ht);
+ Hep3Vector trkatconver = new BasicHep3Vector(TrackUtils.extrapolateTrackUsingFieldMap(trk, 100, BeamlineConstants.HARP_POSITION_TESTRUN, 1, event.getDetector().getFieldMap()).getReferencePoint());
+
+ aida.histogram1D("X (mm) @ Converter using Map").fill(trkatconver.x()); // y tracker frame?
+ aida.histogram1D("Y (mm) @ Converter using Map").fill(trkatconver.y()); // z tracker frame?
if (slt != null) {
aida.histogram1D("X (mm) @ Converter using SLT").fill(slt.getYZAtX(BeamlineConstants.HARP_POSITION_TESTRUN)[0]); // y tracker frame?
aida.histogram1D("Y (mm) @ Converter using SLT").fill(slt.getYZAtX(BeamlineConstants.HARP_POSITION_TESTRUN)[1]); // z tracker frame?
@@ -320,9 +323,9 @@
// HPSTrack hpstrack2 = new HPSTrack(ht2);
// Hep3Vector[] trkatconver2 = hpstrack2.getPositionAtZMap(100, BeamlineConstants.HARP_POSITION, 1);
- HPSTrack hpstrack1 = new HPSTrack(ht1);
+ HpsHelicalTrackFit hpstrack1 = new HpsHelicalTrackFit(ht1);
Hep3Vector[] trkatconver1 = {new BasicHep3Vector(), new BasicHep3Vector(0, 0, 0)};
- HPSTrack hpstrack2 = new HPSTrack(ht2);
+ HpsHelicalTrackFit hpstrack2 = new HpsHelicalTrackFit(ht2);
Hep3Vector[] trkatconver2 = {new BasicHep3Vector(), new BasicHep3Vector(0, 0, 0)};;
if (isMC) {
double[] t1 = slt1.getYZAtX(BeamlineConstants.HARP_POSITION_TESTRUN);
@@ -330,8 +333,9 @@
trkatconver1[0] = new BasicHep3Vector(t1[0], t1[1], BeamlineConstants.HARP_POSITION_TESTRUN);
trkatconver2[0] = new BasicHep3Vector(t2[0], t2[1], BeamlineConstants.HARP_POSITION_TESTRUN);
} else {
- trkatconver1 = hpstrack1.getPositionAtZMap(100, BeamlineConstants.HARP_POSITION_TESTRUN, 1);
- trkatconver2 = hpstrack2.getPositionAtZMap(100, BeamlineConstants.HARP_POSITION_TESTRUN, 1);
+ throw new NotImplementedException("Need to implement using TrackUtils to extrapolate tracks!");
+ //trkatconver1 = hpstrack1.getPositionAtZMap(100, BeamlineConstants.HARP_POSITION_TESTRUN, 1);
+ //trkatconver2 = hpstrack2.getPositionAtZMap(100, BeamlineConstants.HARP_POSITION_TESTRUN, 1);
}
List<TrackerHit> hitsOnTrack1 = trk1.getTrackerHits();
int layer1;
@@ -508,9 +512,10 @@
posvec[1] = slt1.getYZAtX(z)[1];
posvec[2] = z;
} else {
- Hep3Vector[] trk1atz = hpstrack1.getPositionAtZMap(100, z, 1);
- posvec[0] = trk1atz[0].x();
- posvec[1] = trk1atz[0].y();
+ Hep3Vector trk1atz = new BasicHep3Vector(TrackUtils.extrapolateTrackUsingFieldMap(trk1, 100, z, 1, event.getDetector().getFieldMap()).getReferencePoint());
+
+ posvec[0] = trk1atz.x();
+ posvec[1] = trk1atz.y();
posvec[2] = z;
}
Trk1.add(posvec);
@@ -592,9 +597,10 @@
posvec2[1] = slt2.getYZAtX(z)[1];
posvec2[2] = z;
} else {
- Hep3Vector[] trk2atz = hpstrack2.getPositionAtZMap(100, z, 1);
- posvec2[0] = trk2atz[0].x();
- posvec2[1] = trk2atz[0].y();
+ Hep3Vector trk2atz = new BasicHep3Vector(TrackUtils.extrapolateTrackUsingFieldMap(trk2, 100, z, 1, event.getDetector().getFieldMap()).getReferencePoint());
+
+ posvec2[0] = trk2atz.x();
+ posvec2[1] = trk2atz.y();
posvec2[2] = z;
}
Trk2.add(posvec2);
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/CmpGenToFittedTracksDriver.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/CmpGenToFittedTracksDriver.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/CmpGenToFittedTracksDriver.java Tue Jan 26 18:37:32 2016
@@ -12,6 +12,7 @@
import java.util.logging.Level;
import java.util.logging.Logger;
+import org.hps.users.phansson.testrun.TrigRateDriver;
import org.lcsim.event.EventHeader;
import org.lcsim.event.Track;
import org.lcsim.event.TrackerHit;
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/DataTrackerFakeHitDriver.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/DataTrackerFakeHitDriver.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/DataTrackerFakeHitDriver.java Tue Jan 26 18:37:32 2016
@@ -16,10 +16,12 @@
import java.util.List;
import java.util.Set;
+
//===> import org.hps.conditions.deprecated.SvtUtils;
-import org.hps.recon.tracking.HPSTrack;
+import org.hps.recon.tracking.HpsHelicalTrackFit;
import org.hps.recon.tracking.TrackUtils;
import org.hps.recon.tracking.TrackerHitUtils;
+import org.hps.recon.tracking.WTrack;
import org.lcsim.detector.IDetectorElement;
import org.lcsim.detector.ITransform3D;
import org.lcsim.detector.ITranslation3D;
@@ -194,12 +196,12 @@
// Obtain the tracks from the event
- if (!event.hasCollection(HPSTrack.class, trackCollectionName)) {
+ if (!event.hasCollection(HpsHelicalTrackFit.class, trackCollectionName)) {
this.printDebug("No HPSTracks were found, skipping event");
simHits = null;
return;
}
- List<HPSTrack> tracks = event.get(HPSTrack.class, trackCollectionName);
+ List<HpsHelicalTrackFit> tracks = event.get(HpsHelicalTrackFit.class, trackCollectionName);
if (debug) {
System.out.println(this.getClass().getSimpleName() + ": found " + tracks.size() + " tracks (" + this.trackCollectionName + ")");
@@ -223,13 +225,13 @@
System.out.println(this.getClass().getSimpleName() + ": Add hits for " + tracks.size() + " tracks (" + this.trackCollectionName + ")");
}
- for (HPSTrack helix : tracks) {
+ for (HpsHelicalTrackFit helix : tracks) {
if (debug) {
System.out.println(this.getClass().getSimpleName() + ": trying to add hits for this track");
}
// Get the MC Particle associated with this track
- MCParticle mcParticle = helix.getMCParticle();
+ MCParticle mcParticle = helix.getMcParticle();
if (debug) {
System.out.println(this.getClass().getSimpleName() + helix.toString());
@@ -238,8 +240,8 @@
System.out.println(this.getClass().getSimpleName() + ": create a WTrack object");
}
- WTrack wtrack = new WTrack(helix, Math.abs(_bfield.z()), true); //remove sign from B-field (assumed to go along z-direction)
-
+ WTrack wtrack = new WTrack(helix, Math.abs(_bfield.z()));
+
if (debug) {
System.out.println(this.getClass().getSimpleName() + ": " + wtrack.toString());
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/ParticleHelixProducer.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/ParticleHelixProducer.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/ParticleHelixProducer.java Tue Jan 26 18:37:32 2016
@@ -13,7 +13,7 @@
import java.util.List;
import org.hps.analysis.ecal.HPSMCParticlePlotsDriver;
-import org.hps.recon.tracking.HPSTrack;
+import org.hps.recon.tracking.HpsHelicalTrackFit;
import org.hps.recon.tracking.TrackerHitUtils;
import org.lcsim.constants.Constants;
import org.lcsim.event.EventHeader;
@@ -137,7 +137,7 @@
//Make new tracks based on the MC particles
//List<HelicalTrackFit> tracks = new ArrayList<HelicalTrackFit>();
- List<HPSTrack> tracks = new ArrayList<HPSTrack>();
+ List<HpsHelicalTrackFit> tracks = new ArrayList<HpsHelicalTrackFit>();
if (event.hasCollection(MCParticle.class)) {
List<MCParticle> mcparticles = event.get(MCParticle.class).get(0);
@@ -251,7 +251,7 @@
pars[3] = hpc.getZ0();
pars[4] = hpc.getSlopeSZPlane();
//HelicalTrackFit htf = this.trackUtils.makeHelicalTrackFit(pars);
- HPSTrack htf = this.makeHPSTrack(pars, part);
+ HpsHelicalTrackFit htf = this.makeHPSTrack(pars, part);
tracks.add(htf);
if (debug) {
System.out.println(this.getClass().getSimpleName() + ": MC particle created HelicalTrackFit " + htf.toString());
@@ -281,7 +281,7 @@
}
this.printDebug("created " + tracks.size() + " MC particle helix tracks");
- event.put(this.trackOutputCollectionName, tracks, HPSTrack.class, 0);
+ event.put(this.trackOutputCollectionName, tracks, HpsHelicalTrackFit.class, 0);
_totalTracks += tracks.size();
}
@@ -292,9 +292,12 @@
* @param mcParticle : MC particle associated to this HelicalTrackFit
* @return HpsHelicalTrackFit :
*/
- public HPSTrack makeHPSTrack(double[] helixParameters, MCParticle mcParticle) {
- return new HPSTrack(helixParameters, new SymmetricMatrix(5), new double[2], new int[2],
- new HashMap<HelicalTrackHit, Double>(), new HashMap<HelicalTrackHit, MultipleScatter>(), mcParticle);
+ public HpsHelicalTrackFit makeHPSTrack(double[] helixParameters, MCParticle mcParticle) {
+ HpsHelicalTrackFit helicalTrackFit = new HpsHelicalTrackFit(helixParameters, new SymmetricMatrix(5), new double[2], new int[2],
+ new HashMap<HelicalTrackHit, Double>(), new HashMap<HelicalTrackHit, MultipleScatter>());
+ helicalTrackFit.setMcParticle(mcParticle);
+ return helicalTrackFit;
+
}
/**
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TrackingReconstructionPlots.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TrackingReconstructionPlots.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TrackingReconstructionPlots.java Tue Jan 26 18:37:32 2016
@@ -50,7 +50,8 @@
import org.lcsim.util.aida.AIDA;
/**
- *
+ * Analysis class to check recon.
+ *
* @author phansson
*/
public class TrackingReconstructionPlots extends Driver {
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TriggerTurnOnAnalysis.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TriggerTurnOnAnalysis.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TriggerTurnOnAnalysis.java Tue Jan 26 18:37:32 2016
@@ -10,12 +10,15 @@
import hep.aida.IHistogramFactory;
import hep.aida.IPlotter;
import hep.physics.vec.Hep3Vector;
+
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;
import java.util.logging.Level;
import java.util.logging.Logger;
+
import org.hps.analysis.ecal.HPSMCParticlePlotsDriver;
+import org.hps.users.phansson.testrun.TrigRateDriver;
import org.lcsim.event.Cluster;
import org.lcsim.event.EventHeader;
import org.lcsim.event.GenericObject;
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TruthMomentumResolutionDriver.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TruthMomentumResolutionDriver.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/phansson/TruthMomentumResolutionDriver.java Tue Jan 26 18:37:32 2016
@@ -18,6 +18,7 @@
import java.util.logging.Logger;
import org.hps.analysis.ecal.HPSMCParticlePlotsDriver;
+import org.hps.users.phansson.testrun.TrigRateDriver;
import org.lcsim.constants.Constants;
import org.lcsim.event.EventHeader;
import org.lcsim.event.MCParticle;
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/StyleUtil.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/StyleUtil.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/StyleUtil.java Tue Jan 26 18:37:32 2016
@@ -1,10 +1,9 @@
package org.hps.users.spaul;
import java.util.Arrays;
+import java.util.Random;
-import hep.aida.IAnalysisFactory;
-import hep.aida.IFunction;
-import hep.aida.IPlotterRegion;
+import hep.aida.*;
public class StyleUtil {
@@ -27,18 +26,20 @@
r.style().yAxisStyle().setLabel(ly);
r.style().yAxisStyle().labelStyle().setFontSize(16);
r.style().yAxisStyle().tickLabelStyle().setFontSize(14);
- r.style().statisticsBoxStyle().setParameter("backgroundColor", "white");
- r.style().statisticsBoxStyle().boxStyle().backgroundStyle().setColor("White");
- r.style().statisticsBoxStyle().boxStyle().backgroundStyle().setParameter("color", "white");
- r.style().statisticsBoxStyle().boxStyle().backgroundStyle().setOpacity(100);
- System.out.println(Arrays.toString(r.style().dataStyle().availableParameters()));
+ // r.style().statisticsBoxStyle().set;
+ //System.out.println(Arrays.toString());
+ r.style().legendBoxStyle().textStyle().setFontSize(16);
+
+ //r.style().dataStyle().showInLegendBox(false);
+ r.style().legendBoxStyle().boxStyle().foregroundStyle().setOpacity(1.0);
r.style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
r.style().dataStyle().fillStyle().setParameter("showZeroHeightBins", "false");
+ System.out.println(Arrays.toString(r.style().statisticsBoxStyle().boxStyle().backgroundStyle().availableParameters()));
//r.style().dataStyle().setParameter("showDataInStatisticsBox", "false");
r.style().setParameter("hist2DStyle", "colorMap");
//r.style().dataBoxStyle()
}
- static void addHorizontalLine(IPlotterRegion r, double y, String name){
+ public static void addHorizontalLine(IPlotterRegion r, double y, String name){
IAnalysisFactory af = IAnalysisFactory.create();
IFunction f = af.createFunctionFactory(af.createTreeFactory().create()).createFunctionByName(name, "p0");
f.setParameter("p0", y);
@@ -53,4 +54,54 @@
f.setParameter("p2", p2);
region.plot(f);
}
+ public static void noFillHistogramBars(IPlotterRegion region) {
+ region.style().dataStyle().setParameter("fillHistogramBars", "false");
+ region.style().dataStyle().fillStyle().setVisible(false);
+
+
+ region.style().dataStyle().lineStyle().setParameter("colorRotateMethod", "regionOverlayIndex");
+ region.style().dataStyle().lineStyle().setParameter("colorRotate", "black, red, green, blue");
+ region.style().dataStyle().lineStyle().setParameter("thickness", "3");
+ region.style().dataStyle().errorBarStyle().setVisible(false);
+ System.out.println(Arrays.toString(region.style().dataStyle().lineStyle().availableParameterOptions("colorRotateMethod")));
+ }
+
+ public static void main(String arg[]){
+ IAnalysisFactory af = IAnalysisFactory.create();
+ IHistogramFactory hf = af.createHistogramFactory(af.createTreeFactory().create());
+
+ IPlotter p = af.createPlotterFactory().create();
+ p.createRegions(1, 2);
+ IHistogram1D h1 = hf.createHistogram1D("blah", 100, -5, 5);
+ IHistogram1D h2 = hf.createHistogram1D("bleh", 100, -5, 5);
+ Random random = new Random();
+ for(int i = 0; i< 100000; i++){
+ h1.fill(random.nextGaussian());
+ h2.fill(random.nextGaussian()*2);
+ }
+ hideLegendAndStats(p.region(1));
+ noFillHistogramBars(p.region(0));
+ stylize(p.region(0), "title", "x axis label", "y axis label");
+ stylize(p.region(1), "stuff", "x axis label", "y axis label");
+ p.region(0).plot(h1);
+ p.region(0).plot(h2);
+
+
+ IHistogram2D h3 = hf.createHistogram2D("blah", 100, -5, 5, 100, -5,5);
+
+ for(int i = 0; i< 100000; i++){
+ h3.fill(random.nextGaussian(), random.nextGaussian());
+ }
+
+ p.region(1).plot(h3);
+
+
+
+ p.show();
+
+ }
+ public static void hideLegendAndStats(IPlotterRegion r){
+ r.style().statisticsBoxStyle().setVisible(false);
+ r.style().legendBoxStyle().setVisible(false);
+ }
}
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/SumEverything.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/SumEverything.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/SumEverything.java Tue Jan 26 18:37:32 2016
@@ -2,6 +2,7 @@
import java.io.File;
import java.io.IOException;
+import java.util.ArrayList;
import hep.aida.IAnalysisFactory;
import hep.aida.IHistogram1D;
@@ -23,20 +24,44 @@
oneArg(arg[0]);
}
}
- static void oneArg(String indir) throws IllegalArgumentException, IOException{
+ static void oneArg(final String indir) throws IllegalArgumentException, IOException{
File outdir = new File(indir + "/sums");
outdir.mkdir();
- for(File subdirf : new File(indir).listFiles()){
- String subdir = subdirf.getAbsolutePath();
- if(subdir.matches(".*/sums/?"))
- continue;
- String split[] = subdir.split("/");
- String outfile = indir + "/sums/" + split[split.length-1] + ".aida";
- twoArg(subdir, outfile);
+ ArrayList<Thread> threads = new ArrayList<Thread>();
+ for(final File subdirf : new File(indir).listFiles()){
+ Thread t = new Thread(){
+ public void run(){
+
+ String subdir = subdirf.getAbsolutePath();
+ if(subdir.matches(".*/sums/?"))
+ return;
+ String split[] = subdir.split("/");
+ String outfile = indir + "/sums/" + split[split.length-1] + ".aida";
+ try {
+ twoArg(subdir, outfile);
+ } catch (IllegalArgumentException e) {
+ // TODO Auto-generated catch block
+ e.printStackTrace();
+ } catch (IOException e) {
+ // TODO Auto-generated catch block
+ e.printStackTrace();
+ }
+ }
+ };
+ threads.add(t);
+ t.start();
+ }
+ for(Thread t : threads){
+ try {
+ t.join();
+ } catch (InterruptedException e) {
+ // TODO Auto-generated catch block
+ e.printStackTrace();
+ }
}
new File(indir + "/sums/total.aida").delete();
twoArg(outdir.getAbsolutePath(), indir + "/sums/total.aida");
-
+
}
static void twoArg(String indir, String out) throws IllegalArgumentException, IOException{
@@ -45,45 +70,56 @@
new File(out).delete();
ITree tree0 = tf.create(out, "xml", false, true);
IHistogramFactory hf = af.createHistogramFactory(tree0);
-
+
int j = 0;
for(File s : new File(indir).listFiles()){
- ITree tree = af.createTreeFactory().create(s.getAbsolutePath(),"xml");
+ if(!s.getAbsolutePath().endsWith("aida"))
+ continue;
+ try{
+ ITree tree = af.createTreeFactory().create(s.getAbsolutePath(),"xml");
-
- if(j == 0){
- String [] names = tree.listObjectNames();
- tree0.mount("/tmp", tree, "/");
- for(String name : names){
- Object o = tree.find(name);
- if(o instanceof IHistogram1D || o instanceof IHistogram2D)
- tree0.cp("/tmp/" + name, name);
+
+ if(j == 0){
+ String [] names = tree.listObjectNames();
+ tree0.mount("/tmp", tree, "/");
+ for(String name : names){
+ Object o = tree.find(name);
+ if(o instanceof IHistogram1D || o instanceof IHistogram2D)
+ tree0.cp("/tmp/" + name, name);
+ }
+ tree0.unmount("/tmp");
+ tree.close();
+
}
- tree0.unmount("/tmp");
+ else{
+ //tree.
+ String [] names = tree.listObjectNames();
+ tree0.mount("/tmp", tree, "/");
+ for(String name : names){
+ Object o = tree.find(name);
+ if(o instanceof IHistogram1D)
+ ((IHistogram1D)tree0.find(name)).add((IHistogram1D)o);
+ if(o instanceof IHistogram2D)
+ ((IHistogram2D)tree0.find(name)).add((IHistogram2D)o);
+ }
+ tree0.unmount("/tmp");
+ tree.close();
+ }
+
tree.close();
-
+ j++;
+ System.out.println(j + " files have been read");
+
+
+ } catch(IllegalArgumentException e){
+ //print the filename
+ System.out.println("Exception happened at file " + s.getAbsolutePath());
+
+ e.printStackTrace();
}
- else{
- //tree.
- String [] names = tree.listObjectNames();
- tree0.mount("/tmp", tree, "/");
- for(String name : names){
- Object o = tree.find(name);
- if(o instanceof IHistogram1D)
- ((IHistogram1D)tree0.find(name)).add((IHistogram1D)o);
- if(o instanceof IHistogram2D)
- ((IHistogram2D)tree0.find(name)).add((IHistogram2D)o);
- }
- tree0.unmount("/tmp");
- tree.close();
- }
-
- tree.close();
- j++;
- System.out.println(j + " files have been read");
+ tree0.commit();
}
- tree0.commit();
}
}
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java Tue Jan 26 18:37:32 2016
@@ -6,7 +6,7 @@
public class BinGenerator {
public static void main(String arg[]) throws FileNotFoundException{
- int nBins = 20;
+ int nBins = 32;
PrintStream pw = new PrintStream("generatedbins.txt");
pw.println(nBins);
@@ -33,17 +33,21 @@
bins[i] = thetaMin+dTheta*i;
}
return bins; */
- double thetaMax = .145;
- double thetaMin = .035;
+ double thetaMax = .200;
+ double thetaMin = .040;
double[] bins = new double[nBins +1];
- double xMin = 1/(thetaMax*thetaMax);
+ for(int i = 0; i<nBins+1; i++){
+ bins[i] = thetaMin+i*(thetaMax-thetaMin)/nBins;
+ }
+ return bins;
+ /*double xMin = 1/(thetaMax*thetaMax);
double xMax = 1/(thetaMin*thetaMin);
for(int i = 0; i< nBins+1; i++){
double x = xMax - i*(xMax-xMin)/nBins;
bins[i] = Math.pow(x, -.5);
}
- return bins;
+ return bins;*/
}
private static double[] getPhiBounds(double thetaMin, double thetaMax){
double phiBins[] = new double[6];
@@ -52,8 +56,14 @@
boolean prevInRange = false;
for(double phi = 0; phi< 3.14; phi+= dphi){
- boolean inRange = EcalUtil.fid_ECal_spherical(thetaMin, phi) && EcalUtil.fid_ECal_spherical(thetaMax, phi)
- && EcalUtil.fid_ECal_spherical(thetaMin, -phi) && EcalUtil.fid_ECal_spherical(thetaMax, -phi);
+
+ // make the angular cuts on the tracks such that the particles that go into that cut
+ // are expected to be within 4 mm (~= 2 times the angular resolution of 1.5 mrad) of
+ // the ecal cuts.
+ double d = 4;
+
+ boolean inRange = EcalUtil.fid_ECal_spherical_more_strict(thetaMin, phi, d) && EcalUtil.fid_ECal_spherical_more_strict(thetaMax, phi, d)
+ && EcalUtil.fid_ECal_spherical_more_strict(thetaMin, -phi, d) && EcalUtil.fid_ECal_spherical_more_strict(thetaMax, -phi, d);
if(inRange && !prevInRange)
phiBins[edgeNumber++] = phi;
if(prevInRange && !inRange)
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/CustomBinning.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/CustomBinning.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/CustomBinning.java Tue Jan 26 18:37:32 2016
@@ -76,21 +76,5 @@
}
return tot;
}
- /**
- * @param bin
- * @param a = 2E/M
- * @return the integral of 1/sin^4(th/2)*cos^2(th/2)/(1+a*sin^2(th/2)) times dPhi,
- * which appears in the integral of mott scattering.
- */
- public double mottIntegralFactor(double a, int bin){
- double dPhi = 0;
- for(int i = 0; i< phiMax[bin].length; i++)
- dPhi += 2*(phiMax[bin][0] - phiMin[bin][0]); //factor of 2 from top and bottom
-
- double csc2 = Math.pow(Math.sin(thetaMax[bin]/2), -2);
- double Imax = (-csc2+(1+a)*Math.log(a+2*csc2));
- csc2 = Math.pow(Math.sin(thetaMin[bin]/2), -2);
- double Imin = (-csc2+(1+a)*Math.log(a+2*csc2));
- return 2*dPhi*(Imax-Imin);
- }
+
}
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/DisplayHistograms.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/DisplayHistograms.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/DisplayHistograms.java Tue Jan 26 18:37:32 2016
@@ -4,6 +4,7 @@
public class DisplayHistograms {
public static void main(String arg[]) throws IllegalArgumentException, IOException{
+ //System.out.println("dognabo");
MakeHistograms.main(new String[]{arg[0]});
}
}
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/EcalUtil.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/EcalUtil.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/EcalUtil.java Tue Jan 26 18:37:32 2016
@@ -12,6 +12,45 @@
import org.lcsim.event.Cluster;
public class EcalUtil {
+
+ static int[] getCrystalIndex(Cluster c){
+ if(map == null){
+ try{
+ readMap();
+ }catch(Exception e){
+ e.printStackTrace();
+ }
+ }
+ c.getPosition();
+ int bestix = 0, bestiy = 0;
+ double bestdist = 100000;
+
+ double cx = c.getPosition()[0];
+ double cy = c.getPosition()[1];
+ for(int ix = -23; ix<= 23; ix++){
+ if(!map.containsKey(200+ix))
+ continue;
+ double x = map.get(200+ix)[0];
+ if(Math.abs(x-cx)<bestdist){
+ bestdist = Math.abs(x-cx);
+ bestix = ix;
+ }
+ }
+
+ bestdist = 100000;
+ for(int iy = -5; iy<= 5; iy++){
+ if(!map.containsKey(100*iy + bestix))
+ continue;
+ double y = map.get(100*iy + bestix)[1];
+ if(Math.abs(y-cy)<bestdist){
+ bestdist = Math.abs(y-cy);
+ bestiy = iy;
+ }
+ }
+ return new int[]{bestix, bestiy};
+ }
+
+
/**
* Holly's algorithm.
* @param x x position of cluster
@@ -83,6 +122,40 @@
return in_fid;
}
+ /**
+ *
+ * @param x
+ * @param y
+ * @param d the additional distance from the edge of the ecal in addition
+ * to what is required by fid_Cal(double, double)
+ * @return
+ */
+ public static boolean fid_ecal_more_strict(double x, double y, double d){
+ y = Math.abs(y);
+
+ boolean in_fid = false;
+
+ double x_edge_low = -262.74 + d;
+ double x_edge_high = 347.7 - d;
+ double y_edge_low = 33.54 + d;
+ double y_edge_high = 75.18 - d;
+
+ double x_gap_low = -106.66 - d;
+ double x_gap_high = 42.17 + d;
+ double y_gap_high = 47.18 + d;
+
+ y = Math.abs(y);
+
+ if( x > x_edge_low && x < x_edge_high && y > y_edge_low && y < y_edge_high )
+ {
+ if( !(x > x_gap_low && x < x_gap_high && y > y_edge_low && y < y_gap_high) )
+ {
+ in_fid = true;
+ }
+ }
+
+ return in_fid;
+ }
static double[] toSphericalFromBeam(double pxpz, double pypz){
double x = pxpz, y = pypz, z = 1;
double beamTilt = .03057;
@@ -96,10 +169,10 @@
return new double[]{theta, phi};
}
- static Map<Integer, double[]> map = new HashMap();
+ static Map<Integer, double[]> map ;
static void readMap() throws FileNotFoundException{
- Scanner s = new Scanner(new File("ecal_positions.txt"));
-
+ Scanner s = new Scanner(new File(System.getenv("HOME") + "/ecal_positions.txt"));
+ map = new HashMap();
while(s.hasNext()){
int ix =s.nextInt();
@@ -154,6 +227,12 @@
double y = xy[1];
return fid_ECal(x, y);
}
+ public static boolean fid_ECal_spherical_more_strict(double theta, double phi, double d){
+ double[] xy = getXY(theta, phi);
+ double x = xy[0];
+ double y = xy[1];
+ return fid_ecal_more_strict(x, y, d);
+ }
public static void main(String arg[]){
double x = 0, y = 0;
double sp[] = getThetaPhiSpherical(x,y);
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/MakeHistograms.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/MakeHistograms.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/MakeHistograms.java Tue Jan 26 18:37:32 2016
@@ -2,12 +2,15 @@
import java.io.File;
import java.io.IOException;
+import java.util.ArrayList;
+import java.util.Arrays;
import java.util.Collections;
import java.util.Comparator;
import java.util.Date;
import java.util.List;
import org.hps.conditions.ConditionsDriver;
+import org.hps.recon.tracking.TrackType;
import org.hps.record.triggerbank.AbstractIntData;
import org.hps.record.triggerbank.TIData;
import org.lcsim.event.CalorimeterHit;
@@ -15,9 +18,12 @@
import org.lcsim.event.EventHeader;
import org.lcsim.event.GenericObject;
import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackState;
import org.lcsim.lcio.LCIOReader;
import hep.aida.IAnalysisFactory;
+import hep.aida.IDataPointSetFactory;
import hep.aida.IHistogram1D;
import hep.aida.IHistogram2D;
import hep.aida.IHistogramFactory;
@@ -29,6 +35,9 @@
import org.hps.users.spaul.StyleUtil;
public class MakeHistograms {
+
+
+
static boolean display = false;
static CustomBinning cb;
public static void main(String arg[]) throws IllegalArgumentException, IOException{
@@ -40,11 +49,11 @@
path = "temp.aida";
}
IAnalysisFactory af = IAnalysisFactory.create();
- ITreeFactory tf = af.createTreeFactory();
+ ITreeFactory tf = af.createTreeFactory();
ITree tree0 = tf.create(path, "xml");
extractHistograms(tree0);
setupPlotter(af);
-
+
} else{
String input = arg[0];
@@ -86,104 +95,206 @@
}
}
-
+
static IHistogram2D h1, h2, h2a, h2b, h2c;
static IHistogram2D h4,h4a;
- static IHistogram1D h3, h3a;
- static IHistogram1D h5, h6;
-
+ static IHistogram1D h3, h3a, h3_t, h3_b;
+ static IHistogram1D h5, h5a;
+ static IHistogram2D h6, h6a;
+ static IHistogram1D h7, h7a;
+ static IHistogram1D h8;
+ static IHistogram1D h9_t, h9_b;
+ static IHistogram1D h10_t, h10_b;
+
private static void extractHistograms(ITree tree0) {
+ h1 = (IHistogram2D) tree0.find("theta vs energy");
+
h2 = (IHistogram2D) tree0.find("theta vs phi");
h2a = (IHistogram2D) tree0.find("theta vs phi cut");
h2b = (IHistogram2D) tree0.find("theta vs phi cut alt");
h2c = (IHistogram2D) tree0.find("theta vs phi alt");
-
+
h3 = (IHistogram1D) tree0.find("theta");
+ h3a = (IHistogram1D) tree0.find("theta isolated ");
+ h3_t = (IHistogram1D) tree0.find("theta top");
+ h3_b = (IHistogram1D) tree0.find("theta bottom");
+
h4 = (IHistogram2D) tree0.find("px\\/pz vs py\\/pz");
h4a = (IHistogram2D) tree0.find("px\\/pz vs py\\/pz cut");
- h5 = (IHistogram1D) tree0.find("energy");
-
+ System.out.println(h4a.xAxis().bins());
+ h5 = (IHistogram1D) tree0.find("energy top");
+ h5 = (IHistogram1D) tree0.find("energy bottom");
+
+ h6 = (IHistogram2D) tree0.find("cluster");
+ h6a = (IHistogram2D) tree0.find("cluster matched");
+ h7 = (IHistogram1D) tree0.find("y top");
+ h7a = (IHistogram1D) tree0.find("y bottom");
+ h8 = (IHistogram1D) tree0.find("seed energy");
+
+
+ h9_t = (IHistogram1D) tree0.find("pz top");
+ h9_b = (IHistogram1D) tree0.find("pz bottom");
+
+
+ h10_t = (IHistogram1D) tree0.find("clustsize top");
+ h10_b = (IHistogram1D) tree0.find("clustsize bottom");
+
}
static void setupHistograms(IHistogramFactory hf){
//h1 = hf.createHistogram2D("px\\/pz vs py\\/pz", 160, -.16, .24, 160, -.2, .2);
-
-
-
-
+
+
+
+
h2 = hf.createHistogram2D("theta vs phi", 300, 0, .3, 314, -3.14, 3.14);
-
+
h2a = hf.createHistogram2D("theta vs phi cut", 300, 0, .3, 314, -3.14, 3.14);
double thetaBins[] = new double[cb.nTheta+1];
for(int i = 0; i<cb.nTheta; i++){
thetaBins[i] = cb.thetaMin[i];
}
-
+
thetaBins[thetaBins.length-1] = cb.thetaMax[cb.nTheta-1];
double phiBins[] = new double[315];
for(int i = 0; i<315; i++){
phiBins[i] = i/50.-3.14; //every 10 mrad;
}
-
+
+ double eBins[] = new double[66];
+ for(int i = 0; i<66; i++){
+ eBins[i] = i/50.; //every 20 MeV up to 1300 MeV
+ }
+
+
+ h1 = hf.createHistogram2D("theta vs energy", "theta vs energy", thetaBins, eBins);
+
+
//identical to h2a, except different binning
h2b = hf.createHistogram2D("theta vs phi cut alt", "theta vs phi cut alt", thetaBins, phiBins);
h2c = hf.createHistogram2D("theta vs phi alt", "theta vs phi alt", thetaBins, phiBins);
-
+
h3 = hf.createHistogram1D("theta", "theta", thetaBins);
-
- h4 = hf.createHistogram2D("px\\/pz vs py\\/pz", 160, -.16, .24, 160, -.2, .2);
- h4a = hf.createHistogram2D("px\\/pz vs py\\/pz cut", 160, -.16, .24, 160, -.2, .2);
-
- h5 = hf.createHistogram1D("energy", 75, 0, 1.5);
+ h3a = hf.createHistogram1D("theta isolated ", "theta isolated", thetaBins);
+
+ h3_t = hf.createHistogram1D("theta top", "theta top", thetaBins);
+ h3_b = hf.createHistogram1D("theta bottom", "theta bottom", thetaBins);
+
+
+ h4 = hf.createHistogram2D("px\\/pz vs py\\/pz", 300, -.16, .24, 300, -.2, .2);
+ h4a = hf.createHistogram2D("px\\/pz vs py\\/pz cut", 300, -.16, .24, 300, -.2, .2);
+
+ h5 = hf.createHistogram1D("energy top", 75, 0, 1.5);
+ h5a = hf.createHistogram1D("energy bottom", 75, 0, 1.5);
+
+
+ h9_t = hf.createHistogram1D("pz top", 75, 0, 1.5);
+ h9_b = hf.createHistogram1D("pz bottom", 75, 0, 1.5);
+
+ h6 = hf.createHistogram2D("cluster", 47, -23.5, 23.5, 11, -5.5, 5.5);
+ h6a = hf.createHistogram2D("cluster matched", 47, -23.5, 23.5, 11, -5.5, 5.5);
+
+ h7 = hf.createHistogram1D("y top", 500, 0, 100);
+
+ h7a = hf.createHistogram1D("y bottom", 500, 0, 100);
+
+ h8 = hf.createHistogram1D("seed energy", 120, 0, 1.2);
+
+ h10_t = hf.createHistogram1D("clustsize top", 10,0, 10);
+ h10_b = hf.createHistogram1D("clustsize bottom", 10,0, 10);
}
static void setupPlotter(IAnalysisFactory af){
IPlotterFactory pf = af.createPlotterFactory();
IPlotter p = pf.create();
p.createRegions(2,2);
p.region(0).plot(h2);
-
StyleUtil.stylize(p.region(0), "theta", "phi");
- p.region(1).plot(h2a);
- StyleUtil.stylize(p.region(1), "theta", "phi");
+ p.region(1).plot(h3a);
+ StyleUtil.stylize(p.region(1), "theta", "# of particles");
p.region(2).plot(h3);
StyleUtil.stylize(p.region(2), "theta", "# of particles");
p.region(3).plot(h5);
+ p.region(3).plot(h5a);
StyleUtil.stylize(p.region(3), "energy", "# of particles");
-
+
p.show();
+
//new window for the next plot
IPlotter p2 = pf.create();
p2.region(0).plot(h2b);
+ //IDataPointSetFactory dpsf = af.createDataPointSetFactory(af.createTreeFactory().create());
+
StyleUtil.stylize(p2.region(0), "theta", "phi");
-
+
p2.show();
-
+
//new window for the next plot
IPlotter p3 = pf.create();
p3.region(0).plot(h2c);
StyleUtil.stylize(p3.region(0), "theta", "phi");
-
+
p3.show();
-
+
//new window for the next plot
IPlotter p4 = pf.create();
p4.region(0).plot(h4);
StyleUtil.stylize(p4.region(0), "px/pz", "py/pz");
-
+
p4.show();
-
+
//new window for the next plot
IPlotter p5 = pf.create();
p5.region(0).plot(h4a);
StyleUtil.stylize(p5.region(0), "px/pz", "py/pz");
-
+
p5.show();
+
+ IPlotter p6 = pf.create("efficiency");
+ p6.createRegions(1,2);
+ p6.region(0).plot(h6);
+ StyleUtil.stylize(p6.region(0), "ix", "iy");
+ p6.region(1).plot(h6a);
+ StyleUtil.stylize(p6.region(1), "ix", "iy");
+ p6.show();
+
+ IPlotter p7 = pf.create("theta vs energy");
+ //p6.createRegions(1,2);
+ p7.region(0).plot(h1);
+ StyleUtil.stylize(p7.region(0), "theta", "energy");
+ // StyleUtil.stylize(p6.region(1), "ix", "iy");
+ p7.show();
+
+ IPlotter p8 = pf.create("y");
+ //p6.createRegions(1,2);
+ p8.region(0).plot(h7);
+ p8.region(0).plot(h7a);
+ StyleUtil.stylize(p8.region(0), "y", "# of particles");
+ // StyleUtil.stylize(p6.region(1), "ix", "iy");
+ p8.show();
+
+ IPlotter p9 = pf.create("theta: top vs. bottom");
+ //p6.createRegions(1,2);
+ p9.region(0).plot(h3_t);
+ p9.region(0).plot(h3_b);
+ StyleUtil.stylize(p9.region(0), "theta", "theta", "# of particles");
+ StyleUtil.noFillHistogramBars(p9.region(0));
+ StyleUtil.stylize(p6.region(1), "ix", "iy");
+ p9.show();
+
+ IPlotter p10 = pf.create("seed energy");
+ //p6.createRegions(1,2);
+ p10.region(0).plot(h8);
+ StyleUtil.stylize(p10.region(0), "seed energy", "seed energy (GeV)", "# of particles");
+ //StyleUtil.noFillHistogramBars(p10.region(0));
+ StyleUtil.stylize(p6.region(1), "ix", "iy");
+ p10.show();
+
}
private static void processEvent(EventHeader event) {
- if(event.getEventNumber() %1000 == 0)
+ if(event.getEventNumber() %10000 == 0)
System.out.println("event number " + event.getEventNumber());
-
+
for (GenericObject gob : event.get(GenericObject.class,"TriggerBank"))
{
if (!(AbstractIntData.getTag(gob) == TIData.BANK_TAG)) continue;
@@ -194,11 +305,61 @@
}
}
List<ReconstructedParticle> particles = event.get(ReconstructedParticle.class, "FinalStateParticles");
-
- for(ReconstructedParticle p : particles){
+ particles = RemoveDuplicateParticles.removeDuplicateParticles(particles);
+ outer : for(ReconstructedParticle p : particles){
+ //check if this is a duplicate particle (ie, different track same cluster)
+ /*for(ParticleInfo cand : candidates){
+ if(p.getEnergy() == cand.E && cand.isGood == true)
+ continue outer;
+ }*/
boolean isGood = addParticle(p);
-
+ if(!isGood){
+ if(p.getClusters().size()!= 0)
+ candidates.add(new ParticleInfo(0, p.getEnergy(), p.getClusters().get(0).getCalorimeterHits().get(0).getTime(), false));
+
+ }
+
+ }
+ for(ParticleInfo c : candidates){
+ if(c.isGood){
+ boolean isIsolated = true;
+ for(ParticleInfo c2 :candidates){
+ //try to remove events that have possible mollers in them.
+ //correct for this later.
+ if(Math.abs(c2.t - c.t) < 2.5 && c2.E < .3 && c2 != c){
+ isIsolated = false;
+ break;
+ }
+ }
+ if(isIsolated){
+ h3a.fill(c.theta);
+ }
+ }
+ }
+ candidates.clear();
+ processEventEfficiency(particles);
+ }
+
+ private static void processEventEfficiency(List<ReconstructedParticle> parts) {
+ for(ReconstructedParticle p : parts){
+ if(p.getClusters().size() == 0)
+ continue;
+ if(!(p.getEnergy() > eMin && p.getEnergy() < eMax))
+ continue;
+ Cluster c = p.getClusters().get(0);
+ if(!EcalUtil.fid_ECal(c) || c.getCalorimeterHits().size() < 3)
+ continue;
+ //good cluster? now continue
+ int ixiy[] = EcalUtil.getCrystalIndex(c);
+
+ h6.fill(ixiy[0], ixiy[1]);
+ if(p.getTracks().size() == 0)
+ continue;
+ Track t = p.getTracks().get(0);
+ if(t.getChi2()> maxChi2)
+ continue;
+ h6a.fill(ixiy[0], ixiy[1]);
}
}
@@ -209,44 +370,129 @@
static double beamTilt = .03057;
static double maxChi2 = 50;
+ //maximum difference between the reconstructed energy and momentum
+ static double maxdE = .3;
+
+ static double seedEnergyCut = .4;
+
+
+ static ArrayList<ParticleInfo> candidates = new ArrayList();
static boolean addParticle(ReconstructedParticle part){
-
- if(part.getTracks().size() == 0)
+
+
+
+ if(part.getCharge() != -1)
return false;
- if(part.getTracks().get(0).getChi2()>maxChi2){
- return false;
- }
if(part.getClusters().size() == 0)
return false;
Cluster c = part.getClusters().get(0);
double time = c.getCalorimeterHits().get(0).getTime();
+
+ if(!(time>40 && time <50))
+ return false;
+ double seedEnergy = 0;
+ for(CalorimeterHit hit : c.getCalorimeterHits()){
+ if(hit.getCorrectedEnergy() > seedEnergy)
+ seedEnergy = hit.getCorrectedEnergy();
+ }
+ h8.fill(seedEnergy);
+
+
+ if(seedEnergy < seedEnergyCut)
+ return false;
+
+ if(c.getPosition()[1] > 0){
+ h10_t.fill(c.getSize());
+ }
+ else{
+ h10_b.fill(c.getSize());
+ }
+
+
+ if(c.getCalorimeterHits().size() < 3)
+ return false;
+
+
+ if(c.getEnergy() > eMin && c.getEnergy() < eMax){
+ if(c.getPosition()[1] > 0)
+ h7.fill(c.getPosition()[1]);
+ else if(c.getPosition()[1] < 0)
+ h7a.fill(-c.getPosition()[1]);
+ }
+
if(EcalUtil.fid_ECal(c)){
- if(c.getCalorimeterHits().size() < 3)
+
+ if(part.getTracks().size() == 0)
return false;
- if(time>40 && time <48)
+ Track t = part.getTracks().get(0);
+ if(t.getChi2()>maxChi2){
+ return false;
+ }
+ if(!TrackType.isGBL(t.getType()))
+ return false;
+ if(c.getPosition()[1] > 0){
h5.fill(c.getEnergy());
- if(c.getEnergy() > eMin && c.getEnergy() < eMax && (time >40 && time < 48)) {
-
- Hep3Vector p = part.getMomentum();
-
- double px = p.x(), pz = p.z();
- double pxtilt = px*Math.cos(beamTilt)-pz*Math.sin(beamTilt);
- double py = p.y();
- double pztilt = pz*Math.cos(beamTilt)+px*Math.sin(beamTilt);
-
- double theta = Math.atan(Math.hypot(pxtilt, py)/pztilt);
- double phi =Math.atan2(py, pxtilt);
+ }
+ else{
+ h5a.fill(c.getEnergy());
+ }
+
+
+ Hep3Vector p = part.getMomentum();
+
+
+
+ double px = p.x(), pz = p.z();
+ double pxtilt = px*Math.cos(beamTilt)-pz*Math.sin(beamTilt);
+ double py = p.y();
+ double pztilt = pz*Math.cos(beamTilt)+px*Math.sin(beamTilt);
+
+ if(Math.abs(pztilt - c.getEnergy()) > maxdE)
+ return false;
+ if(c.getPosition()[1] > 0)
+ h9_t.fill(pztilt);
+ else
+ h9_b.fill(pztilt);
+
+ double theta = Math.atan(Math.hypot(pxtilt, py)/pztilt);
+ double phi =Math.atan2(py, pxtilt);
+ boolean inRange = cb.inRange(theta, phi);
+ if(inRange)
+ h1.fill(theta, c.getEnergy());
+
+
+
+ if(c.getEnergy() > eMin && c.getEnergy() < eMax) {
+
+
h2.fill(theta, phi);
h2c.fill(theta, phi);
-
+
h4.fill(px/pz, py/pz);
-
- if(cb.inRange(theta, phi)){
+
+ if(inRange){
+
+ //System.out.println(c.getEnergy() + " " + t.getType());
+ /*for(TrackState ts : t.getTrackStates()){
+ if(ts.getLocation() == TrackState.AtIP)
+ System.out.println(Arrays.toString(
+ ts.getReferencePoint()));
+ }*/
h2a.fill(theta, phi);
h2b.fill(theta, phi);
+
h3.fill(theta);
+ if(py > 0)
+ h3_t.fill(theta);
+ else
+ h3_b.fill(theta);
+ //if(h3_t.sumBinHeights()+h3_b.sumBinHeights() != h3.sumBinHeights())
+ //System.out.println("NABO ERROR");
+
+
h4a.fill(px/pz, py/pz);
+ candidates.add(new ParticleInfo(theta, c.getEnergy(), c.getCalorimeterHits().get(0).getTime(), true));
}
@@ -256,5 +502,16 @@
}
return false;
}
-
+ static class ParticleInfo{
+ double theta;
+ double E;
+ double t;
+ boolean isGood;
+ ParticleInfo(double theta, double E, double t, boolean isGood){
+ this.theta = theta;
+ this.E = E;
+ this.t = t;
+ this.isGood = isGood;
+ }
+ }
}
Modified: java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java
=============================================================================
--- java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java (original)
+++ java/branches/jeremy-dev/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java Tue Jan 26 18:37:32 2016
@@ -2,10 +2,16 @@
import java.awt.Canvas;
import java.awt.Color;
+import java.awt.Container;
+import java.awt.Font;
import java.awt.Graphics;
+import java.awt.Graphics2D;
+import java.awt.image.BufferedImage;
import java.io.File;
import java.io.FileNotFoundException;
-
+import java.io.IOException;
+
+import javax.imageio.ImageIO;
import javax.swing.JFrame;
import hep.aida.IAnalysisFactory;
@@ -27,15 +33,47 @@
public static void main(String arg[]) throws FileNotFoundException{
JFrame frame = new JFrame();
frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
- frame.add(new ShowCustomBinning(new File(arg[0])));
- frame.setSize(800, 800);
+ Canvas c = new ShowCustomBinning(new File(arg[0]));
+ String outdir = arg[1];
+ frame.add(c);
+ frame.setSize(1200, 800);
frame.setVisible(true);
+
+ try {
+ BufferedImage im = new BufferedImage(c.getWidth(), c.getHeight(), BufferedImage.TYPE_INT_ARGB);
+ c.paint(im.getGraphics());
+ ImageIO.write(im, "PNG", new File(outdir +"/bins.png"));
+ } catch (IOException e) {
+ // TODO Auto-generated catch block
+ e.printStackTrace();
+ }
+
+ frame = new JFrame();
+ frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
+ c = new ShowCustomBinningXY(new File(arg[0]));
+ frame.add(c);
+ frame.setSize(1200, 615);
+ frame.setVisible(true);
+
+ try {
+ BufferedImage im = new BufferedImage(c.getWidth(), c.getHeight(), BufferedImage.TYPE_INT_ARGB);
+ c.paint(im.getGraphics());
+ ImageIO.write(im, "PNG", new File(outdir + "/bins_xy.png"));
+ } catch (IOException e) {
+ // TODO Auto-generated catch block
+ e.printStackTrace();
+ }
}
public void paint(Graphics g){
+ g.setFont(new Font(Font.DIALOG, Font.PLAIN, 24));
+
drawEcalOutline(g);
drawFidEcalOutline(g);
drawCustomBinRectangles(g);
+ g.setColor(Color.BLACK);
+ drawXAxis(g);
+ drawYAxis(g);
}
void drawFidEcalOutline(Graphics g){
@@ -95,15 +133,24 @@
}
void drawEcalOutline(Graphics g){
+ /*double x_edge_low = -262.74;
+ double x_edge_high = 347.7;
+ double y_edge_low = 33.54;
+ double y_edge_high = 75.18;
+
+ double x_gap_low = -106.66;
+ double x_gap_high = 42.17;
+ double y_gap_high = 47.18;*/
+
g.setColor(Color.BLACK);
- double x_edge_low = -276.50;
- double x_edge_high = 361.55;
- double y_edge_low = 20.17;
- double y_edge_high = 89;
-
- double x_gap_low = -93.30;
- double x_gap_high = 28.93;
- double y_gap_high = 33.12;
+ double x_edge_low = -269.56;
+ double x_edge_high = 354.52;
+ double y_edge_low = 26.72;
+ double y_edge_high = 82;
+
+ double x_gap_low = -99.84;
+ double x_gap_high = 33.35;
+ double y_gap_high = 40.36;
double x1,y1, x2, y2;
double nPoints = 500;
for(int i = 0; i< nPoints-1; i++){
@@ -157,6 +204,8 @@
}
Color altBin1 = new Color(0, 0, 128);
Color altBin2 = new Color(0,128,0);
+ Color fillBin1 = new Color(196, 196, 255);
+ Color fillBin2 = new Color(196,255,196);
void drawCustomBinRectangles(Graphics g){
for(int i = 0; i<binning.nTheta; i++){
g.setColor(i%2 == 0 ? altBin1 : altBin2);
@@ -167,13 +216,45 @@
double theta2 = binning.thetaMax[i];
int x =getX(theta1)+1, y = getY(phi1), w = getX(theta2)-getX(theta1), h = getY(phi2)-getY(phi1);
+
+ g.setColor(i%2 == 0 ? fillBin1 : fillBin2);
+ g.fillRect(x, y, w, h);
+ g.setColor(i%2 == 0 ? altBin1 : altBin2);
g.drawRect(x, y, w, h);
x =getX(theta1)+1; y = getY(-phi2); w = getX(theta2)-getX(theta1); h = getY(-phi1)-getY(-phi2);
-
- g.drawRect(x, y, w, h);
-
-
+ g.setColor(i%2 == 0 ? fillBin1 : fillBin2);
+ g.fillRect(x, y, w, h);
+ g.setColor(i%2 == 0 ? altBin1 : altBin2);
+ g.drawRect(x, y, w, h);
}
+
+ }
+ }
+ void drawXAxis(Graphics g){
+ //x axis
+ g.drawString("θ", getX(.28), getY(0) - 40);
+ g.drawLine(getX(0), getY(0), getX(.28), getY(0));
+ for(int i = 0; i< 28; i++){
+ if(i%5 == 0 && i != 0){
+ g.drawString(i/100.+"", getX(i/100.)-15, getY(0)-20);
+ g.drawLine(getX(i/100.), getY(0), getX(i/100.), getY(0)-15);
+ }
+ g.drawLine(getX(i/100.), getY(0), getX(i/100.), getY(0)-5);
+ }
+ }
+ void drawYAxis(Graphics g){
+ g.drawString("Ï", getX(0)+70, getY(3));
+ g.drawLine(getX(0), getY(-3), getX(0), getY(3));
+ for(int i = -30; i<= 30; i++){
+ if(i%5 == 0 && i != 0){
+ g.drawString(i/10.+"", getX(0)+20, getY(i/10.) + 5);
+
+ g.drawLine(getX(0), getY(i/10.), getX(0) + 15, getY(i/10.));
+ }
+ if(i == 0){
+ //g.drawString(i/10.+"", getX(0)+10, getY(i/10.) - 15);
+ }
+ g.drawLine(getX(0), getY(i/10.), getX(0) + 5, getY(i/10.));
}
}
@@ -186,8 +267,9 @@
}
int getX(double theta){
- return (int)(this.getWidth()*theta/.3);
- }
+ return (int)(this.getWidth()*theta/.3)+left_margin;
+ }
+ int left_margin = 20;
int getY(double phi){
return (int)(this.getHeight()*(3.2-phi)/6.4);
}
|
