 |
 |
| Commit in mcd-analysis/src/main/java/org/lcsim/mcd/analysis on MAIN | |||
| MCbtaggingDriver.java | +559 | added 1.1 | |
Simple driver for looking at MC information in higgs -> bb events. Still has some problems and is probably hard to read.
diff -N MCbtaggingDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MCbtaggingDriver.java 8 Feb 2013 16:45:43 -0000 1.1 @@ -0,0 +1,559 @@
+//package org.lcsim.mcd.analysis;
+
+import java.util.Arrays;
+import java.util.List;
+import java.util.ArrayList;
+import java.util.Set;
+import org.lcsim.spacegeom.SpaceVector;
+import org.lcsim.spacegeom.CartesianVector;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.Track;
+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;
+import org.lcsim.event.TrackState;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+
+import org.lcsim.recon.vertexing.zvtop4.ZvTop;
+import org.lcsim.recon.vertexing.zvtop4.ZvFitter;
+import org.lcsim.recon.vertexing.zvtop4.ZvVertex;
+import org.lcsim.recon.vertexing.zvtop4.ZvTrack;
+
+/**
+ *
+ * @author aconway
+ */
+
+enum DecayChannel {
+ UDS, BBAR, CCBAR, TTBAR, // quarks
+ GLUGLU, GAMGAM, WW, ZZ, ZGAM, // bosons
+ TAUTAUBAR, // leptons
+ OTHER, UNKNOWN // other
+}
+
+public class MCbtaggingDriver extends Driver {
+
+ public MCbtaggingDriver(){
+ /* adds the following collections from MC data:
+ * Clusters
+ * ClustersToMCP
+ * GenFinalStateParticles
+ * MCFastReconstructedParticles
+ * Tracks (of class ReconTrack)
+ * TracksToMCP (LCRelation)
+ */
+ add(new MCFast());
+ }
+
+ private AIDA aida = AIDA.defaultInstance();
+ /*
+ * Bottom Meson PDGID's
+ * B0 511
+ * B+ 521
+ * B_s^0 531
+ * B_c^+ 541
+ */
+ int[] BMesonPDGIDs = { -511,511,
+ -521,521,
+ -531,531,
+ -541,541};
+
+ int[] BBaryonPDGIDs = {-5122,5122, //Lambda_b^0
+ -5112,5112, //Sigmas, Chis and Omegas
+ -5212,5212,
+ -5222,5222,
+ -5114,5114,
+ -5214,5214,
+ -5224,5224,
+ -5132,5132,
+ -5232,5232,
+ -5312,5312,
+ -5322,5322,
+ -5314,5314,
+ -5324,5324,
+ -5332,5332,
+ -5334,5334,
+ -5142,5142,
+ -5242,5242,
+ -5412,5412,
+ -5422,5422,
+ -5414,5414,
+ -5424,5424,
+ -5342,5342,
+ -5432,5432,
+ -5434,5434,
+ -5442,5442,
+ -5444,5444,
+ -5512,5512,
+ -5522,5522,
+ -5514,5514,
+ -5524,5524,
+ -5532,5532,
+ -5534,5534,
+ -5542,5542,
+ -5544,5544,
+ -5554,5554};
+
+
+
+
+ public void process(EventHeader event) {
+ this.processChildren(event);
+ MCParticle Higgs = findHiggs(event);
+ DecayChannel dChan = getMCTag(Higgs);
+
+ List<LCRelation> TracksToMCP = event.get(LCRelation.class,"TracksToMCP");
+ List<Track> tracks = event.get(Track.class,"Tracks");
+
+ /*
+ * DOCUMENTATION = 3
+ * INTERMEDIATE = 2
+ * FINAL_STATE = 1
+ */
+
+
+ if (dChan == DecayChannel.BBAR) {
+ List<MCParticle> HiggsChildren = new ArrayList<MCParticle>(Higgs.getDaughters());
+ for (MCParticle p:Higgs.getDaughters()) {
+ if (Math.abs(p.getPDGID()) != 5) {
+ HiggsChildren.remove(p);
+ }
+ }
+ // This might not be safe for other datasets...
+ MCParticle b = HiggsChildren.get(0);
+ MCParticle bbar = HiggsChildren.get(1);
+ // Check it here:
+ aida.cloud1D("MC charge b").fill(b.getCharge());
+ aida.cloud1D("MC charge bbar").fill(bbar.getCharge());
+
+ //
+ // Plot calculations/data from MC b's
+ //
+ Hep3Vector p_b = b.getMomentum();
+ double p_bx = p_b.x();
+ double p_by = p_b.y();
+ double p_bz = p_b.z();
+ Hep3Vector p_bbar = bbar.getMomentum();
+ double p_bbarx = p_bbar.x();
+ double p_bbary = p_bbar.y();
+ double p_bbarz = p_bbar.z();
+
+ Hep3Vector p_b_trans = new BasicHep3Vector(p_bx,p_by,0.0);
+ Hep3Vector p_bbar_trans = new BasicHep3Vector(p_bbarx,p_bbary,0.0);
+
+ aida.cloud1D("MC momentum b (total)").fill(p_b.magnitude());
+ aida.cloud1D("MC momentum bbar (total)").fill(p_bbar.magnitude());
+ aida.cloud1D("MC momentum b (trans)").fill(p_b_trans.magnitude());
+ aida.cloud1D("MC momentum bbar (trans)").fill(p_bbar_trans.magnitude());
+
+ double mass_b = b.getMass();
+ double ener_b = b.getEnergy();
+ double mass_bbar = bbar.getMass();
+ double ener_bbar = bbar.getEnergy();
+
+ // M_T^2 = m^2 + |p_t|^2
+ double mass_b_trans =
+ Math.sqrt(mass_b*mass_b + p_b_trans.magnitudeSquared());
+ double mass_bbar_trans =
+ Math.sqrt(mass_bbar*mass_bbar + p_bbar_trans.magnitudeSquared());
+
+ /*
+ // E_T = E*sin(theta)
+ double ener_b_trans =
+ ener_b*Math.sqrt(p_bx*p_bx+p_by*p_by)/p_b.magnitude();
+ double ener_bbar_trans =
+ ener_bbar*Math.sqrt(p_bbarx*p_bbarx+p_bbary*p_bbary)/p_bbar.magnitude();
+ */
+ // E_T = sqrt( |p_t|^2 + m^2 )
+ double ener_b_trans =
+ Math.sqrt( p_b_trans.magnitudeSquared() + mass_b*mass_b);
+ double ener_bbar_trans =
+ Math.sqrt ( p_bbar_trans.magnitudeSquared() + mass_bbar*mass_bbar);
+
+ aida.cloud1D("MC mass b").fill(mass_b);
+ aida.cloud1D("MC mass bbar").fill(mass_bbar);
+ aida.cloud1D("MC energy b").fill(ener_b);
+ aida.cloud1D("MC energy bbar").fill(ener_bbar);
+ aida.cloud1D("MC mass b (trans)").fill(mass_b_trans);
+ aida.cloud1D("MC mass bbar (trans)").fill(mass_bbar_trans);
+ aida.cloud1D("MC energy b (trans)").fill(ener_b_trans);
+ aida.cloud1D("MC energy bbar (trans)").fill(ener_bbar_trans);
+
+ double phi_b = VecOp.phi(p_b);
+ double phi_bbar = VecOp.phi(p_bbar);
+ aida.cloud1D("MC b phi").fill(phi_b);
+ aida.cloud1D("MC bba phi").fill(phi_bbar);
+
+ double cosTheta_b = VecOp.cosTheta(p_b);
+ double cosTheta_bbar = VecOp.cosTheta(p_bbar);
+ aida.cloud1D("MC cosTheta b").fill(cosTheta_b);
+ aida.cloud1D("MC cosTheta bbar").fill(cosTheta_bbar);
+
+ // M^2 = m_1^2 + m_2^2 + 2*(E_1*E_2 - p_1 \dot p_2)
+ double inv_mass = Math.sqrt(
+ mass_b*mass_b + mass_bbar*mass_bbar
+ + 2.*(ener_b*ener_bbar-VecOp.dot(p_b, p_bbar)));
+ double inv_mass_trans = Math.sqrt(
+ mass_b*mass_b + mass_bbar*mass_bbar
+ + 2.*(ener_b_trans*ener_bbar_trans-VecOp.dot(p_b_trans, p_bbar_trans)));
+
+ aida.cloud1D("MC Higgs invariant mass").fill(inv_mass);
+ aida.cloud1D("MC Higgs invariant mass (trans)").fill(inv_mass_trans);
+
+ bDecayPlots(b, "b");
+ bDecayPlots(bbar,"bbar");
+ boolean bmeson = false;
+ boolean bbaryon = false;
+
+ List<MCParticle> b_FSChildren = getAllFSChildren(b);
+ List<MCParticle> bbar_FSChildren = getAllFSChildren(bbar);
+
+ bTrackPlots(event,b,"b");
+
+
+ //
+ //
+ //
+ /*
+
+ List<Track> b_jet_tracks = new ArrayList<Track>();
+ List<Track> bbar_jet_tracks = new ArrayList<Track>();
+
+ for (LCRelation lcr:TracksToMCP) {
+ if (b_FSChildren.contains((MCParticle) lcr.getTo())) {
+ b_jet_tracks.add((Track) lcr.getFrom());
+ }
+ else if (bbar_FSChildren.contains((MCParticle) lcr.getTo())) {
+ bbar_jet_tracks.add((Track) lcr.getFrom());
+ }
+ }
+
+
+ ZvFitter myFitter = new ZvFitter(event.getDetector());
+ SpaceVector jetVec = new CartesianVector(p_bx,p_by,p_bz);
+ ZvTop myZvTop = new ZvTop(jetVec,myFitter);
+ List<ZvVertex> b_vertices = myZvTop.findVertices(b_jet_tracks);
+ event.put("bVertices", b_vertices, ZvVertex.class, 0);
+ aida.cloud1D("b jet num vertices").fill(b_vertices.size());
+
+ jetVec = new CartesianVector(p_bbarx,p_bbary,p_bbarz);
+ myZvTop = new ZvTop(jetVec,myFitter);
+ List<ZvVertex> bbar_vertices = myZvTop.findVertices(b_jet_tracks);
+ event.put("bbarVertices", bbar_vertices, ZvVertex.class, 0);
+ aida.cloud1D("bbar jet tracks multiplicity").fill(bbar_vertices.size());
+ */
+ /*
+ List<TrackState> bTracks = new ArrayList<TrackState>();
+ for (ZvVertex vtx:vertices) {
+ for (ZvTrack zvtrk:vtx.getVtxTracks()) {
+ bTracks.addAll(zvtrk.getTrackStates());
+ }
+ }
+ event.put("bTracks", bTracks, TrackState.class,0);
+ */
+
+ }
+
+ if (dChan == DecayChannel.CCBAR) {
+ List<MCParticle> HiggsChildren = new ArrayList<MCParticle>(Higgs.getDaughters());
+ for (MCParticle p:Higgs.getDaughters()) {
+ if (Math.abs(p.getPDGID()) != 4) {
+ HiggsChildren.remove(p);
+ }
+ }
+ // This might not be safe for other datasets...
+ MCParticle c = HiggsChildren.get(0);
+ MCParticle cbar = HiggsChildren.get(1);
+
+ bTrackPlots(event,c,"c");
+ }
+ /*
+ List<MCParticle> particles = event.get(MCParticle.class,"MCParticle");
+
+ for (MCParticle particle:particles) {
+ boolean bmeson = false;
+ boolean bbaryon = false;
+
+ for (int i=0; i<BMesonPDGIDs.length; i++) {
+ if (particle.getPDGID() == BMesonPDGIDs[i]) {
+ bmeson = true;
+ aida.cloud1D("all bmesons").fill(particle.getPDGID());
+ }
+ }
+ if (!bmeson) {
+ for (int i=0; i<BBaryonPDGIDs.length; i++) {
+ if (particle.getPDGID() == BBaryonPDGIDs[i]) {
+ bbaryon = true;
+ aida.cloud1D("all bbaryons").fill(particle.getPDGID());
+ }
+ }
+ }
+ if (bmeson || bbaryon) {
+ List<MCParticle> BFSChildren = getAllFSChildren(particle);
+ List<LCRelation> ChildTracks = new ArrayList<LCRelation>();
+
+ for (LCRelation rel:TracksToMCP) {
+ if (BFSChildren.contains((MCParticle) rel.getTo())) {
+ ChildTracks.add(rel);
+ }
+ }
+ aida.cloud1D("B decay num tracks").fill(ChildTracks.size());
+
+ Hep3Vector BFinalPos = particle.getEndPoint();
+ double transverse_len = Math.sqrt(
+ BFinalPos.x()*BFinalPos.x()+BFinalPos.y()*BFinalPos.y());
+ Hep3Vector BMomentum = particle.getMomentum();
+ double transverse_mom = Math.sqrt(
+ BMomentum.x()*BMomentum.x()+BMomentum.y()*BMomentum.y());
+ aida.cloud1D("B decay length (total)").fill(BFinalPos.magnitude());
+ aida.cloud1D("B decay length (transverse)").fill(transverse_len);
+ aida.cloud1D("B momenta").fill(BMomentum.magnitude());
+ aida.cloud1D("B transverse momenta").fill(transverse_mom);
+
+ double sumE = 0.0;
+ double sumPx = 0.0;
+ double sumPy = 0.0;
+ double sumPz = 0.0;
+
+ for (LCRelation rel:ChildTracks) {
+ ReconTrack trk = (ReconTrack) rel.getFrom();
+ MCParticle prt = (MCParticle) rel.getTo();
+
+ aida.cloud1D("B child D0's (XY)").fill(trk.getTrackParameter(0));
+ aida.cloud1D("B child z impact parameter").fill(trk.getTrackParameter(3));
+ sumE += prt.getEnergy();
+ sumPx += prt.getMomentum().x();
+ sumPy += prt.getMomentum().y();
+ sumPz += prt.getMomentum().z();
+ }
+ double realBinvMass = particle.getEnergy()*particle.getEnergy()
+ - BMomentum.magnitude()*BMomentum.magnitude();
+ double realBinvMassTrans = particle.getEnergy()*particle.getEnergy()
+ - transverse_mom*transverse_mom;
+
+ double sumP = Math.sqrt(sumPx*sumPx+sumPy*sumPy+sumPz*sumPz);
+ double sumPt = Math.sqrt(sumPx*sumPx+sumPy*sumPy);
+
+ double fakeBinvMass = Math.sqrt(sumE*sumE - sumP*sumP);
+ double fakeBinvMassTrans = Math.sqrt(sumE*sumE - sumPt*sumPt);
+
+ aida.cloud1D("B inv mass fraction - tracks").fill(fakeBinvMass/realBinvMass);
+ aida.cloud1D("B inv mass fraction - trans tracks").fill(fakeBinvMassTrans/realBinvMassTrans);
+ }
+ }
+ */
+ }
+ // recursively searches the parent's children for all particles
+ // disregards intermediates with no children
+ public List<MCParticle> getAllChildren(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 if (!child.getDaughters().isEmpty()) {
+ AllChildren.add(child);
+ AllChildren.addAll(getAllChildren(child));
+ }
+ }
+ return AllChildren;
+ }
+
+ // 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;
+ }
+
+ public void bDecayPlots(MCParticle quark, String quark_name) {
+ List<MCParticle> AllChildren = getAllChildren(quark);
+
+ for (MCParticle p:AllChildren) {
+ for (int i=0; i<BMesonPDGIDs.length; i++) {
+ if ( (p.getPDGID() == BMesonPDGIDs[i]) && (p.getOrigin().equals(new BasicHep3Vector(0.0,0.0,0.0))) ) {
+ aida.cloud1D(quark_name +" decay meson PDGIDs").fill(Math.abs(p.getPDGID()));
+ Hep3Vector l = p.getEndPoint();
+ Hep3Vector l_t = new BasicHep3Vector(l.x(),l.y(),0.0);
+ double t = p.getDaughters().get(0).getProductionTime();
+ double m = p.getMass();
+ Hep3Vector mom = p.getMomentum();
+ Hep3Vector mom_t = new BasicHep3Vector(mom.x(),mom.y(),0.0);
+
+ double prop_l = l.magnitude()*m/mom.magnitude();
+ double prop_lt = l_t.magnitude()*m/mom_t.magnitude();
+ double prop_t = 2.99792458*t; // c cm/ns TODO: check this...
+
+ aida.cloud1D("B Meson Mass").fill(m);
+
+ aida.cloud1D("B Meson Decay Length (total)").fill(l.magnitude());
+ aida.cloud1D("B Meson Decay Length (trans)").fill(l_t.magnitude());
+ aida.cloud1D("B Meson Decay Time").fill(t);
+ aida.cloud1D("B Meson Proper Length (total)").fill(prop_l);
+ aida.cloud1D("B Meson Proper Length (trans)").fill(prop_lt);
+ aida.cloud1D("B ctau").fill(prop_t);
+
+ aida.histogram1D("HIST B Meson Proper Decay Length (total)",100,0.0,2.5).fill(prop_l);
+ aida.histogram1D("HIST B Meson Proper Decay Length (trans)",100,0.0,2.5).fill(prop_lt);
+ aida.histogram1D("HIST B Meson Decay Time",100,0.0,25).fill(t);
+
+ }
+ }
+ for (int i=0; i<BBaryonPDGIDs.length; i++) {
+ if (p.getPDGID() == BBaryonPDGIDs[i]) {
+ aida.cloud1D(quark_name +" decay baryon PDGIDs").fill(p.getPDGID());
+ }
+ }
+ if (p.getPDGID()==443) {
+ aida.cloud1D("J/psi mass").fill(p.getMass());
+ }
+ }
+ }
+
+ public void bTrackPlots(EventHeader event, MCParticle q, String qname) {
+ List<LCRelation> TracksToMCP = event.get(LCRelation.class,"TracksToMCP");
+ List<MCParticle> Children = getAllFSChildren(q);
+ List<LCRelation> ChildTracks = new ArrayList<LCRelation>();
+
+ for (LCRelation rel:TracksToMCP) {
+ if (Children.contains((MCParticle) rel.getTo())) {
+ ChildTracks.add(rel);
+ }
+ }
+
+ Hep3Vector p_trans = new BasicHep3Vector(q.getMomentum().x(),q.getMomentum().y(),0.0);
+
+ double sumE = 0.0;
+ double sumPx = 0.0;
+ double sumPy = 0.0;
+ double sumPz = 0.0;
+
+ for (LCRelation rel:ChildTracks) {
+ ReconTrack trk = (ReconTrack) rel.getFrom();
+ MCParticle prt = (MCParticle) rel.getTo();
+
+ aida.histogram1D("HIST "+qname+" child D0's (trans)",100,-40,40).fill(trk.getTrackParameter(0));
+ aida.histogram1D("HIST "+qname+" child z impact parameter",100,-20,20).fill(trk.getTrackParameter(3));
+ sumE += prt.getEnergy();
+ sumPx += prt.getMomentum().x()*2.99;
+ sumPy += prt.getMomentum().y()*2.99;
+ sumPz += prt.getMomentum().z()*2.99;
+
+ }
+ double realBinvMass = (q.getEnergy()*q.getEnergy()
+ - q.getMomentum().magnitudeSquared())/2.99;
+ double realBinvMassTrans = (q.getEnergy()*q.getEnergy()
+ - p_trans.magnitudeSquared())/2.99;
+
+
+
+ double sumP = Math.sqrt(sumPx*sumPx+sumPy*sumPy+sumPz*sumPz);
+ double sumPt = Math.sqrt(sumPx*sumPx+sumPy*sumPy);
+
+ double fakeBinvMass = Math.sqrt(sumE*sumE - sumP*sumP)/8.98;
+ double fakeBinvMassTrans = Math.sqrt(sumE*sumE - sumPt*sumPt)/8.98;
+
+ aida.cloud1D(qname+" inv mass: recon to real - tracks").fill(fakeBinvMass);
+ aida.cloud1D(qname+" inv mass fraction - trans tracks").fill(fakeBinvMassTrans/realBinvMassTrans);
+
+ }
+
+ 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;
+ }
+
+ 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