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;
+ }
+
+}