lcsim/src/org/lcsim/digisim
diff -u -r1.2 -r1.3
--- RandomNoise.java 1 Jun 2005 17:56:19 -0000 1.2
+++ RandomNoise.java 27 Jun 2005 21:35:55 -0000 1.3
@@ -1,43 +1,36 @@
package org.lcsim.digisim;
-import org.apache.commons.math.distribution.DistributionFactory;
-import org.apache.commons.math.distribution.PoissonDistribution;
import java.util.Vector;
import java.util.Map;
+import org.apache.commons.math.distribution.DistributionFactory;
+import org.apache.commons.math.distribution.PoissonDistribution;
/**
- * A modifier for simulation of random noise in calorimeters
+ * An interface for simulation of random noise in calorimeters.
+ * Different noise parametrizations can be implemented extending this
+ * class.
*
* @author Guilherme Lima
- * @version $Id: RandomNoise.java,v 1.2 2005/06/01 17:56:19 lima Exp $
+ * @version $Id: RandomNoise.java,v 1.3 2005/06/27 21:35:55 lima Exp $
*/
-class RandomNoise extends AbstractCalHitModifier {
+abstract class RandomNoise extends AbstractCalHitModifier {
//*** Defining modifier interface ***
+ abstract public double getProbabilityAboveThreshold();
+ abstract public double drawRandomNoise();
+ abstract public double drawRandomNoiseAboveThreshold();
- /** initialize parameters */
+ /** Initialize parameters */
public void init(Vector<Double> params) {
_par = params;
- if(_debug>0) {
- int i=0;
- for( Double x : params ) {
- System.out.println(_name+".init(): _par["+(i++)+"] = "+x);
- }
- }
// decode parameters
- _amplNom = params.get(0);
- _amplSigma = params.get(1);
- _timeNom = params.get(2);
- _timeSigma = params.get(3);
-
- _meanNumNoisyCells = params.get(4);
- _system = params.get(5).intValue();
- _barend = params.get(6).intValue();
+ _system = params.get(0).intValue();
+ _barend = params.get(1).intValue();
+ _threshold = params.get(2);
- // poisson distribution generator
- DistributionFactory df = DistributionFactory.newInstance();
- _poisson = df.createPoissonDistribution( _meanNumNoisyCells );
+ _timeNom = params.get(3);
+ _timeSigma = params.get(4);
}
/**
@@ -49,48 +42,111 @@
System.out.println(" name="+_name+", debug="+_debug);
}
- _selector = _digitizer.getCellSelector();
+ if(_poisson==null) {
+ // estimate mean number of noise-only cells
+ if(_selector==null) _selector = _digitizer.getCellSelector();
+ int ncells = _selector.getCellCount();
+ double probAbove = getProbabilityAboveThreshold();
+ double meanNoisyCount = (probAbove * ncells);
+ System.out.println("probAboveThreshold() = "+probAbove
+ +", meanNoisyCount="+meanNoisyCount);
+
+ // poisson distribution generator
+ DistributionFactory df = DistributionFactory.newInstance();
+ _poisson = df.createPoissonDistribution( meanNoisyCount );
+ }
+
+ addNoiseToExistingHits(hitmap);
+ createNoiseOnlyHits(hitmap);
+ }
+
+ /**
+ * Adds noise to current hits. Noise is drawn according to some
+ * noise parametrization, see drawRandomNoise() method.
+ */
+ public void addNoiseToExistingHits(Map<Long,TempCalHit> hitmap) {
+ for( Long cellid : hitmap.keySet() ) {
+ TempCalHit ihit = hitmap.get(cellid);
+ double noise = drawRandomNoise();
+ double time = drawRandomTiming();
+ // add noise
+ double temp = ihit.getTotalEnergy();
+
+ if(temp+noise>0) {
+ ihit.addContribution( _noiseID, noise, time );
+ if(_debug>0) {
+ System.out.println(" add to existing hit: noise="+noise
+ +", cellid="+ cellid
+ +", energy="+ temp
+ +" --> "+ ihit.getTotalEnergy());
+ }
+ }
+ else {
+ // we might decide to remove these hits, even without a threshold
+ ihit.addContribution( _noiseID, noise, time );
+ if(_debug>0) {
+ System.out.println(" Hit became negative with noise!"
+ +" noise="+noise
+ +", cellid="+ cellid
+ +", energy="+ temp
+ +" --> "+ ihit.getTotalEnergy());
+ }
+ }
+ }
+ }
+
+ /** Timing parametrization: gaussian */
+ private double drawRandomTiming() {
+ return _timeNom + _timeSigma * _random.nextGaussian();
+ }
+
+ /** FIXME: add documentation here
+ */
+ public void createNoiseOnlyHits(Map<Long,TempCalHit> hitmap) {
+
+ if(_selector==null) _selector = _digitizer.getCellSelector();
_selector.setSystemBarrel(_system, _barend);
- // how many noisy cells
+ // get number of noisy cells (use poisson distribution)
int numNoisy = 0;
try {
- numNoisy=_poisson.inverseCumulativeProbability(_random.nextDouble());
+ numNoisy = _poisson.inverseCumulativeProbability(_random.nextDouble());
}
catch(Exception e) {
- System.out.println("RandomNoise: Exception caught: "+e.toString());
+ System.out.println("RandomNoise: Exception caught: "+e.toString());
}
+ if(_debug>0) {
+ System.out.println("RandomNoise: "+numNoisy+" noise only cells");
+ }
for(int i=0; i<numNoisy; ++i) {
// generate noise amplitude and timing
- double ampl = _amplNom + _amplSigma * _random.nextGaussian();
- double time = _timeNom + _timeSigma * _random.nextGaussian();
- if(_debug>1) {
- System.out.println("i="+i+"/"+numNoisy
- +", ampl="+ampl+", time="+time);
- }
+ double noise = drawRandomNoiseAboveThreshold();
+ double time = drawRandomTiming();
- // pick a noisy cell
- long chanid = _selector.getRandomCell();
- TempCalHit ahit = hitmap.get(chanid);
- if( ahit!=null ) {
- // add noise to existing cell
- System.out.println("RandomNoise: adding noise to existing cell"
- +", chanid="+ Long.toHexString(chanid)
- +", noiseID="+_noiseID);
- ahit.addContribution( _noiseID, ampl, time );
- }
- else {
- // create new hit and add it to hitmap
-// System.out.println("RandomNoise: creating new noise cell"
-// +", chanid="+ Long.toHexString(chanid)
-// +", noiseID="+_noiseID);
- ahit = new TempCalHit( chanid, _noiseID, ampl, time );
- hitmap.put( chanid, ahit );
+ // pick an empty cell to become noisy
+ long chanid = 0;
+ TempCalHit ahit = null;
+ for(;;) {
+ chanid = _selector.getRandomCell();
+ ahit = hitmap.get( chanid );
+ // break loop only for an empty hit
+ if( ahit==null ) break;
}
+
+ // create new hit and add it to hitmap
+ if(_debug>1) {
+ System.out.println("RandomNoise: Creating noise only cells"
+ +", cellid="+ Long.toHexString(chanid)
+ +", noise="+ noise+", time="+time);
+ }
+ ahit = new TempCalHit( chanid, _noiseID, noise, time );
+ hitmap.put( chanid, ahit );
}
}
+ //=== Standard modifier interface ===
+
/** debugging printout */
public void print() {
System.out.println( this.toString() );
@@ -105,60 +161,50 @@
System.out.println(ret);
}
- /** Returns a reference to a new instance */
- public RandomNoise newInstance(Digitizer digitizer,
- ModifierParameters modpars) {
- return new RandomNoise(digitizer, modpars);
- }
-
- public static RandomNoise instance() {
- return _me;
- }
+// /** Returns a reference to a new instance */
+// public RandomNoise newInstance(Digitizer digitizer,
+// ModifierParameters modpars) {
+// return new RandomNoise(digitizer, modpars);
+// }
/** private constructor, call newInstance() instead */
- private RandomNoise(Digitizer digitizer, ModifierParameters modpars) {
+ protected RandomNoise(Digitizer digitizer, ModifierParameters modpars) {
super(digitizer, modpars);
}
/** default constructor is only called for global instance */
- private RandomNoise() {
- super( null,
- new ModifierParameters( ModifierParameters.Type.RandomNoise.toString(),
- ModifierParameters.Type.RandomNoise,
- new Vector<Double>()
- )
- );
- }
-
- public static boolean registerMe() {
- String type = ModifierParameters.Type.RandomNoise.toString();
- AbstractCalHitModifier.registerModifier(type, _me);
- return true;
- }
+// protected RandomNoise() {
+// super( null,
+// new ModifierParameters(ModifierParameters.Type.RandomNoise.toString(),
+// ModifierParameters.Type.RandomNoise,
+// new Vector<Double>()
+// )
+// );
+// }
+
+// public static boolean registerMe() {
+// String type = ModifierParameters.Type.RandomNoise.toString();
+// AbstractCalHitModifier.registerModifier(type, _me);
+// return true;
+// }
//=== FIELDS ===
/** Parameters */
Vector<Double> _par;
/** system number to be used in hit encoding */
- private int _system;
+ protected int _system;
/** barrel/endcap number to be used in hit encoding */
- private int _barend;
+ protected int _barend;
/** mean number of noisy cells per event */
- private double _meanNumNoisyCells;
- /** mean amplitude */
- private double _amplNom;
- /** width used for smearing of amplitude */
- private double _amplSigma;
+ protected double _threshold;
/** mean amplitude */
- private double _timeNom;
+ protected double _timeNom;
/** width used for smearing of amplitude */
- private double _timeSigma;
+ protected double _timeSigma;
- /** Global instance, used only for registration */
- private static RandomNoise _me = new RandomNoise();
/** Random cell selector */
- private CellSelector _selector;
+ protected CellSelector _selector;
/** Poisson Distribution */
private PoissonDistribution _poisson;
/** noise ID, used in place of simhit cellID */