diff --git a/recon/src/main/java/org/hps/recon/vertexing/MultipleEventsVertexingDriver.java b/recon/src/main/java/org/hps/recon/vertexing/MultipleEventsVertexingDriver.java index 8ef4452dd4..84cc61dcd4 100644 --- a/recon/src/main/java/org/hps/recon/vertexing/MultipleEventsVertexingDriver.java +++ b/recon/src/main/java/org/hps/recon/vertexing/MultipleEventsVertexingDriver.java @@ -26,6 +26,7 @@ public final class MultipleEventsVertexingDriver extends Driver { private List accumulatedTracksBot = new ArrayList(); private List accumulatedBTracksBot = new ArrayList(); + private String trackCollectionName = "GBLTracks"; protected double[] beamSize = {0.001, 0.130, 0.050}; // rough estimate from harp scans during engineering run private double[] beamPositionToUse = new double[3]; @@ -35,13 +36,26 @@ public final class MultipleEventsVertexingDriver extends Driver { private AIDA aida; private String vtxFold="MultiEventVtx/"; private int _ntrks = 100; + private int _nhits = -1; private String outputPlots = "multiEvt.root"; + public void setNhits(int nhits) { + _nhits = nhits; + } + public void setTrackCollectionName(String val) { + trackCollectionName = val; + } + public void setNtrks( int ntrks) { _ntrks = ntrks; } + public void setVtxFold( String val) { + vtxFold = val; + } + + @Override protected void detectorChanged(Detector detector) { @@ -82,19 +96,23 @@ protected void process(EventHeader event) { } for (Track track : trackCollection) { - accumulatedTracks.add(track); - accumulatedBTracks.add(new BilliorTrack(track)); - - if (track.getTrackStates().get(0).getTanLambda() > 0) { - accumulatedTracksTop.add(track); - accumulatedBTracksTop.add(new BilliorTrack(track)); - } - - else { - accumulatedTracksBot.add(track); - accumulatedBTracksBot.add(new BilliorTrack(track)); - } + if (_nhits < 0 || track.getTrackerHits().size() == _nhits) { + + //accumulatedTracks.add(track); + accumulatedBTracks.add(new BilliorTrack(track)); + + if (track.getTrackStates().get(0).getTanLambda() > 0) { + //accumulatedTracksTop.add(track); + accumulatedBTracksTop.add(new BilliorTrack(track)); + } + + else { + //accumulatedTracksBot.add(track); + accumulatedBTracksBot.add(new BilliorTrack(track)); + } + + } } } @@ -116,7 +134,6 @@ protected void endOfData() { vtxFitter.setDebug(debug); vtxFitter.doBeamSpotConstraint(false); - //Randomize the tracks ? Better to do something more reproducible Collections.shuffle(accumulatedBTracks); @@ -124,6 +141,9 @@ protected void endOfData() { FitMultiVtx(accumulatedBTracksTop,vtxFitter,"_top"); FitMultiVtx(accumulatedBTracksBot,vtxFitter,"_bottom"); + + + } diff --git a/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java b/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java index 3238b43c93..4a0d578b1c 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java +++ b/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java @@ -1931,4 +1931,48 @@ public static boolean detectorElementContainsPoint(Hep3Vector trackPosition, Det return detectorElementContainsPoint( trackPosition, sensor,0.0); } + + //This methods checks if a track has only hole hits in the back of the detector + //return true if all back layers have hole hits, false if all back layers have slot hits + + public static boolean isHoleTrack(Track trk) { + + boolean holeTrack = false; + + TrackState trackState = trk.getTrackStates().get(0); + boolean isTop = true; + if (trackState.getTanLambda()<0) + isTop=false; + + //System.out.println("--------------"); + for (TrackerHit hit : trk.getTrackerHits()) { + + int stripLayer = ((HpsSiSensor) ((RawTrackerHit) hit.getRawHits().get(0)).getDetectorElement()).getLayerNumber(); + int hpslayer = (stripLayer + 1 ) / 2; + String side = ((HpsSiSensor) ((RawTrackerHit) hit.getRawHits().get(0)).getDetectorElement()).getSide(); + + if (isTop) { + if (hpslayer == 5 || hpslayer ==6 || hpslayer==7) + if (side=="ELECTRON") + holeTrack=true; + } + else { + if (hpslayer == 5 || hpslayer ==6 || hpslayer==7) + if (side=="ELECTRON") + holeTrack=true; + } + + String moduleName = ((HpsSiSensor) ((RawTrackerHit) hit.getRawHits().get(0)).getDetectorElement()).getName(); + int hpsmodule = ((HpsSiSensor) ((RawTrackerHit) hit.getRawHits().get(0)).getDetectorElement()).getModuleNumber(); + int MPID = ((HpsSiSensor) ((RawTrackerHit) hit.getRawHits().get(0)).getDetectorElement()).getMillepedeId(); + + //System.out.println("Hit on track :: " + moduleName); + //System.out.println("Layer="+hpslayer+" Module=" + hpsmodule +" MPID=" + MPID+" side="+side +" top=" +isTop); + + } + //System.out.println("Track is hole=" + holeTrack); + //System.out.println("=========="); + + return holeTrack; + } } 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 bb2e8f712e..981e57e16c 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 @@ -65,7 +65,7 @@ public class GBLOutputDriver extends Driver { String resFolder="/res/"; String hitFolder="/hit/"; private boolean b_doGBLkinks = true; - private boolean b_doDetailPlots = true; + private boolean b_doDetailPlots = false; //This should be moved to the GBL Refitter!!! //The field map for extrapolation @@ -157,7 +157,7 @@ public void process(EventHeader event) { i++; } - doBasicGBLtrack(trk); + doBasicGBLtrack(trk,sensorHits); doGBLresiduals(trk, sensorHits,event); //doMTresiduals(matchedTrack, sensorHits); if (b_doGBLkinks) @@ -220,7 +220,14 @@ private void FillGBLTrackPlot(String str, String isTop, String charge, double va aidaGBL.histogram2D(str+isTop+charge).fill(valX,valY); } - private void doBasicGBLtrack(Track trk) { + private void FillGBLTrackPlot(String str, String isTop, String charge, double valX, double valY, double valZ) { + aidaGBL.histogram3D(str+isTop).fill(valX,valY,valZ); + aidaGBL.histogram3D(str+isTop+charge).fill(valX,valY,valZ); + } + + + + private void doBasicGBLtrack(Track trk, Map sensorHits) { TrackState trackState = trk.getTrackStates().get(0); @@ -232,13 +239,14 @@ private void doBasicGBLtrack(Track trk) { if (trk.getType()==1 && trk.getTrackerHits().size() < 12) { return; } + + List missingHits; + missingHits = findMissingLayer(trk); if (trackState.getTanLambda() > 0) { isTop = "_top"; } - - - + String charge = "_pos"; if (trk.getCharge()<0) charge = "_neg"; @@ -252,9 +260,41 @@ private void doBasicGBLtrack(Track trk) { FillGBLTrackPlot(trkpFolder+"p",isTop,charge,trackp); if (trk.getTrackerHits().size()==7) FillGBLTrackPlot(trkpFolder+"p7h",isTop,charge,trackp); - + if (trk.getTrackerHits().size()==6) + FillGBLTrackPlot(trkpFolder+"p6h",isTop,charge,trackp); + if (trk.getTrackerHits().size()==5) + FillGBLTrackPlot(trkpFolder+"p5h",isTop,charge,trackp); + + if (TrackUtils.isHoleTrack(trk)) + FillGBLTrackPlot(trkpFolder+"p_hole",isTop,charge,trackp); + else + FillGBLTrackPlot(trkpFolder+"p_slot",isTop,charge,trackp); + + + //Momentum maps + FillGBLTrackPlot(trkpFolder+"p_vs_phi",isTop,charge,trackState.getPhi(),trackp); + FillGBLTrackPlot(trkpFolder+"p_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackp); + FillGBLTrackPlot(trkpFolder+"p_vs_phi_tanLambda",isTop,charge,trackState.getPhi(),trackState.getTanLambda(),trackp); + + double tanLambda = trackState.getTanLambda(); + double cosLambda = 1. / (Math.sqrt(1+tanLambda*tanLambda)); + + FillGBLTrackPlot(trkpFolder+"pT_vs_phi",isTop,charge,trackState.getPhi(),trackp*cosLambda); + FillGBLTrackPlot(trkpFolder+"pT_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackp*cosLambda); + + + if (trk.getTrackerHits().size()==6) + FillGBLTrackPlot(trkpFolder+"p_Missing1Hit",isTop,charge,missingHits.get(0),trackp); + + if (missingHits.size()==1 && missingHits.get(0)==7) + FillGBLTrackPlot(trkpFolder+"p_MissingLastLayer",isTop,charge,trackp); + + FillGBLTrackPlot(trkpFolder+"Chi2",isTop,charge,trk.getChi2()); - + FillGBLTrackPlot(trkpFolder+"Chi2_vs_p",isTop,charge,trackp,trk.getChi2()); + + + aidaGBL.histogram1D(trkpFolder+"nHits" + isTop).fill(trk.getTrackerHits().size()); aidaGBL.histogram1D(trkpFolder+"nHits" + isTop+charge).fill(trk.getTrackerHits().size()); @@ -507,6 +547,37 @@ private void doGBLresiduals(Track trk, Map sensorHits, } }//doGBLresiduals + private List findMissingLayer(Track trk) { + + List layers = new ArrayList(); + layers.add(1); + layers.add(2); + layers.add(3); + layers.add(4); + layers.add(5); + layers.add(6); + layers.add(7); + + List LayersOnTrack = new ArrayList(); + List missingHits = new ArrayList(); + + for (TrackerHit hit : trk.getTrackerHits()) { + + int stripLayer = ((HpsSiSensor) ((RawTrackerHit) hit.getRawHits().get(0)).getDetectorElement()).getLayerNumber(); + + int hpslayer = (stripLayer + 1 ) / 2; + LayersOnTrack.add(hpslayer); + } + + for (Integer layer : layers) { + if (!LayersOnTrack.contains(layer)) + missingHits.add(layer); + } + + return missingHits; + } + + private void setupPlots() { @@ -579,8 +650,8 @@ private void setupPlots() { if (sens.isTopLayer() && !sens.isAxial()) xmax = 0.001; } - aidaGBL.histogram1D(kinkFolder+"lambda_kink_" + sensor.getName(), 50, -xmax, xmax); - aidaGBL.histogram1D(kinkFolder+"phi_kink_" + sensor.getName(), 50, -xmax, xmax); + aidaGBL.histogram1D(kinkFolder+"lambda_kink_" + sensor.getName(), 250, -xmax, xmax); + aidaGBL.histogram1D(kinkFolder+"phi_kink_" + sensor.getName(), 250, -xmax, xmax); } aidaGBL.histogram2D(kinkFolder+"lambda_kink_mod",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5,nbins,-0.001,0.001); @@ -614,6 +685,11 @@ private void setupPlots() { 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); + aidaGBL.histogram1D(trkpFolder+"p6h"+vol+charge,nbins_p,0.,pmax); + aidaGBL.histogram1D(trkpFolder+"p5h"+vol+charge,nbins_p,0.,pmax); + aidaGBL.histogram1D(trkpFolder+"p_MissingLastLayer"+vol+charge,nbins_p,0.,pmax); + aidaGBL.histogram1D(trkpFolder+"p_hole"+vol+charge,nbins_p,0.,pmax); + aidaGBL.histogram1D(trkpFolder+"p_slot"+vol+charge,nbins_p,0.,pmax); aidaGBL.histogram1D(trkpFolder+"Chi2"+vol+charge,nbins_t*2,0,200); aidaGBL.histogram1D(trkpFolder+"nHits"+vol+charge,15,0,15); @@ -631,8 +707,9 @@ private void setupPlots() { //TH2Ds - + aidaGBL.histogram2D(trkpFolder+"d0_vs_phi"+vol+charge,nbins_t,-0.3,0.3,nbins_t,-5.0,5.0); + aidaGBL.histogram2D(trkpFolder+"Chi2_vs_p"+vol+charge,nbins_p,0.0,pmax,nbins_t*2,0,200); //aidaGBL.histogram2D("d0_vs_phi_bs"+vol+charge,nbins_t,-5.0,5.0,nbins_t,-0.3,0.3); aidaGBL.histogram2D(trkpFolder+"d0_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_t,-5.0,5.0); aidaGBL.histogram2D(trkpFolder+"d0_vs_p"+vol+charge, nbins_p,0.0,pmax,nbins_t,-5.0,5.0); @@ -641,6 +718,16 @@ private void setupPlots() { aidaGBL.histogram2D(trkpFolder+"z0bs_vs_p"+vol+charge,nbins_p,0.0,pmax,nbins_t,-z0bsmax,z0bsmax); aidaGBL.histogram2D(trkpFolder+"z0_vs_tanLambda"+vol+charge, nbins_t,-0.1,0.1,nbins_t,-z0max,z0max); aidaGBL.histogram2D(trkpFolder+"z0bs_vs_tanLambda"+vol+charge,nbins_t,-0.1,0.1,nbins_t,-z0bsmax,z0bsmax); + + aidaGBL.histogram2D(trkpFolder+"p_Missing1Hit"+vol+charge,8,0,8,nbins_p,0.0,pmax); + aidaGBL.histogram2D(trkpFolder+"p_vs_phi"+vol+charge, nbins_t,-0.3,0.3, nbins_p,0.,pmax); + aidaGBL.histogram2D(trkpFolder+"p_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_p,0.,pmax); + aidaGBL.histogram3D(trkpFolder+"p_vs_phi_tanLambda"+vol+charge, 50,-0.3,0.3,50,-0.2,0.2,100,0.,pmax); + + aidaGBL.histogram2D(trkpFolder+"pT_vs_phi"+vol+charge, nbins_t,-0.3,0.3, nbins_p,0.,pmax); + aidaGBL.histogram2D(trkpFolder+"pT_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_p,0.,pmax); + + if (b_doDetailPlots) { //TH2Ds - detail diff --git a/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLRefitterDriver.java b/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLRefitterDriver.java index b1c89dcf57..0ce51a773d 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLRefitterDriver.java +++ b/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLRefitterDriver.java @@ -24,6 +24,7 @@ import org.lcsim.event.RawTrackerHit; import org.lcsim.event.RelationalTable; import org.lcsim.event.Track; +import org.lcsim.event.base.BaseTrack; import org.lcsim.event.TrackerHit; import org.lcsim.event.base.BaseLCRelation; @@ -73,7 +74,7 @@ public class GBLRefitterDriver extends Driver { private List Alignabledes = new ArrayList(); private List sensors = new ArrayList(); private boolean usePoints = true; - + //Calculator for Frame to Frame derivatives private FrameToFrameDers f2fD = new FrameToFrameDers(); @@ -216,7 +217,7 @@ protected void process(EventHeader event) { List refittedTracks = new ArrayList(); List trackRelations = new ArrayList(); - + List kinkDataCollection = new ArrayList(); List kinkDataRelations = new ArrayList(); @@ -256,10 +257,18 @@ protected void process(EventHeader event) { if (momentum.magnitude() < 1) continue; - //At least 6 hits - if (gblTrk.getTrackerHits().size() < 6) - continue; + ////At least 6 hits + //if (gblTrk.getTrackerHits().size() < 6) + //continue; + + double[] trk_prms = track.getTrackParameters(); + double tanLambda = trk_prms[BaseTrack.TANLAMBDA]; + + //Align with tracks without hole on track + if ((tanLambda > 0 && track.getTrackerHits().size() < 5) || (tanLambda < 0 && track.getTrackerHits().size() < 5 )) + continue; + } //System.out.printf("gblTrkNDF %d gblTrkChi2 %f getMaxTrackChisq5 %f getMaxTrackChisq6 %f \n", gblTrk.getNDF(), gblTrk.getChi2(), cuts.getMaxTrackChisq(5), cuts.getMaxTrackChisq(6)); @@ -366,7 +375,7 @@ protected void process(EventHeader event) { int flag = 1 << LCIOConstants.TRBIT_HITS; event.put(outputCollectionName, refittedTracks, Track.class, flag); event.put(trackRelationCollectionName, trackRelations, LCRelation.class, 0); - + if (computeGBLResiduals) { event.put(trackResidualsColName, trackResidualsCollection, TrackResidualsData.class, 0); event.put(trackResidualsRelColName, trackResidualsRelations, LCRelation.class, 0); 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 a807c5bfb1..d4c1da7491 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 @@ -120,6 +120,7 @@ public class SimpleGBLTrajAliDriver extends Driver { private boolean constrainedD0Fit = false; private boolean constrainedZ0Fit = false; private boolean usePoints = true; + private int trackSide = -1; private GblTrajectoryMaker _gblTrajMaker; @@ -129,6 +130,11 @@ public class SimpleGBLTrajAliDriver extends Driver { //Setting 0 is a single refit, 1 refit twice and so on.. private int gblRefitIterations = 5; + //Set -1 for no selection, 0-slot side tracks 1-hole side tracks + public void setTrackSide (int side) { + trackSide = side; + } + public void setCompositeAlign (boolean val) { compositeAlign = val; } @@ -363,24 +369,50 @@ protected void process(EventHeader event) { continue; + if (enableAlignmentCuts) { - //At least 3.5 GeV + //Get the track parameters + double[] trk_prms = track.getTrackParameters(); + double tanLambda = trk_prms[BaseTrack.TANLAMBDA]; + + //Momentum cut: 3.8 - 5.2 Hep3Vector momentum = new BasicHep3Vector(track.getTrackStates().get(0).getMomentum()); - if (momentum.magnitude() < 3.5) + if (momentum.magnitude() < 3.8 || momentum.magnitude() > 5.2) continue; - //At least 6 hits - if (track.getTrackerHits().size() < 6) + + //Align with tracks with at least 6 hits + if ((tanLambda > 0 && track.getTrackerHits().size() < 6) || (tanLambda < 0 && track.getTrackerHits().size() < 6)) continue; + // ask tracks only on a side + if (trackSide >= 0) + { + + if (trackSide > 1) { + System.out.println("SimpleGBLTrajAliDriver:: wrong settings for track side selection"); + continue; + } + + if (trackSide == 0 && TrackUtils.isHoleTrack(track) ) + continue; + + else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) + continue; + + + } + } + - //Track biasing example + //Track biasing example + //Re-fit the track? if (constrainedFit) { double momentum_param = 2.99792458e-04; //Get the track parameters @@ -390,7 +422,8 @@ protected void process(EventHeader event) { int sign = trk_prms[BaseTrack.OMEGA] > 0. ? 1 : -1; //Bias the FEEs to beam energy. Correct the curvature by projecting on X / Y plane double tanLambda = trk_prms[BaseTrack.TANLAMBDA]; - double targetpT = 4.55 * Math.cos(Math.atan(tanLambda)); + double cosLambda = 1. / (Math.sqrt(1+tanLambda*tanLambda)); + double targetpT = 4.55 * cosLambda; //System.out.println("TargetpT: " + targetpT + " tanLambda = " + tanLambda); double pt_bias = targetpT - pt; //System.out.println("pT bias: " + pt_bias); @@ -793,115 +826,6 @@ protected void process(EventHeader event) { }//usePoints }//alitest - - - - /* - //System.out.printf("gblTrkNDF %d gblTrkChi2 %f getMaxTrackChisq5 %f getMaxTrackChisq6 %f \n", gblTrk.getNDF(), gblTrk.getChi2(), cuts.getMaxTrackChisq(5), cuts.getMaxTrackChisq(6)); - if (enableStandardCuts && (gblTrk.getChi2() > cuts.getMaxTrackChisq(gblTrk.getTrackerHits().size()))) - continue; - */ - - - //refittedTracks.add(gblTrk); - //trackRelations.add(new BaseLCRelation(track, gblTrk)); - //PF :: unused - //inputToRefitted.put(track, gblTrk); - //kinkDataCollection.add(newTrack.getSecond()); - //kinkDataRelations.add(new BaseLCRelation(newTrack.getSecond(), gblTrk)); - - - /* - - if (computeGBLResiduals) { - - GblTrajectory gbl_fit_trajectory = newTrackTraj.getSecond().get_traj(); - - List b_residuals = new ArrayList(); - List b_sigmas = new ArrayList(); - //List u_residuals = new ArrayList(); - //List u_sigmas = new ArrayList(); - List r_sensors = new ArrayList(); - - int numData[] = new int[1]; - //System.out.printf("Getting the residuals. Points on trajectory: %d \n",gbl_fit_trajectory.getNpointsOnTraj()); - //The fitted trajectory has a mapping between the MPID and the ilabel. Use that to get the MPID of the residual. - Integer[] sensorsFromMapArray = newTrackTraj.getSecond().getSensorMap().keySet().toArray(new Integer[0]); - //System.out.printf("Getting the residuals. Sensors on trajectory: %d \n",sensorsFromMapArray.length); - - //System.out.println("Check residuals of the original fit"); - //Looping on all the sensors on track - to get the biased residuals. - for (int i_s = 0; i_s < sensorsFromMapArray.length; i_s++) { - //Get the point label - int ilabel = sensorsFromMapArray[i_s]; - //Get the millepede ID - int mpid = newTrackTraj.getSecond().getSensorMap().get(ilabel); - List aResiduals = new ArrayList(); - List aMeasErrors = new ArrayList(); - List aResErrors = new ArrayList(); - List aDownWeights = new ArrayList(); - gbl_fit_trajectory.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(u_dVals,u_iVals,"",ilabel); - gbl_fit_traj_u.getMeasResults(ilabel,numData,u_aResiduals,u_aMeasErrors,u_aResErrors,u_aDownWeights); - for (int i=0; i