lcsim-contrib/src/main/java/org/lcsim/contrib/SteveMagill
diff -N MCqqbarME.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MCqqbarME.java 12 Aug 2009 18:29:28 -0000 1.1
@@ -0,0 +1,125 @@
+package org.lcsim.contrib.SteveMagill;
+
+import java.util.List;
+import java.util.ArrayList;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import org.lcsim.event.util.JetDriver;
+import org.lcsim.util.Driver;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.jet.*;
+
+/*
+ * PFASelect includes fiducial cut to select region for PFA analysis
+ */
+
+public class MCqqbarME extends Driver
+{
+ private String _mcname;
+ private String _rpname;
+ private AIDA aida = AIDA.defaultInstance();
+
+ public MCqqbarME()
+ {
+
+ }
+
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+
+ // get true qqbar E and M and get Nayeli's mass correction
+ List<MCParticle> mcl = event.get(MCParticle.class,"MCParticle");
+ double qqbarM = 0.;
+ double qqbarE = 0.;
+ double mcpx = 0.;
+ double mcpy = 0.;
+ double mcpz = 0.;
+ double Etot = 0.;
+ double pxprime = 0.;
+ double pyprime = 0.;
+ double pzprime = 0.;
+ for(MCParticle p:mcl)
+ {
+ int gs = p.getGeneratorStatus();
+ int pid = Math.abs(p.getPDGID());
+ double pch = p.getCharge();
+// if (gs==1 || (gs==0 && pid==2112)) // only want neutrons which are called FS, so don't use this
+ if (gs==1 && (pid != 12 && pid != 14 && pid != 16)) // all FS particles except neutrinos
+ {
+ Etot += p.getEnergy();
+ mcpx += p.getPX();
+ mcpy += p.getPY();
+ mcpz += p.getPZ();
+ Hep3Vector rv = p.getEndPoint();
+ double rmag = rv.magnitude();
+ pxprime += p.getEnergy()*rv.x()/rmag;
+ pyprime += p.getEnergy()*rv.y()/rmag;
+ pzprime += p.getEnergy()*rv.z()/rmag;
+ }
+ int id = Math.abs(p.getPDGID());
+ if( (id == 1 )||(id == 2 )||(id == 3 ) )
+ {
+ int parid = Math.abs(p.getParents().get(0).getPDGID());
+ if(parid == 24 || parid == 23)
+ {
+ qqbarM = p.getParents().get(0).getMass();
+ qqbarE = p.getParents().get(0).getEnergy();
+// System.out.println("qqbar Mass " + qqbarM);
+ }
+ }
+ }
+ double pptot = Math.sqrt(pxprime*pxprime+pyprime*pyprime+pzprime*pzprime);
+ double mprime = Math.sqrt(Etot*Etot-pptot*pptot);
+ double ptot = Math.sqrt(mcpx*mcpx+mcpy*mcpy+mcpz*mcpz);
+ double minv = Math.sqrt(Etot*Etot-ptot*ptot);
+ double delm = mprime-minv;
+ aida.cloud1D("MC FS Invariant Mass").fill(minv);
+ aida.cloud1D("MC FS Prime Invariant Mass").fill(mprime);
+ aida.cloud1D("Delta M Mp M").fill(delm);
+// System.out.println("MC FS Invariant Mass " + minv);
+
+ aida.cloud1D("qqbar Parent Mass").fill(qqbarM);
+ aida.cloud1D("qqbar Parent Energy").fill(qqbarE);
+
+// make reconstructed particles out of MC FS list
+ double epho = 0.;
+ double ETot = 0.;
+ List<ReconstructedParticle> rpList = new ArrayList<ReconstructedParticle>();
+ List<MCParticle> mcfsp = event.get(MCParticle.class,_mcname);
+ for (MCParticle mcf : mcfsp)
+ {
+ int mcid = Math.abs(mcf.getPDGID());
+// System.out.println("PDGID " + mcid);
+ if (mcid != 12 && mcid != 14 && mcid != 16)
+ {
+ double mcE = mcf.getEnergy();
+ ETot += mcE;
+ if (mcid == 22) epho += mcE;
+ Hep3Vector mcP = mcf.getMomentum();
+ BaseReconstructedParticle rp = new BaseReconstructedParticle(mcE, mcP);
+ rp.setMass(mcf.getMass());
+ rp.setCharge(mcf.getCharge());
+ rpList.add(rp);
+ }
+ }
+ event.put(_rpname,rpList, ReconstructedParticle.class,0);
+ double phorat = epho/ETot;
+ aida.cloud1D("Fraction of photons").fill(phorat);
+
+ }
+
+ public void setMCParticleName(String mcname)
+ {
+ _mcname = mcname;
+ }
+
+ public void setMCRecoParticleName(String rpname)
+ {
+ _rpname = rpname;
+ }
+
+}