Commit in mcd-analysis/src/main/java/org/lcsim/mcd/analysis on MAIN | |||
MCHiggsDriver.java | +133 | added 1.1 |
Changed name to end with Driver for clarity. Haven't tested yet. Just a basic driver for looking at MC information in an stdhep mumu -> higgs file.
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()); + } + +}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1