lcsim/src/org/lcsim/contrib/timb/mc/fast/cluster/ronan
diff -N MCFastRonan.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MCFastRonan.java 26 May 2006 07:22:06 -0000 1.1
@@ -0,0 +1,149 @@
+package org.lcsim.mc.fast.cluster.ronan;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Random;
+
+import org.lcsim.conditions.ConditionsEvent;
+import org.lcsim.conditions.ConditionsListener;
+import org.lcsim.conditions.ConditionsSet;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.Driver;
+
+
+
+
+/**
+ * Fast Monte Carlo cluster simulator
+ * @author M.Ronan Oct 2000 - Added "refined" cluster simulation
+ * @version
+ */
+public class MCFastRonan extends Driver implements ConditionsListener
+{
+ private final static int ElecID = 11;
+ private final static int NuEID = 12;
+ private final static int MuID = 13;
+ private final static int NuMuID = 14;
+ private final static int NuTauID = 16;
+ private final static int PhotonID = 22;
+ private final static int Neutralino1 = 1000022;
+ private final static int Neutralino2 = 1000023;
+ private final static int Neutralino3 = 1000025;
+ private final static int Neutralino4 = 1000035;
+ private static long iprint=0;
+
+
+ private ClusterResolutionTables clusterParm;
+
+ protected void process(EventHeader event)
+ {
+ if (clusterParm == null)
+ {
+ ConditionsSet conditions = getConditionsManager().getConditions("ClusterParameters");
+ conditions.addConditionsListener(this);
+ clusterParm = new ClusterResolutionTables(conditions);
+ }
+ List<Cluster> cl = new ArrayList<Cluster>();
+
+ boolean hist = getHistogramLevel() > 0;
+
+ List<MCParticle> particles = event.getMCParticles();
+ for (MCParticle p : particles)
+ {
+
+ // filter for FINALSTATE
+ if (p.getGeneratorStatus() != p.FINAL_STATE)
+ {
+ continue;
+ }
+
+ int PDGID = p.getPDGID();
+ int absPDGID = Math.abs(PDGID);
+ double charge = p.getCharge();
+
+
+ // filter neutrinos
+ boolean neutrino = absPDGID == NuEID || absPDGID == NuMuID || absPDGID == NuTauID
+ || absPDGID == Neutralino1 || absPDGID == Neutralino2 || absPDGID == Neutralino3 || absPDGID == Neutralino4;
+ if (neutrino)
+ {
+ continue;
+ }
+
+ double E = p.getEnergy();
+ double pt2 = p.getMomentum().magnitudeSquared()-p.getPZ()*p.getPZ();
+ double pt = Math.sqrt(pt2);
+ double ptot = p.getMomentum().magnitude();
+ double cosTheta = p.getPZ() / p.getMomentum().magnitude();
+
+ Random rand = getRandom();
+
+ // Photons
+ if (PDGID == PhotonID || PDGID == ElecID)
+ {
+ // within acceptance
+ // double thing = (1 - 1 / ( 1 + Math.exp( (E-clusterParm.getEMOnset())*clusterParm.getEMSharpness() ) ));
+ //if (rand.nextDouble() > thing)
+ //{
+ // continue;
+ //}
+ if (E < clusterParm.getEMOnset())
+ {
+ continue;
+ }
+ if (Math.abs(cosTheta) > clusterParm.getPolarEMOuter())
+ {
+ continue;
+ }
+
+ cl.add(new ReconEMCluster(clusterParm, rand, p, hist));
+
+ }
+
+ // Neutral hadrons
+ else if (PDGID != MuID)
+ {
+ // within acceptance
+
+ //double thing = (1 - 1 / ( 1 + Math.exp( (E-clusterParm.getHADOnset())*clusterParm.getHADSharpness() ) ));
+ //if (rand.nextDouble() > thing)
+ //{
+ //continue;
+ //}
+ if (E < clusterParm.getHADOnset())
+ {
+ continue;
+ }
+ if (Math.abs(cosTheta) > clusterParm.getPolarHADOuter())
+ {
+ continue;
+ }
+
+ cl.add(new ReconHADCluster(clusterParm, rand, p, hist));
+ }
+ }
+ double neg_energy_total = 0.;
+ double pos_energy_weight_total = 0.;
+ for (Cluster rcl : cl )
+ {
+ if( Math.abs(((ReconCluster)rcl).getMCParticle().getCharge()) > Double.MIN_VALUE ) continue ;
+ neg_energy_total+=((ReconCluster)rcl).getNegEnergy();
+ pos_energy_weight_total+=((ReconCluster)rcl).getNegEnergy() < 0. ? 0. : Math.min(((ReconCluster)rcl).getSigma(),((ReconCluster)rcl).getEnergy());
+ }
+ iprint++;
+ if(iprint < -20) System.out.println(" MCFast neg_energy_total= "+neg_energy_total+" pos_energy_weight_total= "+pos_energy_weight_total);
+ if( neg_energy_total < -Double.MIN_VALUE ) for (Cluster rcl : cl ) if( Math.abs(((ReconCluster)rcl).getMCParticle().getCharge()) < Double.MIN_VALUE && ((ReconCluster)rcl).getNegEnergy() >= 0. )
+ ((ReconCluster)rcl).adjustEnergy(neg_energy_total,pos_energy_weight_total);
+ event.put(EventHeader.CLUSTERS, cl, Cluster.class, 0);
+ }
+
+
+ public void conditionsChanged(ConditionsEvent event)
+ {
+ ConditionsSet conditions = getConditionsManager().getConditions("ClusterParameters");
+ clusterParm = new ClusterResolutionTables(conditions);
+ }
+}
+