mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -u -r1.9 -r1.10
--- MCHiggs.java 19 Dec 2012 18:29:04 -0000 1.9
+++ MCHiggs.java 8 Feb 2013 16:39:52 -0000 1.10
@@ -10,6 +10,7 @@
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.VecOp;
+import org.lcsim.mc.fast.MCFast;
/**
*
@@ -19,6 +20,7 @@
public MCHiggs()
{
+ add(new MCFast());
}
private AIDA aida = AIDA.defaultInstance();
@@ -37,6 +39,8 @@
hfound=true;
}
}
+ MCParticle p1 = children.get(0);
+ MCParticle p2 = children.get(0);
for (MCParticle particle : children) {
/*
@@ -53,58 +57,76 @@
*/
if ((particle.getPDGID() > 0) && (particle.getPDGID() != 25)) {
aida.cloud1D("decay modes").fill(particle.getPDGID());
+ p1 = particle;
}
- if (particle.getPDGID() == 5) {
- Hep3Vector p = particle.getMomentum();
- Hep3Vector p_bar = VecOp.neg(p);
-
- double mass = particle.getMass();
- double en = particle.getEnergy();
-
- Hep3Vector z_unit = new BasicHep3Vector(0,0,1.);
-
- double p_m = p.magnitude();
-
- double pt = VecOp.dot(p, z_unit);
- double pt_bar = VecOp.dot(p_bar, z_unit);
-
- //aida.cloud1D("b momentum").fill(p_m);
- aida.histogram1D("b momentum hist",100,62.3108,62.3126).fill(p_m);
- //aida.cloud1D("b transverse momentum").fill(pt);
- aida.histogram1D("b transverse momentum hist", 100,0,63).fill(pt);
- //aida.cloud1D("b_bar transverse momentum").fill(pt_bar);
- aida.histogram1D("b_bar transverse momentum hist", 100,0,63).fill(pt_bar);
-
- // b and b_bar are produced with equal energy, opposite momenta
- double inv_mass = Math.sqrt(2.*mass*mass + 2.*(en*en-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("b cos theta").fill(VecOp.cosTheta(p));
- }
- if (particle.getPDGID() == 22) {
- double mass = particle.getMass();
- double ener = particle.getEnergy();
- Hep3Vector p = particle.getMomentum();
- Hep3Vector p2 = VecOp.neg(p);
- double p_m = p.magnitude();
- double p2_m = p_m;
-
- Hep3Vector z_unit = new BasicHep3Vector(0,0,1.);
-
- double pt = VecOp.dot(p, z_unit);
- double pt2 = VecOp.dot(p2,z_unit);
-
-
- double inv_mass = Math.sqrt(2*(ener*ener - VecOp.dot(p, p2)));
-
- aida.cloud1D("gamma momentum").fill(p_m);
- aida.cloud1D("gamma transverse momentum").fill(pt);
- aida.cloud1D("gamma higgs inv mass").fill(inv_mass);
+ 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());
}