Commit in mcd-analysis/src/main/java/org/lcsim/mcd/analysis on MAIN | |||
MCWTaggingDriver.java | +200 | added 1.1 |
Haven't tested it yet. Simple driver for looking at MC information in higgs -> WW events.
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; + } +}
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