lcsim-contrib/src/main/java/org/lcsim/contrib/SteveMagill

diff -N CorrRPdijetm.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ CorrRPdijetm.java	14 Feb 2012 19:33:37 -0000	1.1
@@ -0,0 +1,105 @@

+package org.lcsim.contrib.SteveMagill;
+
+import java.util.*;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.aida.AIDA;
+
+/*
+ * PFASelect includes fiducial cut to select region for PFA analysis
+ */
+
+public class CorrRPdijetm extends Driver
+{
+    private String _scinjets;
+    private String _jdlinkmap;
+    private String _stjetmap;
+    private AIDA aida = AIDA.defaultInstance();
+
+    public CorrRPdijetm()
+    {
+
+    }
+    
+    protected void process(EventHeader event)
+    {
+        super.process(event);
+
+        double sEtot = 0;
+        double sjpx = 0;
+        double sjpy = 0;
+        double sjpz = 0;
+        double ujm2 = 0.;
+        //  get E M of qqbar
+        List<Double> qqm = event.get(Double.class, "DJmass");
+        List<Double> qqE = event.get(Double.class, "DJEnergy");
+        //  get corrected jets
+        List<ReconstructedParticle> scjets = event.get(ReconstructedParticle.class,_scinjets);
+        //  get scintillator jet-track jet link map
+        Map<ReconstructedParticle, ReconstructedParticle> jtmap = (Map<ReconstructedParticle, ReconstructedParticle>) event.get(_stjetmap);
+        //  get jet-deltam map
+        Map<ReconstructedParticle, Double> jdmap = (Map<ReconstructedParticle, Double>) event.get(_jdlinkmap);
+        if (scjets.size() == 2)
+        {
+        for (ReconstructedParticle scjet : scjets)
+        {
+            double jdm = 0.;
+            sEtot += scjet.getEnergy();
+            double jm = scjet.getMass();
+            aida.cloud1D("CS Clus Jet Mass before dM Corr").fill(jm);
+            if (jtmap.containsKey(scjet))
+            {
+                ReconstructedParticle jtl = jtmap.get(scjet);
+                if (jdmap.containsKey(jtl)) jdm = jdmap.get(jtl);
+            }
+            double jetm = jm-jdm;
+            if (jetm<0.) jetm = 0.;
+            aida.cloud1D("dM for CS Clus Jet").fill(jdm);
+            aida.cloud1D("CS Clus Jet Mass after dM Corr").fill(jetm);
+            double corrpmag = Math.sqrt(scjet.getEnergy()*scjet.getEnergy()-jetm*jetm);
+            double pmag = scjet.getMomentum().magnitude();
+            double cp = corrpmag/pmag;
+            Hep3Vector jp = scjet.getMomentum();
+            sjpx += jp.x()*cp;
+            sjpy += jp.y()*cp;
+            sjpz += jp.z()*cp;
+//            ujm2 += scjet.getMass()*scjet.getMass();
+            ujm2 += jetm*jetm;
+        }
+        double sptot = Math.sqrt(sjpx*sjpx+sjpy*sjpy+sjpz*sjpz);
+        double sjmass = Math.sqrt(sEtot*sEtot-sptot*sptot+ujm2);
+        aida.cloud1D("CS Clus DiJet Mass dM").fill(sjmass);
+        aida.cloud1D("CS Clus DiJet ETot dM").fill(sEtot);
+        double jm = 0.;
+        for (Double qm : qqm)
+        {
+            jm = qm.doubleValue();
+        }
+        double jE = 0.;
+        for (Double qE : qqE)
+        {
+            jE = qE.doubleValue();
+        }
+        aida.cloud1D("Diff CS Clus DiJets Mass dM").fill(sjmass-jm);
+        aida.cloud1D("Diff CS Clus DiJets Energy dM").fill(sEtot-jE);
+        }
+    }
+
+    public void setInputScinJets(String scinjets)
+  {
+      _scinjets = scinjets;
+  }
+
+    public void setSTJetMap(String stjetmap)
+  {
+      _stjetmap = stjetmap;
+  }
+
+    public void setTJetDMMap(String jdlinkmap)
+  {
+      _jdlinkmap = jdlinkmap;
+  }
+  
+}