diff --git a/detector-model/src/main/java/org/lcsim/geometry/compact/converter/HPSTracker2019GeometryDefinition.java b/detector-model/src/main/java/org/lcsim/geometry/compact/converter/HPSTracker2019GeometryDefinition.java index a9afe72e37..07beb362a1 100644 --- a/detector-model/src/main/java/org/lcsim/geometry/compact/converter/HPSTracker2019GeometryDefinition.java +++ b/detector-model/src/main/java/org/lcsim/geometry/compact/converter/HPSTracker2019GeometryDefinition.java @@ -87,7 +87,7 @@ public void build() { //LOGGER.info("Construct uChannelL14Bottom"); - UChannelL13 uChannelL14Bottom = new UChannelL14Bottom("support_bottom_L14", svtBox, alignmentCorrections, + UChannelL13 uChannelL14Bottom = new UChannelL14Bottom("support_bottom_L13", svtBox, alignmentCorrections, supportRingKinL13Bottom); surveyVolumes.add(uChannelL14Bottom); @@ -109,7 +109,7 @@ public void build() { //System.out.println("PF::Constructed supportRingKinL13Top: " + supportRingKinL13Top.toString()); - UChannelL13 uChannelL14Top = new UChannelL14Top("support_top_L14", svtBox, alignmentCorrections, + UChannelL13 uChannelL14Top = new UChannelL14Top("support_top_L13", svtBox, alignmentCorrections, supportRingKinL13Top); surveyVolumes.add(uChannelL14Top); diff --git a/tracking/src/main/java/org/hps/recon/tracking/gbl/AlignmentStructuresBuilder.java b/tracking/src/main/java/org/hps/recon/tracking/gbl/AlignmentStructuresBuilder.java index 4775163547..ca91cdc053 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/gbl/AlignmentStructuresBuilder.java +++ b/tracking/src/main/java/org/hps/recon/tracking/gbl/AlignmentStructuresBuilder.java @@ -329,7 +329,7 @@ public AlignmentStructuresBuilder(List sensors) { aliVolumeList.add(alignable_volumes.get("module_L7b_halfmodule_stereo_hole_sensor0_AV")); aliVolumeList.add(alignable_volumes.get("module_L7b_halfmodule_stereo_slot_sensor0_AV")); rot = (alignable_volumes.get("module_L7b_halfmodule_stereo_hole_sensor0_AV")).getL2G().getRotation(); - AlignableVolume doublesensor_stereo_L7_Bot_AV = new AlignableVolume("doublesensor_stereo_L7_Bot_AV", aliVolumeList, rot, null,6); + AlignableVolume doublesensor_stereo_L7_Bot_AV = new AlignableVolume("doublesensor_stereo_L7_Bot_AV", aliVolumeList, rot, null,76); alignable_volumes.put("doublesensor_stereo_L7_Bot_AV",doublesensor_stereo_L7_Bot_AV); aliVolumeList.clear(); diff --git a/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java b/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java index 532b828b7e..a70df6fd2b 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java +++ b/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java @@ -68,7 +68,8 @@ public int getValue() { } }; - private GblTrajectory _traj; + private GblTrajectory _traj = null; + private GblTrajectoryJna _traj_jna = null; private double _chi2; private double _lost; private int _ndf; @@ -90,6 +91,14 @@ public FittedGblTrajectory(GblTrajectory traj, double chi2, int ndf, double lost _ndf = ndf; _lost = lost; } + + + public FittedGblTrajectory(GblTrajectoryJna traj_jna, double chi2, int ndf, double lost) { + _traj_jna = traj_jna; + _chi2 = chi2; + _ndf = ndf; + _lost = lost; + } /** * Find the index (or label) of the GBL point on the trajectory from the {@link GBLPOINT}. @@ -101,8 +110,12 @@ 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 if (point.compareTo(GBLPOINT.LAST) == 0) { + if (_traj != null) + gblPointIndex = _traj.getNumPoints(); + else + gblPointIndex = _traj_jna.getNumPoints(); + } else throw new RuntimeException("This GBL point " + point.toString() + "( " + point.name() + ") is not valid"); return gblPointIndex; @@ -135,7 +148,11 @@ public void getResults(GBLPOINT point, Vector locPar, SymMatrix locCov) { public void getResults(int iLabel, Vector locPar, SymMatrix locCov) { // Get the result from the trajectory - int ok = _traj.getResults(iLabel, locPar, locCov); + int ok = 0; + if (_traj != null) + ok = _traj.getResults(iLabel, locPar, locCov); + else + ok = _traj_jna.getResults(iLabel, locPar, locCov); // check that the fit was ok if (ok != 0) @@ -183,6 +200,10 @@ public GblTrajectory get_traj() { return _traj; } + public GblTrajectoryJna get_traj_jna(){ + return _traj_jna; + } + public double get_chi2() { return _chi2; } @@ -326,27 +347,50 @@ public Pair getCorrectedPerigeeParameters(HelicalTrac * @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)); + float[] lambdaKinks ; + double[] phiKinks ; + + if (_traj != null) { + + lambdaKinks = new float[_traj.getNumPoints() - 1]; + 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; } - oldPhi = newPhi; - oldLambda = newLambda; } - + else { + lambdaKinks = new float[_traj_jna.getNumPoints() - 1]; + phiKinks = new double[_traj_jna.getNumPoints() - 1]; + + double oldPhi = 0, oldLambda = 0; + for (int i = 0; i < _traj_jna.getNumPoints(); i++) { + _traj_jna.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); } diff --git a/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java b/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java index 384b50631a..89e4d8d795 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java +++ b/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java @@ -49,6 +49,7 @@ public class GBLOutputDriver extends Driver { private String outputPlots = "GBLplots_ali.root"; private String trackCollectionName = "GBLTracks"; private String trackResidualsRelColName = "TrackResidualsGBLRelations"; + private String dataRelationCollection = GBLKinkData.DATA_RELATION_COLLECTION; private List sensors = new ArrayList(); private double bfield; public boolean debug = false; @@ -72,7 +73,29 @@ public class GBLOutputDriver extends Driver { //Spacing between top and bottom in the 2D histos private int mod = 5; + + private double minMom = 1.; + private double maxMom = 6.; + + private int nHits = 6; + + + public void setDataRelationCollection (String val) { + dataRelationCollection = val; + } + + public void setNHits (int val ) { + nHits = val; + } + public void setMinMom (double val) { + minMom = val; + } + + public void setMaxMom (double val) { + maxMom = val; + } + //Override the Z of the target. public void setBsZ (double input) { bsZ = input; @@ -130,11 +153,10 @@ public void process(EventHeader event) { List tracks = event.get(Track.class, trackCollectionName); int TrackType = 0; - int nHits = 6; if (trackCollectionName.contains("Kalman") || trackCollectionName.contains("KF")) { TrackType = 1; - nHits = 12; + nHits = 2*nHits; //System.out.println("PF:: DEBUG :: Found Kalman Tracks in the event"); } @@ -173,16 +195,16 @@ public void process(EventHeader event) { Hep3Vector momentum = new BasicHep3Vector(trk.getTrackStates().get(0).getMomentum()); - if (momentum.magnitude() < 1.5) + if (momentum.magnitude() < minMom) continue; - if (momentum.magnitude() > 6) + if (momentum.magnitude() > maxMom) continue; - + //System.out.println("Track passed momentum"); TrackState trackState = trk.getTrackStates().get(0); - if (Math.abs(trackState.getTanLambda()) < 0.025) + if (Math.abs(trackState.getTanLambda()) < 0.015) continue; //System.out.println("Track passed tanLambda"); @@ -254,6 +276,13 @@ public void process(EventHeader event) { } private void doGBLkinks(Track trk, GenericObject kink, Map sensorNums) { + + if (kink == null) { + System.out.println("WARNING::Kink object is null"); + return; + } + + String vol = "_top"; int spacing = 0; if (trk.getTrackStates().get(0).getTanLambda() < 0) { @@ -335,8 +364,9 @@ private void doBasicGBLtrack(Track trk, Map sensorHits) isTop = "_top"; } + //There is a sign flip in the charge String charge = "_pos"; - if (trk.getCharge()<0) + if (trk.getCharge()>0) charge = "_neg"; @@ -405,8 +435,8 @@ private void doBasicGBLtrack(Track trk, Map sensorHits) //Get the PathToPlane BaseTrackState ts_bs = TrackUtils.getTrackExtrapAtVtxSurfRK(trackState,bFieldMap,0.,bsZ); - - + + //Get the track parameters wrt the beamline using helix double [] beamLine = new double [] {bsZ,0}; double [] helixParametersAtBS = TrackUtils.getParametersAtNewRefPoint(beamLine, trackState); @@ -424,6 +454,9 @@ private void doBasicGBLtrack(Track trk, Map sensorHits) FillGBLTrackPlot(trkpFolder+"z0_vs_bs_rk",isTop,charge,ts_bs.getZ0()); FillGBLTrackPlot(trkpFolder+"z0_vs_bs_extrap",isTop,charge,helixParametersAtBS[BaseTrack.Z0]); + + FillGBLTrackPlot(trkpFolder+"phi_vs_bs_extrap",isTop,charge,helixParametersAtBS[BaseTrack.PHI]); + //TH2D - Filling FillGBLTrackPlot(trkpFolder+"d0_vs_phi",isTop,charge,trackState.getPhi(),trackState.getD0()); FillGBLTrackPlot(trkpFolder+"d0_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackState.getD0()); @@ -831,7 +864,7 @@ private void setupPlots() { //TH1Ds aidaGBL.histogram1D(trkpFolder+"d0"+vol+charge,nbins_t,-5.0,5.0); aidaGBL.histogram1D(trkpFolder+"z0"+vol+charge,nbins_t,-1.3,1.3); - aidaGBL.histogram1D(trkpFolder+"phi"+vol+charge,nbins_t,-0.3,0.3); + aidaGBL.histogram1D(trkpFolder+"phi"+vol+charge,nbins_t,-0.06,0.06); aidaGBL.histogram1D(trkpFolder+"tanLambda"+vol+charge,nbins_t,-0.2,0.2); aidaGBL.histogram1D(trkpFolder+"p"+vol+charge,nbins_p,0.,pmax); aidaGBL.histogram1D(trkpFolder+"p7h"+vol+charge,nbins_p,0.,pmax); @@ -854,6 +887,7 @@ private void setupPlots() { aidaGBL.histogram1D(trkpFolder+"z0_vs_bs_rk"+vol+charge, 2*nbins_t, -z0bsmax, z0bsmax); aidaGBL.histogram1D(trkpFolder+"z0_vs_bs_extrap"+vol+charge, 2*nbins_t, -z0bsmax, z0bsmax); aidaGBL.histogram1D(trkpFolder+"z0_vs_bs"+vol+charge, 2*nbins_t, -z0bsmax, z0bsmax); + aidaGBL.histogram1D(trkpFolder+"phi_vs_bs_extrap"+vol+charge,2*nbins_t, -0.06,0.06); //TH2Ds diff --git a/tracking/src/main/java/org/hps/recon/tracking/gbl/GblTrajectoryJna.java b/tracking/src/main/java/org/hps/recon/tracking/gbl/GblTrajectoryJna.java index b804c9b35b..0b79a80cd3 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/gbl/GblTrajectoryJna.java +++ b/tracking/src/main/java/org/hps/recon/tracking/gbl/GblTrajectoryJna.java @@ -30,12 +30,14 @@ Pointer GblTrajectoryCtorPtrComposed(Pointer [] points_1, int npoints_1, double void GblTrajectory_fit(Pointer self, DoubleByReference Chi2, IntByReference Ndf, DoubleByReference lostWeight, char [] optionList, int aLabel); void GblTrajectory_addPoint(Pointer self, Pointer point); int GblTrajectory_isValid(Pointer self); + int GblTrajectory_getNumPoints(Pointer self); void GblTrajectory_printTrajectory(Pointer self, int level); void GblTrajectory_printData(Pointer self); void GblTrajectory_printPoints(Pointer self, int level); void GblTrajectory_getResults(Pointer self, int aSignedLabel, double[] localPar, int[] nLocalPar, double[] localCov, int[] sizeLocalCov); void GblTrajectory_milleOut(Pointer self, Pointer millebinary); + void GblTrajectory_getMeasResults(Pointer self, int aLabel, int[] numdata, double[] aResiduals, double[] aMeasErrors, double[] aResErrors, double[] aDownWeights); } @@ -116,13 +118,14 @@ public GblTrajectoryJna(List < Pair , Matrix> > PointsAndTrans ppoints_2, npoints_2, trafo_2); } + //to perform the full fit public void fit(DoubleByReference Chi2, IntByReference Ndf, DoubleByReference lostWeight, String optionList) { char[] c_optionList = optionList.toCharArray(); GblTrajectoryInterface.INSTANCE.GblTrajectory_fit(self, Chi2, Ndf, lostWeight, c_optionList,-999); } - + //To perform the fit removing a particular point public void fit(DoubleByReference Chi2, IntByReference Ndf, DoubleByReference lostWeight, String optionList, int aLabel) { char [] c_optionList = optionList.toCharArray(); @@ -138,6 +141,10 @@ public int isValid () { } + public int getNumPoints() { + return GblTrajectoryInterface.INSTANCE.GblTrajectory_getNumPoints(self); + } + public void printTrajectory(int level) { GblTrajectoryInterface.INSTANCE.GblTrajectory_printTrajectory(self,level); } @@ -150,8 +157,27 @@ public void printPoints(int level) { GblTrajectoryInterface.INSTANCE.GblTrajectory_printPoints(self,level); } + + public void getMeasResults(int aLabel, int numData[], List aResiduals,List aMeasErrors, List aResErrors, List aDownWeights) { + + double[] d_aResiduals = new double[2]; + double[] d_aMeasErrors = new double[2]; + double[] d_aResErrors = new double[2]; + double[] d_aDownWeights = new double[2]; + + GblTrajectoryInterface.INSTANCE.GblTrajectory_getMeasResults(self, aLabel, numData, d_aResiduals, d_aMeasErrors, d_aResErrors, d_aDownWeights); + + for (int i=0; i<2; i++) { + aResiduals.add(d_aResiduals[i]); + aMeasErrors.add(d_aMeasErrors[i]); + aResErrors.add(d_aResErrors[i]); + aDownWeights.add(d_aDownWeights[i]); + } + + } + //!! Only 5-localPar and 5x5 local Cov for the moment - public void getResults(int aSignedLabel, Vector localPar, SymMatrix localCov) { + public int getResults(int aSignedLabel, Vector localPar, SymMatrix localCov) { double[] d_localPar = new double[5]; double[] d_localCov = new double[25]; @@ -169,6 +195,8 @@ public void getResults(int aSignedLabel, Vector localPar, SymMatrix localCov) { localCov.set(i,j,d_localCov[i+5*j]); } } + + return 0; } diff --git a/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblTrajectoryCreator.java b/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblTrajectoryCreator.java index 9d9520346f..359eecdc59 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblTrajectoryCreator.java +++ b/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblTrajectoryCreator.java @@ -8,7 +8,7 @@ import hep.physics.vec.VecOp; import java.util.ArrayList; -import java.util.HashMap; +//import java.util.HashMap; import java.util.List; import java.util.Map; import java.util.logging.Logger; @@ -32,14 +32,16 @@ public HpsGblTrajectoryCreator() { // System.out.println("level " + LOGGER.getLevel().toString()); } - public List MakeGblPointsList(List hits, GBLBeamSpotPoint bs, double bfac) { + public List MakeGblPointsList(List hits, GBLBeamSpotPoint bs, double bfac, Map sensorMap, Map pathLengthMap) { - // Save the association between strip cluster and label, and between label and path length - useless here. - Map pathLengthMap = new HashMap(); - Map sensorMap = new HashMap(); + // Save the association between strip cluster and label, and between label and path length + //Map pathLengthMap = new HashMap(); + //Map sensorMap = new HashMap(); List listOfPoints = new ArrayList(); + + int iLabel = 0; double s = 0.; Matrix jacPointToPoint = new Matrix(5,5); diff --git a/tracking/src/main/java/org/hps/recon/tracking/gbl/SimpleGBLTrajAliDriver.java b/tracking/src/main/java/org/hps/recon/tracking/gbl/SimpleGBLTrajAliDriver.java index 358cc5e5f7..152c9eee83 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/gbl/SimpleGBLTrajAliDriver.java +++ b/tracking/src/main/java/org/hps/recon/tracking/gbl/SimpleGBLTrajAliDriver.java @@ -2,9 +2,11 @@ import java.util.ArrayList; import java.util.List; +import java.util.Collection; import java.util.Map; import java.util.HashMap; import static java.lang.Math.sqrt; +import org.apache.commons.math3.util.Pair; //import hep.physics.vec.VecOp; @@ -50,15 +52,18 @@ import org.lcsim.event.Track; import org.lcsim.event.base.BaseTrack; +//Fiducial cuts on the calorimeter cluster +import org.hps.record.triggerbank.TriggerModule; + //import org.lcsim.event.base.BaseTrackState; //import org.lcsim.fit.helicaltrack.HelixUtils; import org.lcsim.geometry.FieldMap; import org.lcsim.fit.helicaltrack.HelicalTrackFit; import org.lcsim.event.TrackerHit; -//import org.lcsim.event.base.BaseLCRelation; +import org.lcsim.event.base.BaseLCRelation; import org.lcsim.geometry.Detector; -//import org.lcsim.lcio.LCIOConstants; +import org.lcsim.lcio.LCIOConstants; import org.lcsim.util.Driver; @@ -128,18 +133,32 @@ public class SimpleGBLTrajAliDriver extends Driver { private boolean compositeAlign = false; private boolean constrainedFit = false; private boolean constrainedBSFit = false; - private double bsZ = -7.5; private boolean constrainedD0Fit = false; private boolean constrainedZ0Fit = false; + private boolean constrainedTanLFit = false; + private boolean constrainedPhi0Fit = false; + private double seedPhi0 = 1000000.0; // 1 mrad + private double beam_angle = 0.0305; // beam angle + private double seed_precision = 10000; // the constraint on q/p + private double bsZ = -7.5; + private double bsX = 0.0; + private double bsY = 0.0; private int trackSide = -1; private boolean doCOMAlignment = false; - private double seed_precision = 10000; // the constraint on q/p private double momC = 4.55; - private double minMom = 3; - private double maxMom = 6; - private double maxtanL = 0.025; - private int nHitsCut = 6; + private double minMom = -999; + private double maxMom = 999; + private double maxtanL = -999; + private double minPhi = -999; + private double maxPhi = 999; + private int nHitsCut = 4; private boolean useParticles = false; + private double clusterEnergyCutMin = -1; + private double clusterEnergyCutMax = 999; + private double posEoP = -1; + private double eleEoP = -1; + private boolean correctTrack = false; + private GblTrajectoryMaker _gblTrajMaker; @@ -152,6 +171,23 @@ public class SimpleGBLTrajAliDriver extends Driver { //Setting 0 is a single refit, 1 refit twice and so on.. private int gblRefitIterations = 5; + + + public void setCorrectTrack (boolean val) { + correctTrack = val; + } + public void setPosEoP (double val) { + posEoP = val; + } + + public void setEleEoP (double val) { + eleEoP = val; + } + + //Switch the COM alignment + public void setDoCOMAlignment(boolean val) { + doCOMAlignment = val; + } //Set -1 for no selection, 0-slot side tracks 1-hole side tracks public void setTrackSide (int side) { @@ -161,6 +197,15 @@ public void setTrackSide (int side) { public void setMomC (double val) { momC = val; } + + public void setClusterEnergyCutMin (double val) { + clusterEnergyCutMin = val; + } + + public void setClusterEnergyCutMax (double val) { + clusterEnergyCutMax = val; + } + public void setCompositeAlign (boolean val) { compositeAlign = val; } @@ -168,10 +213,27 @@ public void setCompositeAlign (boolean val) { public void setConstrainedFit (boolean val) { constrainedFit = val; } + + public void setConstrainedPhi0Fit (boolean val) { + constrainedPhi0Fit = val; + } + public void setSeedPhi0 (double val) { + seedPhi0 = val; + } + + public void setConstrainedTanLFit (boolean val) { + constrainedTanLFit = val; + } public void setBsZ(double val) { bsZ = val; } + public void setBsX(double val) { + bsX = val; + } + public void setBsY(double val) { + bsY = val; + } public void setSeedPrecision(double val) { seed_precision = val; @@ -291,6 +353,14 @@ public void setMaxtanL(double val) { maxtanL = val; } + public void setMinPhi(double val) { + minPhi = val; + } + + public void setMaxPhi(double val) { + maxPhi = val; + } + public void setMaxMom(double val) { maxMom = val; } @@ -326,7 +396,7 @@ protected void endOfData() { @Override protected void detectorChanged(Detector detector) { - + bFieldMap = detector.getFieldMap(); if (aidaGBL == null) @@ -429,12 +499,15 @@ protected void process(EventHeader event) { //If using Seed Tracker, get the hits from the event if (TrackType == 0) { - hitToStrips = TrackUtils.getHitToStripsTable(event,helicalTrackHitRelationsCollectionName); hitToRotated = TrackUtils.getHitToRotatedTable(event,rotatedHelicalTrackHitRelationsCollectionName); } - //Get the tracks from the particles + + // Create a mapping of matched Tracks to corresonding Clusters. + HashMap TrackClusterPairs = null; + + //Get the tracks from the particles - redundant, just loop on the map if (useParticles) { for (ReconstructedParticle particle : particles) { @@ -442,13 +515,9 @@ protected void process(EventHeader event) { continue; tracks.add(particle.getTracks().get(0)); } - } - - // Create a mapping of matched Tracks to corresonding Clusters. - HashMap TrackClusterPairs = null; - - if (useParticles) + TrackClusterPairs = GetClustersFromParticles(particles); + } //Loop over the tracks for (Track track : tracks) { @@ -497,9 +566,7 @@ protected void process(EventHeader event) { } // ask tracks only on a side - if (trackSide >= 0) - { - + if (trackSide >= 0) { if (trackSide > 1) { System.out.println("SimpleGBLTrajAliDriver:: wrong settings for track side selection"); continue; @@ -511,9 +578,8 @@ protected void process(EventHeader event) { else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) continue; } - } - + //Get the E/p from the cluster if (useParticles) { @@ -528,30 +594,130 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) double trackp = new BasicHep3Vector(trackState.getMomentum()).magnitude(); double e_o_p = em_cluster.getEnergy() / trackp; + //compute the correction momC = em_cluster.getEnergy(); + //Cluster energy cut + if (clusterEnergyCutMin > 0 || clusterEnergyCutMax < 999) { + if (!TriggerModule.inFiducialRegion(em_cluster)) + continue; + if (em_cluster.getEnergy() < clusterEnergyCutMin) + continue; + if (em_cluster.getEnergy() > clusterEnergyCutMax) + continue; + } + + double[] trk_prms = track.getTrackParameters(); - + if (trk_prms[BaseTrack.TANLAMBDA] > 0) { aidaGBL.histogram1D(eopFolder+"Ecluster_top").fill(momC); aidaGBL.histogram1D(eopFolder+"EoP_top").fill(momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_top").fill(trk_prms[BaseTrack.PHI],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_top").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_top").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + + //Track sign is flipped + if (track.getCharge() > 0) { + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_ele_top").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_ele_top").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_ele_top").fill(trk_prms[BaseTrack.PHI],momC/trackp); + } + else { + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_pos_top").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_pos_top").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_pos_top").fill(trk_prms[BaseTrack.PHI],momC/trackp); + } } else { aidaGBL.histogram1D(eopFolder+"Ecluster_bottom").fill(momC); aidaGBL.histogram1D(eopFolder+"EoP_bottom").fill(momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_bottom").fill(trk_prms[BaseTrack.PHI],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_bottom").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_bottom").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + + //Track sign is flipped + if (track.getCharge() > 0) { + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_ele_bottom").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_ele_bottom").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_ele_bottom").fill(trk_prms[BaseTrack.PHI],momC/trackp); + } + else { + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_pos_bottom").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_pos_bottom").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_pos_bottom").fill(trk_prms[BaseTrack.PHI],momC/trackp); + } + } - + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi").fill(trk_prms[BaseTrack.PHI],momC/trackp); + aidaGBL.histogram3D(eopFolder+"EoP_vs_tanLambda_phi").fill(trk_prms[BaseTrack.TANLAMBDA], + trk_prms[BaseTrack.PHI], + momC/trackp); - } + + + if (TriggerModule.inFiducialRegion(em_cluster)) { + + if (trk_prms[BaseTrack.TANLAMBDA] > 0) { + aidaGBL.histogram1D(eopFolder+"Ecluster_top_fid").fill(momC); + aidaGBL.histogram1D(eopFolder+"EoP_top_fid").fill(momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_top_fid").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_top_fid").fill(trk_prms[BaseTrack.PHI],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_top_fid").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + + //Track sign is flipped + if (track.getCharge() > 0) { + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_ele_top_fid").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_ele_top_fid").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_ele_top_fid").fill(trk_prms[BaseTrack.PHI],momC/trackp); + } + else { + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_pos_top_fid").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_pos_top_fid").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_pos_top_fid").fill(trk_prms[BaseTrack.PHI],momC/trackp); + } + } + else { + aidaGBL.histogram1D(eopFolder+"Ecluster_bottom_fid").fill(momC); + aidaGBL.histogram1D(eopFolder+"EoP_bottom_fid").fill(momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_bottom_fid").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_bottom_fid").fill(trk_prms[BaseTrack.PHI],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_bottom_fid").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + + + //Track sign is flipped + if (track.getCharge() > 0) { + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_ele_bottom_fid").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_ele_bottom_fid").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_ele_bottom_fid").fill(trk_prms[BaseTrack.PHI],momC/trackp); + } + else { + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP_pos_bottom_fid").fill(trackp,momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_pos_bottom_fid").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_pos_bottom_fid").fill(trk_prms[BaseTrack.PHI],momC/trackp); + } + } + + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_fid").fill(trk_prms[BaseTrack.TANLAMBDA],momC/trackp); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_fid").fill(trk_prms[BaseTrack.PHI],momC/trackp); + aidaGBL.histogram3D(eopFolder+"EoP_vs_tanLambda_phi_fid").fill(trk_prms[BaseTrack.TANLAMBDA], + trk_prms[BaseTrack.PHI], + momC/trackp); + + } + + } //Track biasing example //Re-fit the track? - //Only active for ST tracks - if (constrainedFit && TrackType == 0) { + //If momC < 0, only add a term in the covariance matrix to fix the momentum + if (constrainedFit) { + double momentum_param = 2.99792458e-04; //Get the track parameters double[] trk_prms = track.getTrackParameters(); @@ -561,18 +727,62 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) //Bias the FEEs to beam energy. Correct the curvature by projecting on X / Y plane double tanLambda = trk_prms[BaseTrack.TANLAMBDA]; double cosLambda = 1. / (Math.sqrt(1+tanLambda*tanLambda)); - double targetpT = momC * cosLambda; - //System.out.println("TargetpT: " + targetpT + " tanLambda = " + tanLambda); - double pt_bias = targetpT - pt; - //System.out.println("pT bias: " + pt_bias); - double corrected_pt = pt+pt_bias; - - double corrected_c = sign*(bfield*momentum_param)/(corrected_pt); - trk_prms[BaseTrack.OMEGA] = corrected_c; + + + if (momC > 0 ) { + double targetpT = momC * cosLambda; + //System.out.println("TargetpT: " + targetpT + " tanLambda = " + tanLambda); + double pt_bias = targetpT - pt; + //System.out.println("pT bias: " + pt_bias); + double corrected_pt = pt+pt_bias; + double corrected_c = sign*(bfield*momentum_param)/(corrected_pt); + trk_prms[BaseTrack.OMEGA] = corrected_c; + ((BaseTrack)track).setTrackParameters(trk_prms,bfield); + } + + //Correct positrons. getCharge is flipped + if (posEoP>0 && (track.getCharge() < 0)) { + + double p = sqrt(track.getPX()*track.getPX() + + track.getPY()*track.getPY() + + track.getPZ()*track.getPZ()); + + + //Provide the energy correction + double targetP = p * posEoP; + double targetPt = targetP * cosLambda; + double corrected_c = sign*(bfield*momentum_param)/(targetPt); + trk_prms[BaseTrack.OMEGA] = corrected_c; + ((BaseTrack)track).setTrackParameters(trk_prms,bfield); + } + + //Correct positrons + if (eleEoP>0 && (track.getCharge() > 0)) { + + double p = sqrt(track.getPX()*track.getPX() + + track.getPY()*track.getPY() + + track.getPZ()*track.getPZ()); + + + //Provide the energy correction + double targetP = p * posEoP; + double targetPt = targetP * cosLambda; + double corrected_c = sign*(bfield*momentum_param)/(targetPt); + trk_prms[BaseTrack.OMEGA] = corrected_c; + ((BaseTrack)track).setTrackParameters(trk_prms,bfield); + } + + }//constrained fit + + if (constrainedPhi0Fit) { + + double [] trk_prms = track.getTrackParameters(); + //Bias the track to target beam angle + trk_prms[BaseTrack.PHI] = beam_angle; ((BaseTrack)track).setTrackParameters(trk_prms,bfield); } - if (constrainedD0Fit && TrackType == 0) { + if (constrainedD0Fit) { double [] trk_prms = track.getTrackParameters(); //Bias the track double d0 = trk_prms[BaseTrack.D0]; @@ -580,12 +790,6 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) //double d0bias = targetd0 - d0; double d0bias = 0.; - if (trk_prms[BaseTrack.TANLAMBDA] > 0) { - d0bias = -0.887; - } - else { - d0bias = 1.58; - } double corrected_d0 = d0+d0bias; trk_prms[BaseTrack.D0] = corrected_d0; //System.out.println("d0" + d0); @@ -594,11 +798,11 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) ((BaseTrack)track).setTrackParameters(trk_prms,bfield); } - if (constrainedZ0Fit && TrackType == 0) { + if (constrainedZ0Fit) { double [] trk_prms = track.getTrackParameters(); //Bias the track double z0 = trk_prms[BaseTrack.Z0]; - double targetz0 = 0.; + double targetz0 = z0; double z0bias = targetz0 - z0; double corrected_z0 = z0+z0bias; trk_prms[BaseTrack.Z0] = corrected_z0; @@ -630,7 +834,7 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) HelicalTrackFit htf = TrackUtils.getHTF(track); double bfac = Constants.fieldConversion * bfield; - GBLBeamSpotPoint bsPoint = FormBSPoint(htf, bsZ); + GBLBeamSpotPoint bsPoint = FormBSPoint(htf, bsZ,bsX,bsY); DoubleByReference Chi2 = new DoubleByReference(0.); DoubleByReference lostWeight = new DoubleByReference(0.); @@ -638,12 +842,14 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) //Create a trajectory with the beamspot List points_on_traj = new ArrayList(); + Map sensorMap = new HashMap(); + Map pathLengthMap = new HashMap(); if (constrainedBSFit) { - points_on_traj = _hpsGblTrajCreator.MakeGblPointsList(trackGblStripClusterData, bsPoint, bfac); + points_on_traj = _hpsGblTrajCreator.MakeGblPointsList(trackGblStripClusterData, bsPoint, bfac,sensorMap, pathLengthMap); } else { - points_on_traj = _hpsGblTrajCreator.MakeGblPointsList(trackGblStripClusterData, null, bfac); + points_on_traj = _hpsGblTrajCreator.MakeGblPointsList(trackGblStripClusterData, null, bfac,sensorMap, pathLengthMap); } if (compositeAlign) { @@ -666,20 +872,28 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) GblTrajectoryJna trajForMPII = null; GblTrajectoryJna trajForMPII_unconstrained = new GblTrajectoryJna(points_on_traj,1,1,1); - if (!constrainedFit) { + //seed matrix q/p, yT', xT', xT, yT + SymMatrix seedPrecision = new SymMatrix(5); + + if (constrainedFit) + seedPrecision.set(0,0,seed_precision); + if (constrainedTanLFit) + seedPrecision.set(1,1,seed_precision); + if (constrainedPhi0Fit) + seedPrecision.set(2,2,seed_precision); + + //z0 constraint + //seedPrecision.set(4,4,seed_precision); + + //d0 constraint + //seedPrecision.set(3,3,1000000); + + if (!constrainedFit && !constrainedTanLFit && !constrainedPhi0Fit && !constrainedD0Fit && !constrainedZ0Fit) { trajForMPII = new GblTrajectoryJna(points_on_traj,1,1,1); } else { - //Seed constrained fit - SymMatrix seedPrecision = new SymMatrix(5); - //seed matrix q/p, yT', xT', xT, yT - - //q/p constraint - seedPrecision.set(0,0,seed_precision); - //d0 constraint - //seedPrecision.set(3,3,1000000); trajForMPII = new GblTrajectoryJna(points_on_traj,1,seedPrecision,1,1,1); } @@ -689,17 +903,131 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) //Fit the trajectory to get the Chi2 trajForMPII_unconstrained.fit(Chi2,Ndf, lostWeight,""); - + //Avoid to use tracks with terrible Chi2 if (Chi2.getValue() / Ndf.getValue() > writeMilleChi2Cut) continue; - trajForMPII.milleOut(mille); + if (writeMilleBinary) + trajForMPII.milleOut(mille); + if (correctTrack) { + + //Form the FittedGblTrajectory for the unconstrained fit + FittedGblTrajectory fitTraj = new FittedGblTrajectory(trajForMPII_unconstrained, Chi2.getValue(), Ndf.getValue(), lostWeight.getValue()); + + fitTraj.setSensorMap(sensorMap); + fitTraj.setPathLengthMap(pathLengthMap); + + Collection hth = track.getTrackerHits(); + List allHthList = TrackUtils.sortHits(hth); + Pair newTrack = MakeGblTracks.makeCorrectedTrack(fitTraj, TrackUtils.getHTF(track), allHthList, 0, bfield); + Track gblTrk = newTrack.getFirst(); + + if(computeGBLResiduals) { + + List b_residuals = new ArrayList(); + List b_sigmas = new ArrayList(); + List r_sensors = new ArrayList(); + + int numData[] = new int[1]; + Integer[] sensorsFromMapArray = fitTraj.getSensorMap().keySet().toArray(new Integer[0]); + for (int i_s = 0; i_s < sensorsFromMapArray.length; i_s++) { + int ilabel = sensorsFromMapArray[i_s]; + int mpid = fitTraj.getSensorMap().get(ilabel); + List aResiduals = new ArrayList(); + List aMeasErrors = new ArrayList(); + List aResErrors = new ArrayList(); + List aDownWeights = new ArrayList(); + trajForMPII_unconstrained.getMeasResults(ilabel,numData,aResiduals,aMeasErrors,aResErrors,aDownWeights); + if (numData[0]>1) { + System.out.printf("GBLRefitterDriver::WARNING::We have SCT sensors. Residuals dimensions should be <=1\n"); + } + for (int i=0; i u_aResiduals = new ArrayList(); + List u_aMeasErrors = new ArrayList(); + List u_aResErrors = new ArrayList(); + List u_aDownWeights = new ArrayList(); + + try { + //Fit removing the measurement + gbl_fit_traj_u.fit(Chi2_u,Ndf_u,lostWeight_u,"",ilabel); + gbl_fit_traj_u.getMeasResults(ilabel,u_numData,u_aResiduals,u_aMeasErrors,u_aResErrors,u_aDownWeights); + for (int i=0; i charges = new ArrayList(); + charges.add(""); + charges.add("_ele"); + charges.add("_pos"); + for (String vol : volumes) { - aidaGBL.histogram1D(eopFolder+"Ecluster"+vol,200,0,8); - aidaGBL.histogram1D(eopFolder+"EoP"+vol,200,0,4); + aidaGBL.histogram1D(eopFolder+"Ecluster"+vol,200,0,6); + aidaGBL.histogram1D(eopFolder+"EoP"+vol,200,0,2); + + double lmin = 0.; + double lmax = 0.08; + if (vol == "_bot") { + lmin = -0.08; + lmax = 0.; + } + + for (String charge : charges) { + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP"+charge+vol,200,0,6,200,0,2); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda"+charge+vol,200,lmin,lmax,200,0,2); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi"+charge+vol,200,-0.2,0.2,200,0,2); + } + + aidaGBL.histogram1D(eopFolder+"Ecluster"+vol+"_fid",200,0,5); + aidaGBL.histogram1D(eopFolder+"EoP"+vol+"_fid",200,0,2); + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP"+vol+"_fid",200,0,6,200,0,2); + + for (String charge : charges) { + aidaGBL.histogram2D(eopFolder+"EoP_vs_trackP"+charge+vol+"_fid",200,0,6,200,0,2); + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda"+charge+vol+"_fid",200,0.01,0.08,200,0,2); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi"+charge+vol+"_fid",200,-0.2,0.2,200,0,2); + } } + + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda",200,-0.1,0.1,200,0,2); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi",200,-0.2,0.2,200,0,2); + aidaGBL.histogram3D(eopFolder+"EoP_vs_tanLambda_phi",200,-0.08,0.08,200,-0.2,0.2,200,0,2); + + aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_fid",200,-0.1,0.1,200,0,2); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_fid",200,-0.2,0.2,200,0,2); + aidaGBL.histogram3D(eopFolder+"EoP_vs_tanLambda_phi_fid",200,-0.08,0.08,200,-0.2,0.2,200,0,2); + + - aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda",200,-0.1,0.1,200,0,4); } @@ -1266,15 +1630,15 @@ else if (TrackType == 1) { return null; } - GBLBeamSpotPoint FormBSPoint(HelicalTrackFit htf, double bsZ) { + GBLBeamSpotPoint FormBSPoint(HelicalTrackFit htf, double bsZ, double bsX, double bsY) { //Form the BeamsSpotPoint double [] center = new double[3]; double [] udir = new double[3]; double [] vdir = new double[3]; center[0] = bsZ; //Z - center[1] = 0.; //X - center[2] = 0.; //Y + center[1] = bsX; //X + center[2] = bsY; //Y udir[0] = 0; udir[1] = 0; @@ -1287,8 +1651,8 @@ GBLBeamSpotPoint FormBSPoint(HelicalTrackFit htf, double bsZ) { //Hard coded uncertainties double[] bserror = new double[2]; - bserror[0]=0.01; - bserror[1]=0.1; + bserror[0]=0.02; + bserror[1]=0.2; return GblUtils.gblMakeBsPoint(htf, center, udir, vdir, bserror); } @@ -1443,7 +1807,11 @@ private HashMap GetClustersFromParticles(List