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

diff -N SRPdijetm.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ SRPdijetm.java	12 Aug 2009 18:33:59 -0000	1.1
@@ -0,0 +1,133 @@

+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 SRPdijetm extends Driver
+{
+    private String _jlinkmap;
+    private String _scinjets;
+    private AIDA aida = AIDA.defaultInstance();
+
+    public SRPdijetm()
+    {
+
+    }
+    
+    protected void process(EventHeader event)
+    {
+        super.process(event);
+
+        // array of ReconstructedParticles to add to the event...
+        List<ReconstructedParticle> rpList = new ArrayList<ReconstructedParticle>();
+
+        double sEtot = 0;
+        double sjpx = 0;
+        double sjpy = 0;
+        double sjpz = 0;
+        double ujm2 = 0.;
+        //  get scintillator jets for uncorrected mass
+        List<ReconstructedParticle> scjets = event.get(ReconstructedParticle.class,_scinjets);
+        for (ReconstructedParticle scjet : scjets)
+        {
+            sEtot += scjet.getEnergy();
+            Hep3Vector jp = scjet.getMomentum();
+            sjpx += jp.x();
+            sjpy += jp.y();
+            sjpz += jp.z();
+            ujm2 += scjet.getMass()*scjet.getMass();
+        }
+        double sptot = Math.sqrt(sjpx*sjpx+sjpy*sjpy+sjpz*sjpz);
+        double sjmass = Math.sqrt(sEtot*sEtot-sptot*sptot+ujm2);
+        aida.cloud1D("Uncorrected Scin Jet Mass").fill(sjmass);
+        aida.cloud1D("Uncorrected Scin Jet ETot").fill(sEtot);
+
+        //  now using linked jets, do Cerenkov correction
+        Map<ReconstructedParticle, ReconstructedParticle> SCjmap = (Map<ReconstructedParticle, ReconstructedParticle>) event.get(_jlinkmap);
+        int nlinks = SCjmap.size();
+        double jpx = 0.;
+        double jpy = 0.;
+        double jpz = 0.;
+        double Etot = 0.;
+        double mjtot2 = 0.;
+        int jind = 0;
+        double[] ejet = new double[10];
+        double[] pzjet = new double[10];
+
+        //  loop over scintillator jets and correct with linked cerenkov jets
+        if (nlinks == 2)
+        {
+//        for(ReconstructedParticle jet : SCjmap.keySet())
+        for(ReconstructedParticle jet : scjets)
+        {
+            ReconstructedParticle cjet = SCjmap.get(jet);
+            double scrat = cjet.getEnergy()/jet.getEnergy();  // C/S ratio for jet
+            if (scrat>1.) scrat = 1.0;
+            aida.cloud1D("C over S ratio per jet").fill(scrat);
+//           System.out.println("CS ratio " + scrat);
+//            double sjetCE = jet.getEnergy()/(.555-.200*scrat+.643*scrat*scrat);
+            double sjetCE = jet.getEnergy()/(.480+.628*scrat-1.085*scrat*scrat+.975*scrat*scrat*scrat);
+//            double sjetCE = jet.getEnergy()/(.68+.31*scrat);
+            aida.cloud1D("Correction Factor for Jets").fill(jet.getEnergy()/sjetCE);
+            Hep3Vector jmom = jet.getMomentum();
+            aida.cloud1D("Uncorrected Jet pT").fill(Math.sqrt(jmom.x()*jmom.x()+jmom.y()*jmom.y()));
+            double ucjm2 = jet.getMass()*jet.getMass();
+//            System.out.println("Uncorrected Jet Mass " + jet.getMass());
+            double jpmag = jmom.magnitude();
+            double cjpmag = jmom.magnitude()*sjetCE/jet.getEnergy();
+            double jm2 = sjetCE*sjetCE-cjpmag*cjpmag;
+//            System.out.println("Corrected Jet Mass " + Math.sqrt(jm2));
+//            System.out.println("jet mass squared " + jm2);
+//            System.out.println("Jet 4V  " + sjetCE + " " +  jmom.x()*sjetCE/jet.getEnergy()  + " " +  jmom.y()*sjetCE/jet.getEnergy()  + " " +  jmom.z()*sjetCE/jet.getEnergy());
+            double cjpx = jmom.x()*sjetCE/jet.getEnergy();
+            double cjpy = jmom.y()*sjetCE/jet.getEnergy();
+            aida.cloud1D("Corrected Jet pT").fill(Math.sqrt(cjpx*cjpx+cjpy*cjpy));
+            jpx += jmom.x()*sjetCE/jet.getEnergy();
+            jpy += jmom.y()*sjetCE/jet.getEnergy();
+            jpz += jmom.z()*sjetCE/jet.getEnergy();
+            Etot += sjetCE;
+            mjtot2 += jm2;
+            ejet[jind] = sjetCE;
+            pzjet[jind] = jmom.z()*sjetCE/jet.getEnergy();
+            jind++;
+        }
+        double emin = 0.;
+        double pzj = 0.;
+        if (jind<3)
+        {
+            emin = ejet[1];
+            pzj = pzjet[1];
+            if (ejet[0] < ejet[1])
+            {
+                emin = ejet[0];
+                pzj = pzjet[0];
+            }
+//            aida.cloud1D("ycut for 2 jets").fill(2.*emin*(emin-Math.abs(pzj))/(Etot*Etot));
+        }
+        double ptot = Math.sqrt(jpx*jpx+jpy*jpy+jpz*jpz);
+        double jinv = Math.sqrt(Etot*Etot-ptot*ptot+mjtot2);
+        aida.cloud1D("Ceren Corrected Scin Jet Mass").fill(jinv);
+        aida.cloud1D("Ceren Corrected Scin Jet ETot").fill(Etot);
+//        System.out.println("Total S Jet Mass " + jinv);
+        }
+    }
+
+  public void setInputScChJetMap(String jlinkmap)
+  {
+      _jlinkmap = jlinkmap;
+  }
+
+    public void setInputScinJets(String scinjets)
+  {
+      _scinjets = scinjets;
+  }
+  
+}