mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N MCHiggsDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MCHiggsDriver.java 8 Feb 2013 16:43:36 -0000 1.1
@@ -0,0 +1,133 @@
+package org.lcsim.mcd.analysis;
+
+import java.util.List;
+import java.util.ArrayList;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+import org.lcsim.mc.fast.MCFast;
+
+/**
+ *
+ * @author agias
+ */
+public class MCHiggsDriver extends Driver{
+
+ public MCHiggsDriver()
+ {
+ add(new MCFast());
+ }
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ protected String mcParticlesSkimmedName = "MCParticle";
+
+ public void process(EventHeader event) {
+ List<MCParticle> particles = event.get(MCParticle.class,mcParticlesSkimmedName);
+ List<MCParticle> children = new ArrayList<MCParticle>();
+
+ boolean hfound = false;
+ for (MCParticle particle : particles) {
+ // 25 == h0/H01
+ if ((particle.getPDGID()==25)&&(!hfound)) {
+ children.addAll(particle.getDaughters());
+ hfound=true;
+ }
+ }
+ MCParticle p1 = children.get(0);
+ MCParticle p2 = children.get(0);
+
+ for (MCParticle particle : children) {
+ /*
+ * Only count positive PDGID's because negatives
+ * are just antiparticles
+ * 4 == c
+ * 5 == b
+ * 15 == tau-
+ * 21 == gluon
+ * 22 == gamma
+ * 23 == Zo
+ * 24 == W
+ * 25 == h0/H01
+ */
+ if ((particle.getPDGID() > 0) && (particle.getPDGID() != 25)) {
+ aida.cloud1D("decay modes").fill(particle.getPDGID());
+ p1 = particle;
+ }
+ if ((particle.getPDGID() < 0)) {
+ p2 = particle;
+ }
+ }
+
+ if (p1.getPDGID() == 5) {
+ Hep3Vector p = p1.getMomentum();
+ Hep3Vector p_bar = p2.getMomentum();
+
+ double mass1 = p1.getMass();
+ double en1 = p1.getEnergy();
+ double mass2 = p2.getMass();
+ double en2 = p2.getEnergy();
+
+ Hep3Vector z_unit = new BasicHep3Vector(0,0,1.);
+
+ double p1_m = p.magnitude();
+ double p2_m = p_bar.magnitude();
+
+ double pt = VecOp.dot(p, z_unit);
+ double pt_bar = VecOp.dot(p_bar, z_unit);
+
+ aida.cloud1D("b momentum").fill(p1_m);
+ aida.cloud1D("b_bar momentum").fill(p2_m);
+ aida.histogram1D("b momentum hist",100,62.3108,62.3126).fill(p1_m);
+ aida.histogram1D("b_bar momentum hist",100,62.3,62.4).fill(p2_m);
+
+ aida.cloud1D("b transverse momentum").fill(pt);
+ aida.cloud1D("b_bar transverse momentum").fill(pt_bar);
+ aida.histogram1D("b transverse momentum hist", 100,0,63).fill(pt);
+ aida.histogram1D("b_bar transverse momentum hist", 100,0,63).fill(pt_bar);
+
+
+ double inv_mass = Math.sqrt(mass1*mass1 + mass2*mass2 + 2.*(en1*en2-VecOp.dot(p, p_bar)));
+ aida.cloud1D("H invariant mass").fill(inv_mass);
+ aida.histogram1D("H inv mass hist",100,124.990,124.995).fill(inv_mass);
+ aida.histogram1D("H inv mass hist big",200,100,130).fill(inv_mass);
+
+ aida.cloud1D("b phi").fill(VecOp.phi(p));
+ aida.cloud1D("bbar phi").fill(VecOp.phi(p_bar));
+ aida.cloud1D("b cos theta").fill(VecOp.cosTheta(p));
+ aida.cloud1D("bbar cos theta").fill(VecOp.cosTheta(p_bar));
+ }
+ if (p1.getPDGID() == 22) {
+ double mass1 = p1.getMass();
+ double ener1 = p1.getEnergy();
+ double mass2 = p2.getMass();
+ double ener2 = p2.getEnergy();
+
+ Hep3Vector pt1 = p1.getMomentum();
+ Hep3Vector pt2 = p2.getMomentum();
+ double p1_m = pt1.magnitude();
+ double p2_m = pt2.magnitude();
+
+ Hep3Vector z_unit = new BasicHep3Vector(0,0,1.);
+
+ double p1trans = VecOp.dot(pt1, z_unit);
+ double p2trans = VecOp.dot(pt2,z_unit);
+
+
+ double inv_mass = Math.sqrt(2*(ener1*ener2 - VecOp.dot(pt1, pt2)));
+
+ aida.cloud1D("gamma1 momentum").fill(p1_m);
+ aida.cloud1D("gamma2 momentum").fill(p2_m);
+ aida.cloud1D("gamma1 transverse momentum").fill(p1trans);
+ aida.cloud1D("gamma2 transverse momentum").fill(p2trans);
+ aida.cloud1D("gamma higgs inv mass").fill(inv_mass);
+ }
+ aida.cloud1D("num children").fill(children.size());
+ }
+
+}