lcsim/src/org/lcsim/util
diff -u -r1.4 -r1.5
--- OverlayDriver.java 12 Mar 2011 12:13:18 -0000 1.4
+++ OverlayDriver.java 29 Apr 2011 16:23:37 -0000 1.5
@@ -1,5 +1,10 @@
package org.lcsim.util;
+import hep.physics.particle.properties.ParticleType;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.BasicHepLorentzVector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.HepLorentzVector;
import hep.physics.vec.SpacePoint;
import java.io.File;
@@ -14,6 +19,7 @@
import org.apache.commons.math.distribution.DistributionFactory;
import org.apache.commons.math.distribution.PoissonDistribution;
import org.freehep.record.source.EndOfSourceException;
+import org.lcsim.detector.IDetectorElement;
import org.lcsim.event.EventHeader;
import org.lcsim.event.GenericObject;
import org.lcsim.event.HitWithPosition;
@@ -21,6 +27,7 @@
import org.lcsim.event.SimCalorimeterHit;
import org.lcsim.event.SimTrackerHit;
import org.lcsim.event.EventHeader.LCMetaData;
+import org.lcsim.event.base.BaseMCParticle;
import org.lcsim.event.base.BaseSimCalorimeterHit;
import org.lcsim.event.base.BaseSimTrackerHit;
import org.lcsim.geometry.Detector;
@@ -43,8 +50,8 @@
*/
public class OverlayDriver extends Driver {
- protected double c = 299.792458; // speed of light in mm/ns
- protected SpacePoint interactionPoint = new SpacePoint(); // assuming 0 0 0 as IP
+ static protected double c = 299.792458; // speed of light in mm/ns
+ static protected SpacePoint interactionPoint = new SpacePoint(); // assuming 0 0 0 as IP
protected DistributionFactory df;
protected double tofCaloOffset = -0.25; // tolerance for keeping calo hits: tof is calculated to center of cell, but interaction can happen earlier
protected int bunchCrossings;
@@ -62,6 +69,9 @@
protected PoissonDistribution backgroundDistribution;
protected Map<String, Double> readoutTimes;
protected Map<String, Map<Long,SimCalorimeterHit>> caloHitMap;
+ protected List<MCParticle> overlayMcParticles;
+ protected List<MCParticle> allMcParticles;
+ protected Map<MCParticle, MCParticle> mcParticleReferences;
// -------------------- Constructors --------------------
/**
@@ -80,7 +90,9 @@
overlayList = new ArrayList<Integer>();
readoutTimes = new HashMap<String, Double>();
caloHitMap = new HashMap<String, Map<Long,SimCalorimeterHit>>();
+ overlayMcParticles = new ArrayList<MCParticle>();
shuffleBackground = true;
+ mcParticleReferences = new HashMap<MCParticle, MCParticle>();
}
// -------------------- Steering Parameters --------------------
@@ -222,9 +234,12 @@
@Override
protected void process(EventHeader event) {
-
- // reset list of hit calorimeter cells
- caloHitMap.clear();
+
+ // reset the overlay mcParticle lists
+ overlayMcParticles.clear();
+
+ // always keep all mc particles from signal, no need for reset then
+ allMcParticles = event.getMCParticles();
// shift the signal event in time according to its BX
if (randomSignal) {
@@ -234,11 +249,10 @@
}
double signalTime = 0;
if (!signalAtZero) signalTime = signalBunchCrossing * bunchSpacing;
- if (this.getHistogramLevel() > HLEVEL_OFF) System.out.println("Moving signal event to BX: "+signalBunchCrossing);
+ if (this.getHistogramLevel() > HLEVEL_DEFAULT) System.out.println("Moving signal event to BX: "+signalBunchCrossing);
this.moveEventToTime(event, signalTime);
// building a list of all bunch crossings in this train
- overlayList.clear();
for (int bX = 0; bX != bunchCrossings; bX++) {
int nBackgroundEvts = 1;
if (overlayWeight != 0.) {
@@ -256,17 +270,30 @@
}
// shuffle the list
if (shuffleBackground) Collections.shuffle(overlayList, this.getRandom());
+ int bxCounter = 0;
for (int bX : overlayList) {
- if (this.getHistogramLevel() > HLEVEL_OFF) System.out.println("Overlaying BX "+bX);
+ bxCounter++;
double overlayTime = (bX - signalBunchCrossing) * bunchSpacing;
if (!signalAtZero) overlayTime = bX * bunchSpacing;
- EventHeader overlayEvent = this.getNextEvent(overlayEvents);
- if (this.getHistogramLevel() > HLEVEL_NORMAL) System.out.println("Memory free: "+100*Runtime.getRuntime().freeMemory()/Runtime.getRuntime().totalMemory()+(" %"));
+ if (this.getHistogramLevel() > HLEVEL_DEFAULT) {
+ System.out.println("Overlaying background event "+bxCounter+" / "+overlayList.size()+" at BX "+bX+" ("+overlayTime+"ns)");
+ }
+
+ if (this.getHistogramLevel() > HLEVEL_OFF) {
+ int toMB = 1024*1024;
+ long freeMemory = Runtime.getRuntime().freeMemory();
+ long totalMemory = Runtime.getRuntime().totalMemory();
+ System.out.println("Memory free: "+freeMemory/toMB+"MB / "+totalMemory/toMB+"MB ("+100*freeMemory/totalMemory+"%)");
+ }
+
+ EventHeader overlayEvent = getNextEvent(overlayEvents);
if (overlayEvent != null) {
if (event.getDetector().equals(overlayEvent.getDetector())) {
+ // clear the mc particle references, which are only needed within one background event
+ mcParticleReferences.clear();
this.mergeEvents(event, overlayEvent, overlayTime);
} else {
System.err.println("Unable to merge events simulated in different detectors");
@@ -275,6 +302,15 @@
System.err.println("Error reading from overlay event list");
}
}
+
+ // put the overlay mc particles into the event
+ int flags = event.getMetaData(event.getMCParticles()).getFlags();
+ flags = LCIOUtil.bitSet(flags, LCIOConstants.BITSubset, true);
+ event.put(mcOverlayName, overlayMcParticles, MCParticle.class, flags);
+
+ // reset all lists
+ caloHitMap.clear();
+ overlayList.clear();
}
@Override
@@ -296,7 +332,7 @@
* @param lcio The LCIO source
* @return The next event in the LCIO file
*/
- protected EventHeader getNextEvent(LCIOEventSource lcio) {
+ static protected EventHeader getNextEvent(LCIOEventSource lcio) {
EventHeader event = null;
try {
lcio.next();
@@ -320,7 +356,7 @@
* the position of the given hit along a straight line.
* @param hit
*/
- protected double getLosTof(HitWithPosition hit) {
+ static protected double getLosTof(HitWithPosition hit) {
return SpacePoint.distance(new SpacePoint(hit.getPositionVec()), interactionPoint)/c;
}
@@ -393,18 +429,15 @@
}
} else if (collectionType.isAssignableFrom(SimTrackerHit.class)) {
// SimTrackerHits
- movedCollection = new ArrayList<SimTrackerHit>();
- for (SimTrackerHit hit : event.get(SimTrackerHit.class, collectionName)) {
+ movedCollection = event.get(SimTrackerHit.class, collectionName);
+ for (SimTrackerHit hit : (List<SimTrackerHit>)movedCollection) {
// check if hit falls into relevant readout time window
double hitTime = hit.getTime() + time;
- double tofCorr = this.getLosTof(hit);
+ double tofCorr = getLosTof(hit);
if (timeWindow > 0) {
if (hitTime < signalTime + tofCorr + tofCaloOffset || hitTime > signalTime + tofCorr + timeWindow) continue;
}
- SimTrackerHit movedHit = new BaseSimTrackerHit(hit.getPosition(),
- hit.getdEdx(), hit.getMomentum(), hit.getPathLength(), hit.getTime()+time,
- hit.getCellID(), hit.getMCParticle(), hit.getMetaData(), hit.getDetectorElement());
- movedCollection.add(movedHit);
+ ((BaseSimTrackerHit)hit).setTime(hit.getTime()+time);
}
} else if (collectionType.isAssignableFrom(SimCalorimeterHit.class)) {
// SimCalorimeterHits
@@ -416,8 +449,10 @@
int nHitsMoved = 0;
for (SimCalorimeterHit hit : event.get(SimCalorimeterHit.class, collectionName)) {
// check if earliest energy deposit is later than relevant time window
- double tofCorr = this.getLosTof(hit);
+ double tofCorr = getLosTof(hit);
BaseSimCalorimeterHit movedHit = null;
+ nHitsMoved++;
+ if (this.getHistogramLevel() > HLEVEL_HIGH && nHitsMoved%100 == 0) System.out.print("Moved "+nHitsMoved+" / "+nSimCaloHits+" hits\n");
if (fullCaloProcessing) {
if (hit.getTime() > signalTime + tofCorr + timeWindow) continue;
nHitsMoved++;
@@ -467,9 +502,7 @@
movedHit.shiftTime(time);
}
movedCollection.add(movedHit);
- if (this.getHistogramLevel() > HLEVEL_HIGH && nHitsMoved%100 == 0) System.out.print("Moved "+nHitsMoved+" / "+nSimCaloHits+" hits\r");
}
- if (this.getHistogramLevel() > HLEVEL_HIGH) System.out.print("\n");
} else if (collectionType.isAssignableFrom(GenericObject.class)) {
// nothing to do for GenericObjects
return event.get(GenericObject.class, collectionName);
@@ -489,12 +522,6 @@
* @param overlayTime the time offset for the overlay event
*/
protected void mergeEvents(EventHeader event, EventHeader overlayEvent, double overlayTime) {
- // create a collection for the background mc particles
- if (!event.hasCollection(MCParticle.class, mcOverlayName)) {
- int flags = event.getMetaData(event.getMCParticles()).getFlags();
- flags = LCIOUtil.bitSet(flags, LCIOConstants.BITSubset, true);
- event.put(mcOverlayName, new ArrayList<MCParticle>(), MCParticle.class, flags);
- }
// need to copy list of collections to avoid concurrent modification
List<LCMetaData> overlayCollections = new ArrayList<LCMetaData>(overlayEvent.getMetaData());
@@ -504,6 +531,8 @@
this.mergeCollections(event.getMetaData((List)event.get(overlayCollectionName)), overlayCollection, overlayTime);
} else {
// event does not contain corresponding collection from overlayEvent, just put it there
+ // First move hits and apply timing cuts
+ List collection = this.moveCollectionToTime(overlayCollection, overlayTime);
this.putCollection(overlayCollection, (List)overlayEvent.get(overlayCollectionName), event);
}
}
@@ -515,20 +544,105 @@
* @param event
* @param particle
*/
- protected void addOverlayMcParticle(EventHeader event, MCParticle particle) {
- List<MCParticle> overlayMcParticles = event.get(MCParticle.class, mcOverlayName);
- List<MCParticle> mcParticles = event.getMCParticles();
- if (!overlayMcParticles.contains(particle)) {
- overlayMcParticles.add(particle);
- mcParticles.add(particle);
+ protected void addOverlayMcParticle(MCParticle particle) {
+ if (!mcParticleReferences.containsKey(particle)) {
+ // keep a copy of the mc particle instead of the original in order to close the background event
+ MCParticle mcp = copyMcParticle(particle);
+ mcParticleReferences.put(particle, mcp);
+ overlayMcParticles.add(mcp);
+ allMcParticles.add(mcp);
List<MCParticle> parents = particle.getParents();
+ // keep the parents as well and set the parent daughter relations
for (MCParticle parent : parents) {
- this.addOverlayMcParticle(event, parent);
+ this.addOverlayMcParticle(parent);
+ ((BaseMCParticle)mcParticleReferences.get(parent)).addDaughter(mcp);
}
}
}
/**
+ *
+ * @param mcParticle
+ * @return
+ */
+ static public MCParticle copyMcParticle(MCParticle mcParticle) {
+ Hep3Vector origin = new BasicHep3Vector(mcParticle.getOriginX(), mcParticle.getOriginY(), mcParticle.getOriginZ());
+ HepLorentzVector p = new BasicHepLorentzVector(mcParticle.getEnergy(), new double[] {mcParticle.getPX(), mcParticle.getPY(), mcParticle.getPZ()});
+ ParticleType ptype = mcParticle.getType().getParticlePropertyProvider().get(mcParticle.getPDGID());
+ int status = mcParticle.getGeneratorStatus();
+ double time = mcParticle.getProductionTime();
+ BaseMCParticle copyMcP = new BaseMCParticle(origin, p, ptype, status, time);
+ copyMcP.setMass(mcParticle.getMass());
+ copyMcP.setSimulatorStatus(mcParticle.getSimulatorStatus().getValue());
+ return copyMcP;
+ }
+
+ /**
+ *
+ * @param hit
+ * @param meta
+ * @return
+ */
+ protected SimTrackerHit copySimTrackerHit(SimTrackerHit hit, LCMetaData meta) {
+
+ double[] position = new double[3];
+ double[] momentum = new double[3];
+ double[] hitp = hit.getPosition();
+ double[] hitm = hit.getMomentum();
+ for (int i = 0; i != 3; i++ ) {
+ position[i] = hitp[i];
+ momentum[i] = hitm[i];
+ }
+ double dEdx = hit.getdEdx();
+ double pathLength = hit.getPathLength();
+ double time = hit.getTime();
+ int cellID = hit.getCellID();
+ MCParticle hitMC = hit.getMCParticle();
+ this.addOverlayMcParticle(hitMC);
+ MCParticle mcParticle = mcParticleReferences.get(hitMC);
+ IDetectorElement de = null;
+ try {
+ de = hit.getDetectorElement();
+ } catch (Exception e) {
+ // nothing to do
+ }
+
+ return new BaseSimTrackerHit(position, dEdx, momentum, pathLength, time, cellID, mcParticle, meta, de);
+ }
+
+ protected SimCalorimeterHit copySimCalorimeterHit(SimCalorimeterHit hit, LCMetaData meta, boolean hasPDG) {
+ long id = hit.getCellID();
+ double rawEnergy = hit.getRawEnergy();
+ double time = hit.getTime();
+ int nMCP = hit.getMCParticleCount();
+ Object[] mcparts = new Object[nMCP];
+ float[] energies = new float[nMCP];
+ float[] times = new float[nMCP];
+ int[] pdgs = null;
+ if (hasPDG) pdgs = new int[nMCP];
+ // fill arrays with values from hit
+ for (int i = 0; i != nMCP; i++) {
+ MCParticle hitMC = hit.getMCParticle(i);
+ this.addOverlayMcParticle(hitMC);
+ mcparts[i] = mcParticleReferences.get(hitMC);
+ energies[i] = (float)hit.getContributedEnergy(i);
+ times[i] = (float)hit.getContributedTime(i);
+ if (hasPDG) pdgs[i] = hit.getPDG(i);
+ }
+
+ BaseSimCalorimeterHit copyHit = new BaseSimCalorimeterHit(id, rawEnergy, time, mcparts, energies, times, pdgs);
+
+ copyHit.setMetaData(meta);
+ try {
+ copyHit.setDetectorElement(hit.getDetectorElement());
+ } catch (Exception e) {
+ // nothing to do
+ }
+
+ return copyHit;
+ }
+
+ /**
* Merges two collections and applies a time offset to all entries in
* the overlay collection.
* @param collection the collection where the overlay collection is merged into
@@ -550,7 +664,7 @@
List overlayEntries = this.moveCollectionToTime(overlayCollection, overlayTime);
//List overlayEntries = overlayCollection.getEvent().get(overlayCollectionType, overlayCollection.getName());
// Check if there are actually entries to overlay
- if (overlayEntries.size() == 0) return true;
+ if (overlayEntries.isEmpty()) return true;
EventHeader event = collection.getEvent();
if (collectionType.isAssignableFrom(MCParticle.class)) {
@@ -562,22 +676,23 @@
} else if (collectionType.isAssignableFrom(SimTrackerHit.class)) {
// SimTrackerHits: just append all hits from overlayEvents
- event.get(SimTrackerHit.class, collectionName).addAll(overlayEntries);
+ List<SimTrackerHit> signalTrackerHits = event.get(SimTrackerHit.class, collectionName);
// add contributing mc particles to lists
for (SimTrackerHit hit : (List<SimTrackerHit>)overlayEntries) {
- this.addOverlayMcParticle(event, hit.getMCParticle());
+ SimTrackerHit overlayHit = copySimTrackerHit(hit, collection);
+ signalTrackerHits.add(overlayHit);
}
} else if (collectionType.isAssignableFrom(SimCalorimeterHit.class)) {
// SimCalorimeterHits: need to merge hits in cells which are hit in both events
// check if map has already been filled
Map<Long,SimCalorimeterHit> hitMap;
- List<SimCalorimeterHit> hits = event.get(SimCalorimeterHit.class, collectionName);
+ List<SimCalorimeterHit> signalCaloHits = event.get(SimCalorimeterHit.class, collectionName);
if (!caloHitMap.containsKey(collectionName)) {
// build map of cells which are hit in signalEvent
hitMap = new HashMap<Long, SimCalorimeterHit>();
- for (SimCalorimeterHit hit : hits) {
+ for (SimCalorimeterHit hit : signalCaloHits) {
hitMap.put(hit.getCellID(), hit);
}
caloHitMap.put(collectionName, hitMap);
@@ -587,18 +702,13 @@
boolean hasPDG = LCIOUtil.bitTest(collection.getFlags(),LCIOConstants.CHBIT_PDG);
// loop over the hits from the overlay event
- for (SimCalorimeterHit overlayHit : (List<SimCalorimeterHit>)overlayEntries) {
- long cellID = overlayHit.getCellID();
-
- // add contributing mc particles to lists
- for (int i = 0; i != overlayHit.getMCParticleCount(); i++) {
- this.addOverlayMcParticle(event, overlayHit.getMCParticle(i));
- }
+ for (SimCalorimeterHit hit : (List<SimCalorimeterHit>)overlayEntries) {
+ long cellID = hit.getCellID();
if (hitMap.containsKey(cellID)) {
- SimCalorimeterHit hit = hitMap.get(overlayHit.getCellID());
- int nHitMcP = hit.getMCParticleCount();
- int nOverlayMcP = overlayHit.getMCParticleCount();
+ SimCalorimeterHit oldHit = hitMap.get(hit.getCellID());
+ int nHitMcP = oldHit.getMCParticleCount();
+ int nOverlayMcP = hit.getMCParticleCount();
int nMcP = nHitMcP + nOverlayMcP;
// arrays of mc particle contributions to the hit
Object[] mcpList = new Object[nMcP];
@@ -609,37 +719,42 @@
double rawEnergy = 0.;
// fill arrays with values from hit
for (int i = 0; i != nHitMcP; i++) {
- mcpList[i] = hit.getMCParticle(i);
- eneList[i] = (float)hit.getContributedEnergy(i);
- timeList[i] = (float)hit.getContributedTime(i);
- if (hasPDG) pdgList[i] = hit.getPDG(i);
+ mcpList[i] = oldHit.getMCParticle(i);
+ eneList[i] = (float)oldHit.getContributedEnergy(i);
+ timeList[i] = (float)oldHit.getContributedTime(i);
+ if (hasPDG) pdgList[i] = oldHit.getPDG(i);
rawEnergy += eneList[i];
}
// add values of overlay hit
for (int i = 0; i != nOverlayMcP; i++) {
int j = nHitMcP + i;
- mcpList[j] = overlayHit.getMCParticle(i);
- eneList[j] = (float)overlayHit.getContributedEnergy(i);
- timeList[j] = (float)overlayHit.getContributedTime(i);
- if (hasPDG) pdgList[j] = overlayHit.getPDG(i);
+ MCParticle hitMC = hit.getMCParticle(i);
+ if (!mcParticleReferences.containsKey(hitMC)) {
+ this.addOverlayMcParticle(hitMC);
+ }
+ mcpList[j] = mcParticleReferences.get(hitMC);
+ eneList[j] = (float)hit.getContributedEnergy(i);
+ timeList[j] = (float)hit.getContributedTime(i);
+ if (hasPDG) pdgList[j] = hit.getPDG(i);
rawEnergy += eneList[j];
}
// need to set time to 0 so it is recalculated from the timeList
- SimCalorimeterHit mergedHit = new BaseSimCalorimeterHit(hit.getCellID(),
+ SimCalorimeterHit mergedHit = new BaseSimCalorimeterHit(oldHit.getCellID(),
rawEnergy, 0., mcpList, eneList, timeList, pdgList);
mergedHit.setMetaData(collection);
// set detector elements for collections which actually implement this method
try {
- mergedHit.setDetectorElement(hit.getDetectorElement());
+ mergedHit.setDetectorElement(oldHit.getDetectorElement());
} catch (Exception e) {
// nothing to do
}
// replace old hit with merged hit
- hits.remove(hit);
- hits.add(mergedHit);
+ signalCaloHits.remove(oldHit);
+ signalCaloHits.add(mergedHit);
hitMap.put(cellID, mergedHit);
} else {
- hits.add(overlayHit);
+ SimCalorimeterHit overlayHit = copySimCalorimeterHit(hit, collection, hasPDG);
+ signalCaloHits.add(overlayHit);
hitMap.put(cellID, overlayHit);
}
