Skip to content

new track truth matching tool, and truth analysis driver #882

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
cbravo135 merged 1 commit into from
Apr 6, 2022
Merged

new track truth matching tool, and truth analysis driver #882

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
266 changes: 266 additions & 0 deletions recon/src/main/java/org/hps/recon/utils/TrackTruthMatcher.java
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,266 @@
package org.hps.recon.utils;

import java.util.ArrayList;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Set;

import org.lcsim.event.MCParticle;
import org.lcsim.event.RawTrackerHit;
import org.lcsim.event.RelationalTable;
import org.lcsim.event.SimTrackerHit;
import org.lcsim.event.Track;
import org.lcsim.event.TrackerHit;

import hep.physics.vec.BasicHep3Vector;

public class TrackTruthMatcher {

private MCParticle _mcp = null;
private int _nhits;
private int _ngoodhits;
private int _nbadhits;
private double _purity;
private Map<MCParticle, int[]> mcParticleStripHits = new HashMap<MCParticle, int[]>();
private Map<MCParticle, Set<Integer>> mcParticleHitsByLayer = new HashMap<MCParticle, Set<Integer>>();
private Map<Integer, List<RawTrackerHit>> _stripHitsOnLayerMap = new HashMap<Integer, List<RawTrackerHit>>();
private Map<Integer, Set<MCParticle>> _mcpsOnLayerMap = new HashMap<Integer, Set<MCParticle>>();
private Set<Integer> _layersOnTrack = new HashSet<Integer>();
private Map<RawTrackerHit, List<MCParticle>> _stripHitsToMCPsMap = new HashMap<RawTrackerHit, List<MCParticle>>();

public TrackTruthMatcher(Track track, RelationalTable rawtomc){
doAnalysis(track, rawtomc);
}

//Purity cut is a number less than 1
//enforces minimum ratio of bestMCPLayerHits / totalNLayerHits
//For example: 10 layer hits / 12 total layer hits is purity of 0.83
public TrackTruthMatcher(Track track, RelationalTable rawtomc, double purityCut){
doAnalysis(track, rawtomc, purityCut);
}

public TrackTruthMatcher(Track track, RelationalTable rawtomc, double purityCut, int nHitsRequired){
doAnalysis(track, rawtomc, purityCut, nHitsRequired);
}

private void doAnalysis(Track track, RelationalTable rawtomc){
getMCPsOnTrack(track, rawtomc);
matchTrackToMCP(this.mcParticleHitsByLayer);
}

private void doAnalysis(Track track, RelationalTable rawtomc, double purityCut){
getMCPsOnTrack(track, rawtomc);
matchTrackToMCP(this.mcParticleHitsByLayer);
if( _purity < purityCut)
this._mcp = null;
}

private void doAnalysis(Track track, RelationalTable rawtomc, double purityCut, int nHitsRequired){
getMCPsOnTrack(track, rawtomc);
matchTrackToMCP(this.mcParticleHitsByLayer);
if( _purity < purityCut || _ngoodhits < nHitsRequired)
this._mcp = null;
}

//Return MCParticle matched to Track
public MCParticle getMCParticle(){
return this._mcp;
}

//Return number of Layers with hits on track
public int getNHits() {
return _nhits;
}

//Return number of Layer hits for best MCP
public int getNGoodHits() {
return _ngoodhits;
}

//Return number of Layer hits from other MCPs (not best MCP)
public int getNBadHits() {
return _nbadhits;
}

//Return _ngoodhits / _nhits
public double getPurity() {
return _purity;
}

//Return Map of MCParticles and the layers they leave hits
public Map<MCParticle, Set<Integer>> getLayerHitsForAllMCPs(){
return this.mcParticleHitsByLayer;
}

//Return Map of number of strip hits left by all MCParticles
public Map<MCParticle, int[]> getNStripHitsForAllMCPs(){
return this.mcParticleStripHits;
}

//Return number of Strip Hits on layer
public List<RawTrackerHit> getStripHitsOnLayer(int layer){
return this._stripHitsOnLayerMap.get(layer);
}

//Return number of MCPs that hit layer
public Set<MCParticle> getMCPsOnLayer(int layer){
return this._mcpsOnLayerMap.get(layer);
}

//Return layers that Track leaves hits on
public Set<Integer> getLayersOnTrack(){
return this._layersOnTrack;
}

//Return MCParticles that left RawTrackerHit
public List<MCParticle> getMCPsOnRawTrackerHit(RawTrackerHit rawhit){
return this._stripHitsToMCPsMap.get(rawhit);
}

public Set<Integer> getBadHitLayers(){
Set<Integer> badLayers = new HashSet<Integer>();
Set<Integer> trackLayers = getLayersOnTrack();
Set<Integer> goodLayers = getLayerHitsForAllMCPs().get(this._mcp);
for(Integer layer : trackLayers){
if(!goodLayers.contains(layer)){
badLayers.add(layer);
}
}
return badLayers;
}

public Set<Integer> getGoodHitLayers(){
return getLayerHitsForAllMCPs().get(this._mcp);
}

public void getMCPsOnTrack(Track track, RelationalTable rawtomc){

double trackPmag = new BasicHep3Vector(track.getTrackStates().get(0).getMomentum()).magnitude();
int _nhits = track.getTrackerHits().size();

//loop over tracker hits
for (TrackerHit hit: track.getTrackerHits()) {
//Get layer by layer truth information
Map<Integer, Map<RawTrackerHit, List<MCParticle>>> hitTruth_byLayer = getMCParticlesOnTrackerHit(hit, rawtomc);
//loop over each layer that this trackerhit exists on
for(Map.Entry<Integer, Map<RawTrackerHit, List<MCParticle>>> entry : hitTruth_byLayer.entrySet()){
int layer = entry.getKey();
this._layersOnTrack.add(layer);
//loop over each strip hit that was clustered into this layer
//hit
List<RawTrackerHit> rawhitsOnLayer = new ArrayList<RawTrackerHit>();
Set<MCParticle> mcpsOnLayer = new HashSet<MCParticle>();
for(Map.Entry<RawTrackerHit, List<MCParticle>> subentry : entry.getValue().entrySet()){
rawhitsOnLayer.add(subentry.getKey());
//loop over each unique MCParticle that contributed to this
//strip hit
for(MCParticle particle : subentry.getValue()){
//Only count a MCParticle hit once per layer
Set<Integer> hitsOnLayer = new HashSet<Integer>();
hitsOnLayer.add(entry.getKey());
mcpsOnLayer.add(particle);
if(!mcParticleHitsByLayer.containsKey(particle)){
mcParticleHitsByLayer.put(particle, hitsOnLayer);
}
else{
Set<Integer> tmp = mcParticleHitsByLayer.get(particle);
tmp.add(entry.getKey());
mcParticleHitsByLayer.put(particle, tmp);
}
if(!mcParticleStripHits.containsKey(particle)){
mcParticleStripHits.put(particle, new int[1]);
mcParticleStripHits.get(particle)[0] = 0;
}
mcParticleStripHits.get(particle)[0]++;
}
}
this._mcpsOnLayerMap.put(layer, mcpsOnLayer);
this._stripHitsOnLayerMap.put(layer, rawhitsOnLayer);
}
}
}

//match track to MCP by matching to MCP that leaves hits on the most layers
private void matchTrackToMCP(Map<MCParticle, Set<Integer>> mcParticleHitsByLayer){

MCParticle bestMCP = null;
int bestNhits = 0;

for(Map.Entry<MCParticle, Set<Integer>> entry : mcParticleHitsByLayer.entrySet()){
int nhits = entry.getValue().size();
if(nhits > bestNhits){
bestNhits = nhits;
bestMCP = entry.getKey();
}
}

this._purity = (double) bestNhits/ (double) _layersOnTrack.size();
this._mcp = bestMCP;
this._ngoodhits = bestNhits;
this._nbadhits = _layersOnTrack.size() - bestNhits;
this._nhits = _layersOnTrack.size();
}

public Map<Integer, Map<RawTrackerHit, List<MCParticle>>> getMCParticlesOnTrackerHit(TrackerHit hit, RelationalTable rawtomc){

//Get all 1d Strip Hits on TrackerHit
List<RawTrackerHit> rawhits = hit.getRawHits();
//GBL TrackerHits are 3d clusters of 2d hits from axial and stereo pair
//Separate these 3d clusters into their independent layer 2d constituents
Set<Integer> layers = new HashSet<Integer>();
//Define which layers this TrackerHit is composed of
for(RawTrackerHit rawhit : rawhits){
int layer = rawhit.getLayerNumber();
layers.add(layer);
}

//Map StripHits to their respective layers
Map<Integer, Map<RawTrackerHit, List<MCParticle>>> rawHitsMCPMap_byLayer = new HashMap<Integer,Map<RawTrackerHit, List<MCParticle>>>();
//Loop over Striphits, layer by layer
for(Integer layer : layers){
Map<RawTrackerHit, List<MCParticle>> rawHitsMCPMap = new HashMap<RawTrackerHit, List<MCParticle>>();
for(RawTrackerHit rawhit : rawhits){
if(rawhit.getLayerNumber() != layer)
continue;
//Get list of unique MCPs that make up this Strip hit
List<MCParticle> rawhitMCPs = getMCParticlesOnRawTrackerHit(rawhit, rawtomc);
rawHitsMCPMap.put(rawhit,rawhitMCPs);
}
rawHitsMCPMap_byLayer.put(layer, rawHitsMCPMap);
}

return rawHitsMCPMap_byLayer;
}

public List<MCParticle> getMCParticlesOnRawTrackerHit(RawTrackerHit rawhit, RelationalTable rawtomc){

Set<MCParticle> mcps = new HashSet<MCParticle>();
List<MCParticle> mcpList = new ArrayList<MCParticle>();
Set<SimTrackerHit> simhits = rawtomc.allFrom(rawhit);
//loop over simhits on rawhit
for(SimTrackerHit simhit : simhits){
//get mcp that left simhit
MCParticle particle = simhit.getMCParticle();
int pdgid = particle.getPDGID();
//only add electron or positron mcps
if(Math.abs(pdgid) != 11)
continue;
/*
if(particle.getOriginZ() > 0)
continue;
*/
mcps.add(particle);
}

//Convert set to list, so that list contains only 1 entry per MCP
mcpList.addAll(mcps);
this._stripHitsToMCPsMap.put(rawhit, mcpList);

return mcpList;
}
}



Loading