Print

Print


Commit in lcsim/src/org/lcsim/mc/fast on MAIN
cluster/ronan/ClusterResolutionTables.java-141.2 -> 1.3
             /MCFastRonan.java+15-111.6 -> 1.7
             /ReconCluster.java+17-161.3 -> 1.4
reconstructedparticle/MCFastReconstructedParticleDriver.java+35-81.7 -> 1.8
tracking/MCFastTracking.java+18-11.5 -> 1.6
        /ReconTrack.java+15-61.5 -> 1.6
        /SmearTrack.java+11.1 -> 1.2
MCFast.java+7-21.2 -> 1.3
+108-58
8 modified files
Daniel's updates.

lcsim/src/org/lcsim/mc/fast/cluster/ronan
ClusterResolutionTables.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- ClusterResolutionTables.java	15 Jul 2005 22:09:27 -0000	1.2
+++ ClusterResolutionTables.java	9 Aug 2005 18:34:45 -0000	1.3
@@ -8,7 +8,6 @@
     private double EMConstantTerm;
     private double EMPositionError;
     private double EMResolution;
-    private double EMmin;
     private double EMOnset;
     private double EMSharpness;
     
@@ -16,7 +15,6 @@
     private double HADConstantTerm;
     private double HADPositionError;
     private double HADResolution;
-    private double HADmin;
     private double HADOnset;
     private double HADSharpness;
     
@@ -27,7 +25,6 @@
     
     ClusterResolutionTables(ConditionsSet set)
     {
-        EMmin = set.getDouble("EMmin");
         EMOnset = set.getDouble("EMOnset");
         EMSharpness = set.getDouble("EMSharpness");
         PolarEMInner = set.getDouble("PolarEMInner");
@@ -38,7 +35,6 @@
         EMPositionError = set.getDouble("EMPositionError");
         EMAlignmentError = set.getDouble("EMAlignmentError");
         
-        HADmin = set.getDouble("HADmin");
         HADOnset = set.getDouble("HADOnset");
         HADSharpness = set.getDouble("HADSharpness");
         PolarHADInner = set.getDouble("PolarHADInner");
@@ -70,11 +66,6 @@
         return EMResolution;
     }
     
-    public double getEMmin()
-    {
-        return EMmin;
-    }
-    
     public double getEMOnset()
     {
         return EMOnset;
@@ -105,11 +96,6 @@
         return HADResolution;
     }
     
-    public double getHADmin()
-    {
-        return HADmin;
-    }
-    
     public double getHADOnset()
     {
         return HADOnset;

lcsim/src/org/lcsim/mc/fast/cluster/ronan
MCFastRonan.java 1.6 -> 1.7
diff -u -r1.6 -r1.7
--- MCFastRonan.java	4 Aug 2005 02:23:23 -0000	1.6
+++ MCFastRonan.java	9 Aug 2005 18:34:45 -0000	1.7
@@ -9,10 +9,14 @@
 import org.lcsim.conditions.ConditionsEvent;
 import org.lcsim.conditions.ConditionsListener;
 import org.lcsim.conditions.ConditionsSet;
+import org.lcsim.event.MCParticle;
 import org.lcsim.event.Cluster;
 import org.lcsim.event.EventHeader;
 import org.lcsim.util.Driver;
 
+
+
+
 /**
  * Fast Monte Carlo cluster simulator
  * @author M.Ronan  Oct 2000 - Added "refined" cluster simulation
@@ -34,19 +38,20 @@
             clusterParm = new ClusterResolutionTables(conditions);
         }
         List cl = new ArrayList();
-        List p1 = new ArrayList();
+
         boolean hist = getHistogramLevel() > 0;
         
-        for (Particle p : event.getMCParticles())
-        {     
+        List<MCParticle> particles = event.getMCParticles();
+        for (MCParticle p : particles)
+        {
+            
             // filter for FINALSTATE
             if (p.getGeneratorStatus() != p.FINAL_STATE)
             {
                 continue;
             }
             
-            ParticleType ptype = p.getType();
-            int PDGID = ptype.getPDGID();
+            int PDGID = p.getPDGID();
             double charge = p.getCharge();
             
             // filter neutrinos
@@ -57,10 +62,10 @@
             }
             
             double E = p.getEnergy();
-            double pt2 = (p.getPX() * p.getPX()) + (p.getPY() * p.getPY());
+            double pt2 = p.getMomentum().magnitudeSquared()-p.getPZ()*p.getPZ();
             double pt = Math.sqrt(pt2);
-            double ptot = Math.sqrt(pt2 + (p.getPZ() * p.getPZ()));
-            double cosTheta = p.getPZ() / ptot;
+            double ptot = p.getMomentum().magnitude();
+            double cosTheta = p.getPZ() / p.getMomentum().magnitude();
             
             Random rand = getRandom();
             
@@ -68,7 +73,6 @@
             if (PDGID == PhotonID)
             {
                 // within acceptance
-                double t = clusterParm.getEMOnset();
                 double thing = (1 - 1 / ( 1 + Math.exp( (E-clusterParm.getEMOnset())*clusterParm.getEMSharpness() ) ));
                 if (rand.nextDouble() > thing)
                 {
@@ -91,11 +95,11 @@
                 double thing = (1 - 1 / ( 1 + Math.exp( (E-clusterParm.getHADOnset())*clusterParm.getHADSharpness() ) ));
                 if (rand.nextDouble() > thing)
                 {
-                    continue;
+                    //continue;
                 }
                 if (Math.abs(cosTheta) > clusterParm.getPolarHADOuter())
                 {
-                    continue;
+                    //continue;
                 }
                 
                 cl.add(new ReconHADCluster(clusterParm, rand, p, hist));

lcsim/src/org/lcsim/mc/fast/cluster/ronan
ReconCluster.java 1.3 -> 1.4
diff -u -r1.3 -r1.4
--- ReconCluster.java	4 Aug 2005 02:23:23 -0000	1.3
+++ ReconCluster.java	9 Aug 2005 18:34:45 -0000	1.4
@@ -13,7 +13,6 @@
    protected double a = 0;
    protected double b = 0;
    protected double c = 0;
-   protected double cosTh;
    protected double d = 0;
    protected double energy;
    protected double phi;
@@ -49,10 +48,11 @@
    protected void smearEnergy(Random rand, double E, boolean hist)
    {
       double sigma = ((a / Math.sqrt(E)) + b) * E;
-      energy = E + (sigma * rand.nextGaussian());
-      if (energy < mcp.getMass())
+      
+      energy = 0;
+      while (energy <= mcp.getMass())
       {
-         energy = mcp.getMass() + 0.05;
+         energy = E + (sigma * rand.nextGaussian());
       }
    }
 
@@ -63,17 +63,14 @@
       double Py = mcp.getPY();
       double Pz = mcp.getPZ();
 
-      double Pt2 = (Px * Px) + (Py * Py);
-      double Pt = Math.sqrt(Pt2);
-      double P = Math.sqrt(Pt2 + (Pz * Pz));
+      double P = Math.sqrt((Px * Px) + (Py * Py) + (Pz * Pz));
       double Phi = Math.atan2(Py, Px);
       if (Phi < 0)
       {
          Phi += (2 * Math.PI);
       }
 
-      double CosTh = Pz / P;
-      double Theta = Math.acos(CosTh);
+      double Theta = Math.acos(Pz / P);
 
       // Simulate position smearing on a sphere of radius 2 meters
       radius = 2000.0;
@@ -82,19 +79,23 @@
       double y = (radius * Py) / P;
       double z = (radius * Pz) / P;
 
-      double drPhi = transDist * CosTh;
-      double transPhi = Math.PI * rand.nextDouble();
-      x = x + (drPhi * Math.sin(transPhi));
-      y = y + (drPhi * Math.cos(transPhi));
-      z = z + (transDist * Math.sin(Theta));
+      //these vectors vt and vs (orthonorm.) span a plane perpendicular to the momentum vector,
+      //so smearing with a transdist will involve a lin. comb. of these
+      double[] vt = {-Math.cos(Theta)*Math.cos(Phi),-Math.cos(Theta)*Math.sin(Phi),Math.sin(Theta)};
+      double[] vs = {Math.sin(Phi),-Math.cos(Phi),0};
+      
+      //restricted to [0,PI] since transdist can be negative
+      double alpha = rand.nextDouble() * Math.PI;
+      x = x + transDist * ( Math.cos(alpha) * vt[0] + Math.sin(alpha) * vs[0] );
+      y = y + transDist * ( Math.cos(alpha) * vt[1] + Math.sin(alpha) * vs[1] );
+      z = z + transDist * ( Math.cos(alpha) * vt[2] + Math.sin(alpha) * vs[2] );
 
       phi = Math.atan2(y, x);
       if (phi < 0)
       {
          phi += (2 * Math.PI);
       }
-      cosTh = z / radius;
-      theta = Math.acos(cosTh);
+      theta = Math.acos(z / radius);
    }
    
    public Particle getMCParticle()

lcsim/src/org/lcsim/mc/fast/reconstructedparticle
MCFastReconstructedParticleDriver.java 1.7 -> 1.8
diff -u -r1.7 -r1.8
--- MCFastReconstructedParticleDriver.java	18 Jul 2005 19:18:26 -0000	1.7
+++ MCFastReconstructedParticleDriver.java	9 Aug 2005 18:34:45 -0000	1.8
@@ -14,6 +14,7 @@
 import java.util.ArrayList;
 import java.util.List;
 import java.util.Random;
+import org.lcsim.event.MCParticle;
 import org.lcsim.event.Cluster;
 import org.lcsim.event.EventHeader;
 import org.lcsim.event.ReconstructedParticle;
@@ -55,6 +56,8 @@
     protected void process(EventHeader event)
     {
         
+        boolean hist = getHistogramLevel() > 0;
+        
         if (IDEff == null)
         {
            ConditionsSet conditions = getConditionsManager().getConditions("IDEfficiency");
@@ -80,23 +83,37 @@
                 // electrons and muons are special
                 if((abs(pdgid)== 11) && (rand.nextDouble() < IDEff.getElectronEff()))
                 {
-                    type = new Electron((int)Math.signum(pdgid));
+                    type = new Electron(rt.getCharge());
                 }
                 else if((abs(pdgid)== 13) && (rand.nextDouble() < IDEff.getMuonEff()))
                 {
-                    type = new Muon((int)Math.signum(pdgid));
+                    type = new Muon(rt.getCharge());
                 }
                 else
                 {
                     type = new Pion(rt.getCharge());
                 }
                 
-                aida.histogram1D("track-particle", 150, -10, 10).fill( (t.getPX()*t.getPX() + t.getPY()*t.getPY() + t.getPZ()*t.getPZ() + type.mass()*type.mass()) - p.getEnergy());
+                if ((p.getEnergy() > 10) && (hist)) 
+                {
+                    if (Math.abs(t.getTrackParameter(4)) < 1)
+                    {
+                        aida.histogram1D("track-particle", 150, -0.0003, 0.0003).fill( ( Math.sqrt(t.getPX()*t.getPX() + t.getPY()*t.getPY() + t.getPZ()*t.getPZ()) - Math.sqrt(p.getPX()*p.getPX() + p.getPY()*p.getPY() + p.getPZ()*p.getPZ()) ) / ( p.getPX()*p.getPX() + p.getPY()*p.getPY() + p.getPZ()*p.getPZ()) );
+                    }  
+                    if ((Math.abs(t.getTrackParameter(4)) < 2.16) && (Math.abs(t.getTrackParameter(4)) > 1.96) && (p.getMomentum().magnitude() > 80) && (p.getMomentum().magnitude() < 120))
+                    {
+                        aida.histogram1D("track-particle-cut", 150, -0.01, 0.01).fill( ( Math.sqrt(t.getPX()*t.getPX() + t.getPY()*t.getPY() + t.getPZ()*t.getPZ()) - Math.sqrt(p.getPX()*p.getPX() + p.getPY()*p.getPY() + p.getPZ()*p.getPZ()) ) / (Math.sqrt( p.getPX()*p.getPX() + p.getPY()*p.getPY() + p.getPZ()*p.getPZ())) );
+                    }
+                }
                 
                 // assume pion for remaining charged tracks
             MCFastReconstructedParticle rp = new MCFastReconstructedParticle(t, type, p);
             rpList.add(rp);
-            aida.histogram1D("recon-particle", 150, -10, 10).fill(rp.getEnergy()-p.getEnergy());
+            
+                if (hist)
+                {
+                    aida.histogram1D("recon-particle", 150, -5, 5).fill((rp.getEnergy()-p.getEnergy())/(Math.sqrt(p.getEnergy())));
+                }
             }
        }
         
@@ -104,15 +121,18 @@
         List<Cluster> clusters = event.getClusters();
         for(Cluster c : clusters)
         {
-            ParticleType type = null;
             Particle p = null;
+            ParticleType type = null;
             // photons for EM
             if( c instanceof ReconEMCluster)
             {
                 ReconEMCluster emc = (ReconEMCluster) c;
                 p = emc.getMCParticle();
                 type = new Photon();
-                aida.histogram1D("cluster-particle", 150, -10, 10).fill(emc.getEnergy()-emc.getMCParticle().getEnergy());
+                if (hist)
+                {
+                    aida.histogram1D("photonCLS-particle", 150, -3, 3).fill((emc.getEnergy()-emc.getMCParticle().getEnergy())/(Math.sqrt(emc.getMCParticle().getEnergy())));
+                }
             }
             // assume a KZeroLong here for had cluster
             else if(c instanceof ReconHADCluster)
@@ -120,7 +140,10 @@
                 ReconHADCluster emc = (ReconHADCluster) c;
                 p = emc.getMCParticle();
                 int pdgid = p.getPDGID();
-                aida.histogram1D("cluster-particle", 150, -10, 10).fill(emc.getEnergy()-emc.getMCParticle().getEnergy());
+                if (hist)
+                {
+                    aida.histogram1D("hadronCLS-particle", 150, -3, 3).fill((emc.getEnergy()-emc.getMCParticle().getEnergy())/(Math.sqrt(emc.getMCParticle().getEnergy())));
+                }
                 
                 if ((abs(pdgid)==2112) && (rand.nextDouble() < IDEff.getNeutronEff()))
                 {
@@ -133,7 +156,11 @@
             }
             MCFastReconstructedParticle rp = new MCFastReconstructedParticle(c, type, p);
             rpList.add(rp);
-            aida.histogram1D("recon-particle", 150, -10, 10).fill(rp.getEnergy()-p.getEnergy());
+            if (hist)
+            {
+                aida.histogram1D("recon-particle", 150, -10, 10).fill(rp.getEnergy()-p.getEnergy());
+            }
+            
         }
         // add the reconstructedparticles to the event
         event.put(event.RECONSTRUCTEDPARTICLES, rpList, ReconstructedParticle.class, 0);

lcsim/src/org/lcsim/mc/fast/tracking
MCFastTracking.java 1.5 -> 1.6
diff -u -r1.5 -r1.6
--- MCFastTracking.java	18 Jul 2005 19:18:26 -0000	1.5
+++ MCFastTracking.java	9 Aug 2005 18:34:45 -0000	1.6
@@ -20,7 +20,9 @@
 public class MCFastTracking extends Driver implements ConditionsListener
 {
    private TrackResolutionTables parm;
+   private SimpleTables SmTbl;
    private boolean beamSpotConstraint;
+   private boolean simple;
    private final static double[] IP = { 0, 0, 0 }; 
 
    public MCFastTracking()
@@ -32,6 +34,12 @@
    {
       this.beamSpotConstraint = beamSpotConstraint;
    }
+   
+   public MCFastTracking(boolean beamSpotConstraint, boolean simple)
+   {
+       this.beamSpotConstraint = beamSpotConstraint;
+       this.simple = simple;
+   }
 
    public void setBeamSpotConstraint(boolean beamSpotConstraint)
    {
@@ -68,6 +76,13 @@
          parm = setTrackResolutionTables(conditions, beamSpotConstraint);
       }
       
+      if (SmTbl == null)
+      {
+         ConditionsSet conditions = getConditionsManager().getConditions("SimpleTrack");
+         conditions.addConditionsListener(this);
+         SmTbl = new SimpleTables(conditions);
+      }
+      
       double bField = event.getDetector().getFieldMap().getField(IP)[2];
       boolean hist = getHistogramLevel() > 0;
 
@@ -103,7 +118,7 @@
             continue;
          }
 
-         trackList.add(new ReconTrack(bField, parm, getRandom(), p, hist));
+         trackList.add(new ReconTrack(bField, parm, SmTbl, getRandom(), p, hist, simple));
       }
       event.put(EventHeader.TRACKS, trackList, Track.class, 0);
    }
@@ -111,6 +126,8 @@
    public void conditionsChanged(ConditionsEvent event)
    {
       ConditionsSet conditions = getConditionsManager().getConditions("TrackParameters");
+      ConditionsSet simpleconditions = getConditionsManager().getConditions("SimpleTrack");
       parm = setTrackResolutionTables(conditions, beamSpotConstraint);
+      SmTbl = new SimpleTables(simpleconditions);
    }
 }

lcsim/src/org/lcsim/mc/fast/tracking
ReconTrack.java 1.5 -> 1.6
diff -u -r1.5 -r1.6
--- ReconTrack.java	2 Jul 2005 06:40:45 -0000	1.5
+++ ReconTrack.java	9 Aug 2005 18:34:45 -0000	1.6
@@ -9,6 +9,7 @@
 import java.util.Collections;
 import java.util.List;
 import java.util.Random;
+import org.lcsim.mc.fast.tracking.SimpleTables;
 
 
 /**
@@ -17,7 +18,7 @@
  * are provided. <br>
  *
  * @author Tony Johnson, Wolfgang Walkowiak
- * @version $Id: ReconTrack.java,v 1.5 2005/07/02 06:40:45 ngraf Exp $
+ * @version $Id: ReconTrack.java,v 1.6 2005/08/09 18:34:45 ngraf Exp $
  */
 public class ReconTrack implements Track
 {
@@ -40,20 +41,21 @@
    transient private Particle mc;
    private int m_tcharge;
 
-   ReconTrack(double bField, TrackResolutionTables parm, Random rand, Particle mc, boolean hist)
+   ReconTrack(double bField, TrackResolutionTables parm, SimpleTables SmTbl, Random rand, Particle mc, boolean hist, boolean simple)
    {
       this.mc = mc;
 
       // get original momentum from MCParticle
       // convert to helical parameters
       m_nosmear = new DocaTrackParameters(bField, mc);
+      
+      double pt = m_nosmear.getPt();
 
       if (hist)
       {
 
          double r  = Math.abs(m_nosmear.getD0());
-         double pt = m_nosmear.getPt();
-
+         
          AIDA aida = AIDA.defaultInstance();
          aida.cloud1D("ptsqr").fill(pt*pt);
          aida.cloud1D("pt").fill(pt);
@@ -72,8 +74,15 @@
       m_nosmear.fillErrorMatrix(getErrorMatrixFromTable(table, abscth, ptot));
 
       // smear tracks according to error matrix
-      m_smear = (DocaTrackParameters) SmearTrack.smearTrack(bField, m_nosmear, rand);
-
+      if (simple == true)
+      {
+          m_smear = (DocaTrackParameters) SmearTrackSimple.smearTrackSimple(bField, m_nosmear, rand, SmTbl, pt, hist);
+          //double[] slice = {0, 0, 1, 0, 0};
+          //m_smear = (DocaTrackParameters) SmearTrackSimpleII.SmearTrackSimpleII(bField, m_nosmear, rand, SmTbl, slice, hist);
+      } else {
+          m_smear = (DocaTrackParameters) SmearTrack.smearTrack(bField, m_nosmear, rand);
+      }
+      
       if (hist)
       {
          AIDA aida = AIDA.defaultInstance();

lcsim/src/org/lcsim/mc/fast/tracking
SmearTrack.java 1.1 -> 1.2
diff -u -r1.1 -r1.2
--- SmearTrack.java	1 Feb 2005 19:42:48 -0000	1.1
+++ SmearTrack.java	9 Aug 2005 18:34:45 -0000	1.2
@@ -18,6 +18,7 @@
     */
    static DocaTrackParameters smearTrack(double bField, TrackParameters noSmear, Random rand)
    {
+
       // get error matrix and do a sanity check
       Matrix M = Matrix.constructWithCopy(noSmear.getErrorMatrix());
       if (M.det() <= 0.)

lcsim/src/org/lcsim/mc/fast
MCFast.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- MCFast.java	1 Jul 2005 23:51:57 -0000	1.2
+++ MCFast.java	9 Aug 2005 18:34:46 -0000	1.3
@@ -12,10 +12,15 @@
 public class MCFast extends Driver
 {
    /** Creates a new instance of MCFast */
-   public MCFast()
+   public MCFast(boolean beamSpotConstraint, boolean simple)
    {
-      add(new MCFastTracking());
+      add(new MCFastTracking(beamSpotConstraint, simple));
       add(new MCFastRonan());
       add(new MCFastReconstructedParticleDriver());
    }
+   
+   public MCFast()
+   {
+       this(false,false);
+   }
 }
CVSspam 0.2.8