Print

Print


Commit in lcsim/src/org/lcsim/contrib/timb/mc/fast/cluster/ronan on MAIN
MCFastRonan.java+149added 1.1
fastmc modifications Oct 2005 - Feb 2006 

lcsim/src/org/lcsim/contrib/timb/mc/fast/cluster/ronan
MCFastRonan.java added at 1.1
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);
+    }
+}
+
CVSspam 0.2.8