 |
 |
| Commit in mcd-analysis/src/main/java/org/lcsim/mcd/analysis on MAIN | |||
| HiggsAnalysisDriver.java | +416 | added 1.1 | |
Simple exploratory driver for looking at Z/gamma background and Higgs samples at 125.0 GeV.
diff -N HiggsAnalysisDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ HiggsAnalysisDriver.java 19 Mar 2013 19:58:49 -0000 1.1 @@ -0,0 +1,416 @@
+package org.lcsim.mcd.analysis;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.mc.fast.MCFast;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.event.util.JetDriver;
+import hep.physics.jet.JetFinder;
+import hep.physics.jet.DurhamJetFinder;
+import hep.physics.jet.EventShape;
+import hep.physics.vec.BasicHepLorentzVector;
+import hep.physics.vec.HepLorentzVector;
+import org.lcsim.contrib.Cassell.recon.Cheat.MakePerfectReconParticles;
+/**
+ *
+ * @author aconway
+ */
+public class HiggsAnalysisDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ MCFast MCFast;
+ JetDriver jDrive;
+ JetFinder jFind;
+ double _ctcut = 0.9;
+ double _jetcut = 0.2;
+
+ public HiggsAnalysisDriver() {
+ MCFast = new MCFast();
+ }
+
+ public void process(EventHeader event) {
+ this.add(MCFast);
+ this.processChildren(event);
+ this.remove(MCFast);
+
+ MCParticle boson = getBoson(event);
+ String MCTag = getMCTag(boson);
+ double bosonInvMass = Math.sqrt(boson.getEnergy()*boson.getEnergy()
+ - boson.getMomentum().magnitudeSquared());
+ aida.cloud1D("boson masses").fill(boson.getMass());
+ aida.cloud1D("boson energies").fill(boson.getEnergy());
+ aida.cloud1D("boson invariant masses")
+ .fill(bosonInvMass);
+
+
+ List<MCParticle> FSChildren = event.get(MCParticle.class,"GenFinalStateParticles");
+
+ List<MCParticle> ctCutList = cosThetaCut(FSChildren, _ctcut);
+ List<MCParticle> nuCutList = neutrinoCut(FSChildren);
+ List<MCParticle> ctnuCutList = neutrinoCut(ctCutList);
+
+ event.put("noCut", FSChildren, MCParticle.class,0);
+ event.put("ctCut", ctCutList, MCParticle.class, 0);
+ event.put("nuCut", nuCutList, MCParticle.class, 0);
+ event.put("ctnuCut", ctnuCutList, MCParticle.class, 0);
+
+ MakePerfectReconParticles mprp1 = new MakePerfectReconParticles("noCut","rcp noCut");
+ MakePerfectReconParticles mprp2 = new MakePerfectReconParticles("ctCut","rcp ctCut");
+ MakePerfectReconParticles mprp3 = new MakePerfectReconParticles("nuCut","rcp nuCut");
+ MakePerfectReconParticles mprp4 = new MakePerfectReconParticles("ctnuCut","rcp ctnuCut");
+
+ this.add(mprp1);
+ this.add(mprp2);
+ this.add(mprp3);
+ this.add(mprp4);
+ this.processChildren(event);
+ this.remove(mprp1);
+ this.remove(mprp2);
+ this.remove(mprp3);
+ this.remove(mprp4);
+
+ String evt;
+ if (boson.getMass() >= 124.9) {
+ evt = "Z*/";
+ }
+// else if (boson.getMass() >= 110) {
+// evt = "High Z*/";
+// }
+ else if (boson.getMass() >= 70) {
+ evt = "Z+gamma/";
+ } else {
+ evt = "Low Z*/";
+ }
+ makePlots(event, FSChildren, "noCut", MCTag, evt);
+ makePlots(event, ctCutList, "ctCut", MCTag, evt);
+ makePlots(event, nuCutList, "nuCut", MCTag, evt);
+ makePlots(event, ctnuCutList, "ctnuCut", MCTag, evt);
+ }
+
+ private void makePlots(EventHeader event, List<MCParticle> mcpl, String prefix, String decay, String evt) {
+
+ makeJetPlots(event,prefix,"/ALL",evt);
+ makeThrustPlots(mcpl,evt+prefix+"/ALL");
+ makeGammaPlots(mcpl,evt+prefix+"/ALL");
+ makeBosonPlots(mcpl,evt+prefix+"/ALL", event);
+
+ makeJetPlots(event,prefix,decay,evt);
+ makeThrustPlots(mcpl,evt+prefix+"/"+decay);
+ makeGammaPlots(mcpl,evt+prefix+"/"+decay);
+ if(decay.equals("WW")) makeWWPlots(mcpl,evt+prefix+"/"+decay, event);
+ }
+
+ private void makeWWPlots(List<MCParticle> mcpl, String prefix, EventHeader event) {
+ List<MCParticle> Ws = new ArrayList<MCParticle>();
+ MCParticle boson = getBoson(event);
+ List<MCParticle> FS = cosThetaCut(neutrinoCut(
+ event.get(MCParticle.class,"GenFinalStateParticles")),_ctcut);
+ for (MCParticle p:boson.getDaughters()) {
+ if(Math.abs(p.getPDGID()) == 24) {
+ Ws.add(p);
+ }
+ }
+ String decay = getWWDecay(Ws);
+ prefix = prefix + "/" + decay + "/";
+ aida.cloud1D(prefix+"Num W's").fill(Ws.size());
+ List<MCParticle> dslist = new ArrayList<MCParticle>();
+ for(MCParticle p:Ws) {
+ List<MCParticle> ds = p.getDaughters();
+ for (MCParticle d:ds) {
+ if (d.getGeneratorStatus() == 3) {
+ dslist.add(d);
+ aida.cloud1D(prefix+"W decay PDGIDs").fill(d.getPDGID());
+ }
+ }
+ Hep3Vector pt = p.getMomentum();
+ aida.cloud1D(prefix+"W Masses").fill(p.getMass());
+ aida.cloud1D(prefix+"W Momenta").fill(pt.magnitude());
+ aida.cloud1D(prefix+"W Cos Theta").fill(VecOp.cosTheta(pt));
+ Hep3Vector pt_t = new BasicHep3Vector(pt.x(),pt.y(),0.0);
+ aida.cloud1D(prefix+"W Transverse Momenta").fill(pt_t.magnitude());
+ }
+ aida.cloud1D(prefix+"num FSP's").fill(FS.size());
+ makeThrustPlots(FS,prefix);
+ }
+
+ private void makeBosonPlots(List<MCParticle> mcpl, String prefix, EventHeader event) {
+ MCParticle boson = getBoson(event);
+ Hep3Vector p = boson.getMomentum();
+ prefix = prefix+"/boson/";
+ aida.cloud1D(prefix+"boson mass").fill(boson.getMass());
+ aida.cloud1D(prefix+"boson energy").fill(boson.getEnergy());
+ aida.cloud1D(prefix+"boson momentum").fill(p.magnitude());
+ double invMass = Math.sqrt(boson.getEnergy()*boson.getEnergy() - p.magnitudeSquared());
+ aida.cloud1D(prefix+"boson invariant mass").fill(invMass);
+ aida.cloud1D(prefix+"boson cos theta").fill(VecOp.cosTheta(p));
+ Hep3Vector p_t = new BasicHep3Vector(p.x(),p.y(),0);
+ aida.cloud1D(prefix+"boson transverse momentum").fill(p_t.magnitude());
+
+ for (MCParticle mcp:mcpl) {
+ if( (mcp.getPDGID()==22) && (mcp.getEnergy()>20.0)) {
+ Hep3Vector pg = mcp.getMomentum();
+ double imass = Math.sqrt(
+ (mcp.getEnergy()+boson.getEnergy())*(mcp.getEnergy()+boson.getEnergy())
+ - VecOp.add(p, pg).magnitudeSquared());
+ aida.cloud1D(prefix+"boson+gamma inv mass").fill(imass);
+ aida.cloud1D(prefix+"boson+gamma total energy").fill(mcp.getEnergy()+boson.getEnergy());
+ aida.cloud2D(prefix+"boson+gamma energy correl.").fill(boson.getEnergy(), mcp.getEnergy());
+ aida.cloud2D(prefix+"boson+gamma momentum correl.").fill(p.magnitude(),pg.magnitude());
+ aida.cloud2D(prefix+"boson+gamma cos theta correl.").fill(VecOp.cosTheta(p), VecOp.cosTheta(pg));
+ aida.cloud2D(prefix+"boson+gamma mass-energy correl.").fill(boson.getMass(),mcp.getEnergy());
+ aida.cloud1D(prefix+"boson inv mass ct cut").fill(invMass);
+ }
+ }
+ aida.histogram1D(prefix+"boson mass hist 1",50,115,125).fill(boson.getMass());
+ aida.histogram1D(prefix+"boson mass hist 2",50,120,125).fill(boson.getMass());
+ aida.histogram1D(prefix+"boson mass hist 3",50,124,125).fill(boson.getMass());
+ aida.histogram1D(prefix+"boson mass hist 4",50,124.9,125).fill(boson.getMass());
+ }
+
+ private void makeGammaPlots(List<MCParticle> mcpl, String prefix) {
+ for (MCParticle p:mcpl) {
+ if ( (p.getPDGID()==22) && (p.getEnergy()>20.0)) {
+ aida.cloud2D(prefix+"/gamma/costheta vs en"
+ ).fill(p.getEnergy(),VecOp.cosTheta(p.getMomentum()));
+ aida.cloud1D(prefix+"/gamma/costheta").fill(VecOp.cosTheta(p.getMomentum()));
+ aida.cloud1D(prefix+"/gamma/energy").fill(p.getEnergy());
+ Hep3Vector p_t = new BasicHep3Vector(p.getMomentum().x(),p.getMomentum().y(),0);
+ aida.cloud1D(prefix+"/gamma/transverse momentum").fill(p_t.magnitude());
+ }
+ }
+ }
+
+ private void makeThrustPlots(List<MCParticle> mcpl, String prefix) {
+ EventShape es = getEventShape(mcpl);
+ Hep3Vector thrust = es.thrust();
+ Hep3Vector majAxis = es.majorAxis();
+ Hep3Vector minAxis = es.minorAxis();
+ double ob = es.oblateness();
+
+ if(thrust.x()>0) aida.cloud1D(prefix+"/Thrust").fill(thrust.x());
+ aida.cloud1D(prefix+"/Major Axis").fill(thrust.y());
+ aida.cloud1D(prefix+"/Minor Axis").fill(thrust.z());
+ aida.cloud1D(prefix+"/Oblateness").fill(ob);
+ }
+
+ private void makeJetPlots(EventHeader event, String prefix, String decay,String evt) {
+ jFind = new DurhamJetFinder(_jetcut);
+ jDrive = new JetDriver();
+ jDrive.setFinder(jFind);
+ jDrive.setInputCollectionName("rcp "+prefix);
+ jDrive.setOutputCollectionName("Jets "+prefix);
+ this.add(jDrive);
+ this.processChildren(event);
+ this.remove(jDrive);
+ List<ReconstructedParticle> jets = event.get(ReconstructedParticle.class, "Jets "+prefix);
+ List<MCParticle> allFSPs = event.get(MCParticle.class,prefix);
+
+ int nJets = jFind.njets();
+ aida.cloud1D(evt+prefix+"/"+decay+"/num Jets").fill(nJets);
+ prefix = evt+prefix+"/"+decay+"/"+nJets+"jets";
+
+
+ int ps_in_jets = 0;
+ for (int i=0; i<jFind.njets(); i++) {
+ aida.cloud1D(prefix+"/paricles per jet").fill(jFind.nParticlesPerJet(i));
+ ps_in_jets += jFind.nParticlesPerJet(i);
+ }
+ aida.cloud1D(prefix+"/fewest tracks in a jet").fill(jFind.fewestTracks());
+ aida.cloud1D(prefix+"/MCPs not in jets").fill(allFSPs.size()-ps_in_jets);
+
+
+ double sumE = 0.0;
+ double sumM = 0.0;
+ double sumP = 0.0;
+ Hep3Vector sumPv = new BasicHep3Vector(0.0,0.0,0.0);
+ for(ReconstructedParticle jet:jets) {
+ Hep3Vector pj = jet.getMomentum();
+ aida.cloud1D(prefix+"/jet energy").fill(jet.getEnergy());
+ aida.cloud1D(prefix+"/jet mass").fill(jet.getMass());
+ aida.cloud1D(prefix+"/jet momentum").fill(pj.magnitude());
+ aida.cloud1D(prefix+"/jet cosTheta").fill(VecOp.cosTheta(pj));
+ aida.cloud1D(prefix+"/jet charge").fill(jet.getCharge());
+ double invMass = Math.sqrt(jet.getEnergy()*jet.getEnergy()
+ - pj.magnitudeSquared());
+ aida.cloud1D(prefix+"/jet invariant mass").fill(invMass);
+ Hep3Vector pj_t = new BasicHep3Vector(pj.x(),pj.y(),0);
+ aida.cloud1D(prefix+"/jet transverse momentum").fill(pj_t.magnitude());
+ sumE += jet.getEnergy();
+ sumM += jet.getMass();
+ sumP += jet.getMomentum().magnitude();
+ VecOp.add(sumPv,jet.getMomentum());
+ }
+ double invMass = Math.sqrt(sumE*sumE - sumPv.magnitudeSquared());
+ aida.cloud1D(prefix+"/total jet invariant mass").fill(invMass);
+ aida.cloud1D(prefix+"/total jet energy").fill(sumE);
+ aida.cloud1D(prefix+"/total jet mass").fill(sumM);
+ aida.cloud1D(prefix+"/total jet momentum").fill(sumP);
+ }
+
+ private String getWWDecay(List<MCParticle> Ws) {
+ String rv = "Unknown";
+ boolean jjjj = true;
+ boolean llnn = true;
+ boolean jj = true;
+ boolean ll = true;
+ List<MCParticle> dslist = new ArrayList<MCParticle>();
+ for (MCParticle p:Ws) {
+ List<MCParticle> ds = p.getDaughters();
+ for(MCParticle d:ds) {
+ if(d.getGeneratorStatus() == 3) {
+ dslist.add(d);
+ int pdgid = Math.abs(d.getPDGID());
+ if (pdgid <= 5) llnn = false;
+ else if(pdgid < 15) jjjj = false;
+ }
+ }
+ }
+ if((!jjjj) && (!llnn)) rv = "jjlnu";
+ else if((jjjj) && (!llnn)) rv = "jjjj";
+ else if((llnn) && (!jjjj)) rv = "ll";
+ return rv;
+ }
+
+ private EventShape getEventShape(List<MCParticle> mcpl) {
+ List<HepLorentzVector> hlv = new ArrayList<HepLorentzVector>();
+ for(MCParticle p:mcpl) {
+ hlv.add(new BasicHepLorentzVector(p.getEnergy(),p.getMomentum()));
+ }
+ EventShape rv = new EventShape();
+ rv.setEvent(hlv);
+ return rv;
+ }
+
+ private List<MCParticle> cosThetaCut(List<MCParticle> mcpl, double ctcut) {
+ List<MCParticle> rv = new ArrayList<MCParticle>();
+ for (MCParticle p:mcpl) {
+ double cosT = Math.abs(VecOp.cosTheta(p.getMomentum()));
+ if (cosT < ctcut) {
+ rv.add(p);
+ }
+ }
+ return rv;
+ }
+
+ private List<MCParticle> neutrinoCut(List<MCParticle> mcpl) {
+ List<MCParticle> rv = new ArrayList<MCParticle>();
+ for(MCParticle p:mcpl) {
+ if( !( (p.getPDGID() == 12)
+ || (p.getPDGID() == 14)
+ || (p.getPDGID() == 16)
+ || (p.getPDGID() == 18) ) )
+ rv.add(p);
+ }
+ return rv;
+ }
+
+ private MCParticle getBoson(EventHeader event) {
+ MCParticle boson = event.getMCParticles().get(0);
+ for (MCParticle p:event.getMCParticles()) {
+ if ( ((p.getPDGID() == 25) || (p.getPDGID() == 23))
+ && (p.getGeneratorStatus() == 3)
+ && (Math.abs(p.getParents().get(0).getPDGID()) == 13)) {
+ //&& (p.getMass() > 40.0) ){
+ boson = p;
+ }
+ }
+ return boson;
+ }
+
+ private String getMCTag(MCParticle boson) {
+ List<MCParticle> children = boson.getDaughters();
+
+ String chan = "Unknown";
+ if (!children.isEmpty()) {
+ int pdgid1 = 0;
+ int pdgid2 = 0;
+ boolean ffound = false;
+ boolean sfound = false;
+ for (MCParticle p:children) {
+ 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 = "UDS";
+ break;
+ case 2 :
+ chan = "UDS";
+ break;
+ case 3 :
+ chan = "UDS";
+ break;
+ case 4 :
+ chan = "CCBAR";
+ break;
+ case 5 :
+ chan = "BBBAR";
+ break;
+ case 6 :
+ chan = "TTBAR";
+ break;
+ case 15 :
+ chan = "TAUTAU";
+ break;
+ case 21 :
+ chan = "GLUGLU";
+ break;
+ case 22 :
+ chan = "GAMGAM";
+ break;
+ case 23 :
+ chan = "ZZ";
+ break;
+ case 24 :
+ chan = "WW";
+ break;
+ default :
+ chan = "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;
+ }
+
+ // recursively searches the parent's children for final state particles
+ public List<MCParticle> getAllFSChildren(MCParticle parent) {
+ List<MCParticle> AllChildren = new ArrayList<MCParticle>();
+
+ List<MCParticle> children = parent.getDaughters();
+
+ for (MCParticle child:children) {
+ if (child.getGeneratorStatus() == 1) {
+ AllChildren.add(child);
+ }
+ else {
+ AllChildren.addAll(getAllFSChildren(child));
+ }
+ }
+ return AllChildren;
+ }
+}
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