mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N MCWTaggingDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MCWTaggingDriver.java 8 Feb 2013 16:44:39 -0000 1.1
@@ -0,0 +1,200 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.mcd.analysis;
+
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.event.util.JetDriver;
+import org.lcsim.mc.fast.MCFast;
+import org.lcsim.util.Driver;
+import org.lcsim.event.LCRelation;
+import org.lcsim.digisim.MyLCRelation;
+
+/**
+ *
+ * @author aconway
+ *
+ */
+
+enum DecayChannel {
+ UDS, BBAR, CCBAR, TTBAR, // quarks
+ GLUGLU, GAMGAM, WW, ZZ, ZGAM, // bosons
+ TAUTAUBAR, // leptons
+ OTHER, UNKNOWN // other
+}
+
+public class MCWTaggingDriver extends Driver {
+
+ public MCWTaggingDriver() {
+ add(new MCFast());
+ add(new JetDriver());
+ }
+
+
+
+ // Array of W decay products by PDGID <-> index
+ int[] W_DECAY_COUNTS = new int[20];
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ public void process(EventHeader event) {
+ this.processChildren(event);
+ MCParticle Higgs = findHiggs(event);
+ DecayChannel dChan = getMCTag(Higgs);
+
+ if (dChan == DecayChannel.WW) {
+ List<MCParticle> HiggsChildren = new ArrayList<MCParticle>(Higgs.getDaughters());
+ for (MCParticle p:Higgs.getDaughters()) {
+ if (Math.abs(p.getPDGID()) != 24) {
+ HiggsChildren.remove(p);
+ } else if (p.getPDGID() == 24) {
+
+ }
+ }
+ MCParticle Wp = HiggsChildren.get(0);
+ MCParticle Wm = HiggsChildren.get(1);
+
+ List<ReconstructedParticle> jets = event.get(ReconstructedParticle.class,"Jets");
+
+ List<ReconstructedParticle> Wpjets = dijetRPMatch(jets,Wp);
+ List<ReconstructedParticle> Wmjets = dijetRPMatch(jets,Wm);
+
+ List<LCRelation> lcrlist = new ArrayList<LCRelation>();
+ lcrlist.add(new MyLCRelation((List<ReconstructedParticle>) Wpjets, (MCParticle) Wp));
+ lcrlist.add(new MyLCRelation((List<ReconstructedParticle>) Wmjets, (MCParticle) Wm));
+ event.put("DijetToMCP", lcrlist, ReconstructedParticle.class,0);
+ }
+ }
+
+ // Given a list of jets and a particle, find the two jets that best reconstruct the particle
+ public List<ReconstructedParticle> dijetRPMatch(List<ReconstructedParticle> jets, MCParticle p){
+ List<ReconstructedParticle> rvlist = new ArrayList<ReconstructedParticle>();
+ rvlist.add(jets.get(0));
+ rvlist.add(jets.get(1));
+ double best_mass_diff = 10000000.0;
+ double dijet_inv_mass = 0;
+
+ for (ReconstructedParticle jet1:jets) {
+ for (ReconstructedParticle jet2:jets) {
+ double p1Charge = p.getCharge();
+ double j1Charge = jet1.getCharge();
+ double j2Charge = jet2.getCharge();
+ if (j1Charge + j2Charge != p1Charge) {
+ break;
+ }
+ double mj1 = jet1.getMass();
+ double mj2 = jet2.getMass();
+ double ej1 = jet1.getEnergy();
+ double ej2 = jet2.getEnergy();
+ Hep3Vector pj1 = jet1.getMomentum();
+ Hep3Vector pj2 = jet2.getMomentum();
+ dijet_inv_mass = Math.sqrt(
+ mj1*mj1 + mj2*mj2 + 2.0*(ej1*ej2 - VecOp.dot(pj1, pj2)));
+
+ double diff_dj_mass = Math.abs(dijet_inv_mass - p.getMass());
+
+ if (diff_dj_mass < best_mass_diff) {
+ rvlist.clear();
+ rvlist.add(jet1);
+ rvlist.add(jet2);
+ best_mass_diff = diff_dj_mass;
+ }
+ }
+ }
+ aida.cloud1D("W-matched dijet masses").fill(dijet_inv_mass);
+ aida.cloud1D("W-matched dijet mass minus MC W mass").fill(best_mass_diff);
+
+ return rvlist;
+ }
+
+ public MCParticle findHiggs(EventHeader event) {
+ MCParticle Higgs = event.getMCParticles().get(0);
+ for (MCParticle p:event.getMCParticles()) {
+ if ( (p.getPDGID() == 25) && (p.getGeneratorStatus() == 3) ){
+ Higgs = p;
+ }
+ }
+ return Higgs;
+ }
+
+ // Determine the type of higgs decay by identifying its children by PDGID
+ // Return as DecayChannel enum
+ public DecayChannel getMCTag(MCParticle higgs) {
+ List<MCParticle> HiggsChildren = higgs.getDaughters();
+
+ DecayChannel chan = DecayChannel.UNKNOWN;
+ if (!HiggsChildren.isEmpty()) {
+ int pdgid1 = 0;
+ int pdgid2 = 0;
+ boolean ffound = false;
+ boolean sfound = false;
+ for (MCParticle p:HiggsChildren) {
+ if (p.getGeneratorStatus()==3) { // Documentation state
+ int pdgid = p.getPDGID();
+
+ if ( !ffound && !sfound ) {
+ pdgid1 = pdgid;
+ ffound = true;
+ }
+ if ( ffound && !sfound ) {
+ pdgid2 = pdgid;
+ sfound = true;
+ }
+ }
+ }
+ if (ffound && sfound) {
+ if (Math.abs(pdgid1)==Math.abs(pdgid2)) {
+ aida.cloud1D("Decay Channels").fill(Math.abs(pdgid1));
+ switch (Math.abs(pdgid1)) {
+ case 1 :
+ chan = DecayChannel.UDS;
+ break;
+ case 2 :
+ chan = DecayChannel.UDS;
+ break;
+ case 3 :
+ chan = DecayChannel.UDS;
+ break;
+ case 4 :
+ chan = DecayChannel.CCBAR;
+ break;
+ case 5 :
+ chan = DecayChannel.BBAR;
+ break;
+ case 6 :
+ chan = DecayChannel.TTBAR;
+ break;
+ case 21 :
+ chan = DecayChannel.GLUGLU;
+ break;
+ case 22 :
+ chan = DecayChannel.GAMGAM;
+ break;
+ case 23 :
+ chan = DecayChannel.ZZ;
+ break;
+ case 24 :
+ chan = DecayChannel.WW;
+ break;
+ default :
+ chan = DecayChannel.OTHER;
+ break;
+ }
+ }
+ else {
+ aida.cloud1D("Unknown decay channels by child PDGID").fill(pdgid1);
+ aida.cloud1D("Unknown decay channels by child PDGID").fill(pdgid2);
+ }
+ }
+ }
+ return chan;
+ }
+}