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 4289f22e5..60e72ff85 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 @@ -58,6 +58,7 @@ public class GBLOutputDriver extends Driver { private String trackCollectionName = "GBLTracks"; private String inputCollectionName = "FinalStateParticles_KF"; private String trackResidualsRelColName = "TrackResidualsGBLRelations"; + private String trackRelationsColName = "GBLToKFRelations"; private String dataRelationCollection = GBLKinkData.DATA_RELATION_COLLECTION; private List sensors = new ArrayList(); private double bfield; @@ -97,7 +98,11 @@ public class GBLOutputDriver extends Driver { private boolean useParticles = true; - + + + public void setDebug(boolean val) { + debug = val; + } public void setUseParticles(boolean val) { useParticles = val; @@ -218,6 +223,27 @@ public void process(EventHeader event) { // Create a mapping of matched Tracks to corresponding Clusters. HashMap TrackClusterPairs = new HashMap(); + + int TrackType = 0; + if (!useParticles) { + if (debug) + System.out.println("PF:: DEBUG :: NOT Using particles" + trackCollectionName); + if (trackCollectionName.contains("Kalman") || trackCollectionName.contains("KF")) { + TrackType = 1; + } + } + else { + if (debug) + System.out.println("PF:: DEBUG :: Using particles" + inputCollectionName); + if (inputCollectionName.contains("Kalman") || inputCollectionName.contains("KF")) { + + TrackType = 1 ; + } + + } + if (debug) + System.out.println("PF:: DEBUG :: Track Type=" + TrackType); + if (!useParticles) tracks = event.get(Track.class,trackCollectionName); @@ -228,9 +254,37 @@ public void process(EventHeader event) { continue; Track track = particle.getTracks().get(0); Cluster cluster = particle.getClusters().get(0); - tracks.add(track); - TrackClusterPairs.put(track,cluster); - } + + // If track is a Kalman Track, get the GBL refitted track from the relations + if (TrackType==1) { + RelationalTable trackRelationsTable = null; + if (event.hasCollection(LCRelation.class, trackRelationsColName)) { + trackRelationsTable = new BaseRelationalTable(RelationalTable.Mode.ONE_TO_ONE, RelationalTable.Weighting.UNWEIGHTED); + List trackRelation = event.get(LCRelation.class, trackRelationsColName); + for (LCRelation relation : trackRelation) { + if (relation != null && relation.getFrom() != null && relation.getTo() != null) { + trackRelationsTable.add(relation.getFrom(), relation.getTo()); + } + } + } // has the table + + Track gblTrack = (Track) trackRelationsTable.from(track); + + if (gblTrack == null) { + if (debug) + System.out.println("ERROR RETRIEVING GBLTRACK!"); + return; + } + tracks.add(gblTrack); + TrackClusterPairs.put(gblTrack,cluster); + }// it's a kalman track + else { + if (debug) + System.out.println("Particle not made with a kalman track. Adding track directly"); + tracks.add(track); + TrackClusterPairs.put(track,cluster); + } + } } @@ -241,15 +295,8 @@ public void process(EventHeader event) { } System.out.println(); */ - - int TrackType = 0; - if (trackCollectionName.contains("Kalman") || trackCollectionName.contains("KF")) { - TrackType = 1; - //System.out.println("PF:: DEBUG :: Found Kalman Tracks in the event"); - - } - - //System.out.println("Running on "+trackCollectionName); + + //System.out.println("Running on "+trackCollectionName); //RelationalTable trackMatchTable = null; //trackMatchTable = new BaseRelationalTable(RelationalTable.Mode.ONE_TO_ONE, RelationalTable.Weighting.UNWEIGHTED); @@ -308,9 +355,10 @@ public void process(EventHeader event) { if (Math.abs(trackState.getPhi()) > maxPhi) continue; - - //System.out.println("Track passed tanLambda"); - + + //System.out.println("Track passed tanLambda"); + + GenericObject gblKink = GBLKinkData.getKinkData(event, trk); //if (gblKink == null) { @@ -409,7 +457,7 @@ private void doEoPPlots(Track track, Cluster cluster) { aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_"+charge+"_"+vol).fill(phi,eop); aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda").fill(tanL,eop); - aidaGBL.histogram2D(eopFolder+"EoP_vs_phi").fill(tanL,eop); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi").fill(phi,eop); aidaGBL.histogram3D(eopFolder+"EoP_vs_tanLambda_phi").fill(tanL, phi, eop); @@ -429,12 +477,71 @@ private void doEoPPlots(Track track, Cluster cluster) { aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_"+charge+"_"+vol+"_fid").fill(phi,eop); aidaGBL.histogram2D(eopFolder+"EoP_vs_tanLambda_fid").fill(tanL,eop); - aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_fid").fill(tanL,eop); + aidaGBL.histogram2D(eopFolder+"EoP_vs_phi_fid").fill(phi,eop); aidaGBL.histogram3D(eopFolder+"EoP_vs_tanLambda_phi_fid").fill(tanL, phi, eop); + + // Cluster positions + + double clusterX = cluster.getPosition()[0]; + double clusterY = cluster.getPosition()[1]; + TrackState ts_ecal = TrackUtils.getTrackStateAtECal(track); + + if(ts_ecal == null){ + return; + } + + double[] ts_ecalPos = ts_ecal.getReferencePoint(); + double trkX = ts_ecalPos[1]; + double trkY = ts_ecalPos[2]; + + aidaGBL.histogram1D(eopFolder+"Xcluster_"+vol+"_fid").fill(clusterX); + aidaGBL.histogram1D(eopFolder+"Ycluster_"+vol+"_fid").fill(clusterY); + + aidaGBL.histogram1D(eopFolder+"trk_clu_resX_"+vol+"_fid").fill(trkX-clusterX); + aidaGBL.histogram1D(eopFolder+"trk_clu_resY_"+vol+"_fid").fill(trkY-clusterY); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vsX_"+vol+"_fid").fill(trkX,trkX-clusterX); + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vsY_"+vol+"_fid").fill(trkY,trkX-clusterX); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vsX_"+vol+"_fid").fill(trkX,trkY-clusterY); + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vsY_"+vol+"_fid").fill(trkY,trkY-clusterY); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vstrkP_"+vol+"_fid").fill(trackp,trkY-clusterY); + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vstrkP_"+vol+"_fid").fill(trackp,trkX-clusterX); + + aidaGBL.histogram2D(eopFolder+"trkY_vs_tanL_"+vol+"_fid").fill(tanL,trkY); + + + aidaGBL.histogram1D(eopFolder+"Xcluster_"+charge+"_"+vol+"_fid").fill(clusterX); + aidaGBL.histogram1D(eopFolder+"Ycluster_"+charge+"_"+vol+"_fid").fill(clusterY); + + aidaGBL.histogram1D(eopFolder+"trk_clu_resX_"+charge+"_"+vol+"_fid").fill(trkX-clusterX); + aidaGBL.histogram1D(eopFolder+"trk_clu_resY_"+charge+"_"+vol+"_fid").fill(trkY-clusterY); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vsX_"+charge+"_"+vol+"_fid").fill(trkX,trkX-clusterX); + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vsY_"+charge+"_"+vol+"_fid").fill(trkY,trkX-clusterX); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vsX_"+charge+"_"+vol+"_fid").fill(trkX,trkY-clusterY); + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vsY_"+charge+"_"+vol+"_fid").fill(trkY,trkY-clusterY); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vstrkP_"+charge+"_"+vol+"_fid").fill(trackp,trkY-clusterY); + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vstrkP_"+charge+"_"+vol+"_fid").fill(trackp,trkX-clusterX); + + aidaGBL.histogram2D(eopFolder+"trkY_vs_tanL_"+charge+"_"+vol+"_fid").fill(tanL,trkY); + + + // + + // As function of incident angle at ECAL, inclusive and in bin of momentum. + + + + + } @@ -580,7 +687,9 @@ private void doBasicGBLtrack(Track trk, Map sensorHits) // FillGBLTrackPlot(trkpFolder+"p_MissingLastLayer",isTop,charge,trackp); - FillGBLTrackPlot(trkpFolder+"Chi2",isTop,charge,trk.getChi2()); + + FillGBLTrackPlot(trkpFolder+"Chi2",isTop,charge,trk.getChi2()); + FillGBLTrackPlot(trkpFolder+"Chi2oNDF",isTop,charge,trk.getChi2() / trk.getNDF()); FillGBLTrackPlot(trkpFolder+"Chi2_vs_p",isTop,charge,trackp,trk.getChi2()); // deduce multiplication factor for ST-started GBL tracks @@ -595,7 +704,7 @@ private void doBasicGBLtrack(Track trk, Map sensorHits) Hep3Vector beamspot = CoordinateTransformations.transformVectorToDetector(TrackUtils.extrapolateHelixToXPlane(trackState, 0)); if (debug) - System.out.printf("beamspot %s transformed %s \n", beamspot.toString()); + System.out.printf("beamspot %s transformed \n", beamspot.toString()); FillGBLTrackPlot(trkpFolder+"trk_extr_or_x",isTop,charge,beamspot.x()); FillGBLTrackPlot(trkpFolder+"trk_extr_or_y",isTop,charge,beamspot.y()); @@ -770,15 +879,20 @@ private void doGBLresiduals(Track trk, Map sensorHits, trackResidualsTable.add(relation.getFrom(), relation.getTo()); } } + if (debug) + System.out.println("Loaded track Residuals Table"); } else { - //System.out.println("null TrackResidualsGBL Data Relations."); + if (debug) { + System.out.println("null TrackResidualsGBL Data Relations."); + } //Failed finding TrackResidualsGBL return; } GenericObject trackRes = (GenericObject) trackResidualsTable.from(trk); if (trackRes == null) { - //System.out.println("null TrackResidualsGBL Data."); + if (debug) + System.out.println("null TrackResidualsGBL Data."); return; } @@ -944,8 +1058,8 @@ private void setupEoPPlots() { aidaGBL.histogram1D(eopFolder+"Ecluster"+vol,200,0,6); aidaGBL.histogram1D(eopFolder+"EoP"+vol,200,0,2); - - double lmin = 0.; + + double lmin = 0.; double lmax = 0.08; if (vol == "_bot") { lmin = -0.08; @@ -961,11 +1075,55 @@ private void setupEoPPlots() { 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); + + + double cxrange = 20; + double cyrange = 20; + double ecalX = 400; + + aidaGBL.histogram1D(eopFolder+"Xcluster"+vol+"_fid",200,-ecalX,ecalX); + aidaGBL.histogram1D(eopFolder+"Ycluster"+vol+"_fid",200,-ecalX,ecalX); + aidaGBL.histogram1D(eopFolder+"trk_clu_resX"+vol+"_fid",200,-cxrange,cxrange); + aidaGBL.histogram1D(eopFolder+"trk_clu_resY"+vol+"_fid",200,-cyrange,cyrange); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vsX"+vol+"_fid",200,-ecalX,ecalX,200,-cxrange,cxrange); + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vsY"+vol+"_fid",200,-ecalX,ecalX,200,-cxrange,cxrange); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vsX"+vol+"_fid",200,-ecalX,ecalX,200,-cyrange,cyrange); + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vsY"+vol+"_fid",200,-ecalX,ecalX,200,-cyrange,cyrange); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vstrkP"+vol+"_fid",100,0.,5,200,-cyrange,cyrange); + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vstrkP"+vol+"_fid",100,0.,5,200,-cyrange,cyrange); + + aidaGBL.histogram2D(eopFolder+"trkY_vs_tanL"+vol+"_fid",200,-0.2,0.2,200,-100,100); + 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.histogram1D(eopFolder+"Xcluster"+charge+vol+"_fid",200,-ecalX,ecalX); + aidaGBL.histogram1D(eopFolder+"Ycluster"+charge+vol+"_fid",200,-ecalX,ecalX); + aidaGBL.histogram1D(eopFolder+"trk_clu_resX"+charge+vol+"_fid",200,-cxrange,cxrange); + aidaGBL.histogram1D(eopFolder+"trk_clu_resY"+charge+vol+"_fid",200,-cyrange,cyrange); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vsX"+charge+vol+"_fid",200,-ecalX,ecalX,200,-cxrange,cxrange); + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vsY"+charge+vol+"_fid",200,-ecalX,ecalX,200,-cxrange,cxrange); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vsX"+charge+vol+"_fid",200,-ecalX,ecalX,200,-cyrange,cyrange); + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vsY"+charge+vol+"_fid",200,-ecalX,ecalX,200,-cyrange,cyrange); + + aidaGBL.histogram2D(eopFolder+"trk_clu_resY_vstrkP"+charge+vol+"_fid",100,0.,5,200,-cyrange,cyrange); + aidaGBL.histogram2D(eopFolder+"trk_clu_resX_vstrkP"+charge+vol+"_fid",100,0.,5,200,-cyrange,cyrange); + + aidaGBL.histogram2D(eopFolder+"trkY_vs_tanL"+charge+vol+"_fid",200,-0.2,0.2,200,-100,100); + + + + } } @@ -1120,6 +1278,7 @@ private void setupPlots() { 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+"Chi2oNDF"+vol+charge,nbins_t*2,0,50); aidaGBL.histogram1D(trkpFolder+"nHits"+vol+charge,15,0,15); aidaGBL.histogram1D(trkpFolder+"trk_extr_or_x"+vol+charge,nbins_t,-3,3); aidaGBL.histogram1D(trkpFolder+"trk_extr_or_y"+vol+charge,nbins_t,-3,3); 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 cae334d3e..05f014e85 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 @@ -51,6 +51,7 @@ import org.lcsim.event.RelationalTable; import org.lcsim.event.Track; import org.lcsim.event.base.BaseTrack; +import org.lcsim.event.base.BaseTrackState; //Fiducial cuts on the calorimeter cluster import org.hps.record.triggerbank.TriggerModule; @@ -100,8 +101,7 @@ public class SimpleGBLTrajAliDriver extends Driver { private AIDA aidaGBL; String derFolder = "/gbl_derivatives/"; - String eopFolder = "/EoP/"; - + private String inputCollectionName = "MatchedTracks"; private String outputCollectionName = "GBLTracks"; private String trackRelationCollectionName = "MatchedToGBLTrackRelations"; @@ -130,6 +130,7 @@ public class SimpleGBLTrajAliDriver extends Driver { private List Alignabledes = new ArrayList(); private List sensors = new ArrayList(); private boolean debugAlignmentDs = false; + private boolean debug_ = false; private boolean compositeAlign = false; private boolean constrainedFit = false; private boolean constrainedBSFit = false; @@ -261,6 +262,9 @@ public void setDebugAlignmentDs (boolean val) { debugAlignmentDs = val; } + public void setDebug(boolean val) { + debug_ = val; + } public void setEnableAlignmentCuts (boolean val) { enableAlignmentCuts = val; @@ -440,19 +444,21 @@ protected void detectorChanged(Detector detector) { //Alignment Manager - Get the composite structures. IDetectorElement detectorElement = detector.getDetectorElement(); Alignabledes = detectorElement.findDescendants(AlignableDetectorElement.class); - - for (AlignableDetectorElement ade : Alignabledes) { - if (ade.getName().contains("alignable")) { - System.out.printf("Alignable Detector Elements informations: %s \n", ade.getName()); - //System.out.printf(((AlignableDetectorElement)ade).getlocalToGlobal().toString()+"\n"); - if (ade.getParent() != null) { - System.out.printf("The parent is: %s\n", ade.getParent().getName()); - } - else { - System.out.printf("No parent. \n"); - } - } - } + + if (debug_) { + for (AlignableDetectorElement ade : Alignabledes) { + if (ade.getName().contains("alignable")) { + System.out.printf("Alignable Detector Elements informations: %s \n", ade.getName()); + //System.out.printf(((AlignableDetectorElement)ade).getlocalToGlobal().toString()+"\n"); + if (ade.getParent() != null) { + System.out.printf("The parent is: %s\n", ade.getParent().getName()); + } + else { + System.out.printf("No parent. \n"); + } + } + } + } // Get the sensors subcomponents // This should be only HpsSiSensors sensors = detectorElement.findDescendants(SiSensor.class); @@ -460,7 +466,8 @@ protected void detectorChanged(Detector detector) { if (!doCOMAlignment) { //Assign the mothers to the sensors - //TODO FIX this part. For the moment the mother of the sensors are chosen by string parsing. + //TODO FIX this part. For the moment the mother of the sensors are chosen by string parsing. + MakeAlignmentTree("alignable_fullmodule"); //Dump the constrain file @@ -474,6 +481,7 @@ protected void detectorChanged(Detector detector) { @Override protected void process(EventHeader event) { + int runNumber = event.getRunNumber(); //Track collection @@ -517,10 +525,11 @@ protected void process(EventHeader event) { if (inputCollectionName.contains("Kalman") || inputCollectionName.contains("KF")) { TrackType = 1; - //System.out.println("PF:: DEBUG :: Found Kalman Tracks in the event"); - } + } - //System.out.println("DEBUG::Tom::Deduced a track type of "+TrackType); + if (debug_) { + System.out.println("DEBUG::Tom::Deduced a track type of "+TrackType); + } //If using Seed Tracker, get the hits from the event if (TrackType == 0) { @@ -621,7 +630,7 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) //Cluster cuts //FEE Clusters cuts - + /* if (clusterEnergyCutMin > 0 ) { if (clusters.size() != 1 ) continue; @@ -629,7 +638,7 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) if (clusters.get(0).getEnergy() < clusterEnergyCutMin ) continue; } - + */ } @@ -866,11 +875,18 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) Collection hth = track.getTrackerHits(); List allHthList = TrackUtils.sortHits(hth); Pair newTrack = MakeGblTracks.makeCorrectedTrack(fitTraj, TrackUtils.getHTF(track), allHthList, 0, bfield); - - // Track refit failed + // Track refit failed if (newTrack == null) continue; + + + // Extrapolate to the ECAL and make a new trackState there. + + Track gblTrk = newTrack.getFirst(); + BaseTrackState ts_ecal = new BaseTrackState(); + ts_ecal = TrackUtils.getTrackExtrapAtEcalRK(gblTrk, bFieldMap, runNumber); + gblTrk.getTrackStates().add(ts_ecal); //To make sure that the track fit converged if (writeMilleBinary) { @@ -883,7 +899,7 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) - Track gblTrk = newTrack.getFirst(); + //System.out.println("DEBUG::Tom::Correct GBL track has "+gblTrk.getTrackerHits().size()+" hits"); @@ -965,7 +981,10 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) //System.out.println("Refitted track chi2 " + gblTrk.getChi2()); + //Store the track and the relation refittedTracks.add(gblTrk); + trackRelations.add(new BaseLCRelation(gblTrk,track)); + kinkDataCollection.add(newTrack.getSecond()); kinkDataRelations.add(new BaseLCRelation(newTrack.getSecond(), gblTrk)); } @@ -992,15 +1011,16 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) // Put the tracks back into the event and exit 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); - } + event.put(GBLKinkData.DATA_COLLECTION, kinkDataCollection, GBLKinkData.class, 0); event.put(GBLKinkData.DATA_RELATION_COLLECTION, kinkDataRelations, LCRelation.class, 0); + } } } @@ -1318,16 +1338,17 @@ private void MakeAlignmentTree(String regEx) { }//loop on sensors - - for (SiSensor sensor : sensors) { - if (((HpsSiSensor)sensor).getAdeMother() != null) - System.out.printf("DEBUG::PF::MakeAlignmentTree sensor %s has mother %s \n", sensor.getName(), ((HpsSiSensor)sensor).getAdeMother().getName()); - } - - for (AlignableDetectorElement ade : Alignabledes) { - System.out.printf("DEBUG::PF::MakeAlignmentTree ade %s has children \n %s \n", ade.getName(), ade.getChildren().toString()); - - } + if (debug_) { + for (SiSensor sensor : sensors) { + if (((HpsSiSensor)sensor).getAdeMother() != null) + System.out.printf("DEBUG::PF::MakeAlignmentTree sensor %s has mother %s \n", sensor.getName(), ((HpsSiSensor)sensor).getAdeMother().getName()); + } + + for (AlignableDetectorElement ade : Alignabledes) { + System.out.printf("DEBUG::PF::MakeAlignmentTree ade %s has children \n %s \n", ade.getName(), ade.getChildren().toString()); + + } + } } //Matching by name