lcsim-contrib/src/main/java/org/lcsim/contrib/homer/vvHratio
diff -u -r1.1 -r1.2
--- MyDBDvvH.java 11 Sep 2012 23:41:08 -0000 1.1
+++ MyDBDvvH.java 27 Sep 2012 09:00:41 -0000 1.2
@@ -1,22 +1,56 @@
import java.util.List;
+//import org.lcsim.event.*;
import org.lcsim.event.EventHeader;
import org.lcsim.event.MCParticle;
import org.lcsim.event.ReconstructedParticle;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
+import org.lcsim.recon.postrecon.leptonID.electron.*;
+
+import java.util.*;
+import java.io.IOException;
+
+
public class MyDBDvvH extends Driver {
+
+ String rplnew = "eehReconstructedParticles";
+ String rpl = "ReconstructedParticles";
+
+
private AIDA aida = AIDA.defaultInstance();
public MyDBDvvH(){
+//
+// Overwrite "ReconstructedParticles" with extra electron ID
+//
+ add(new IsolatedHighPElectronIdentifier());
+//
+// Make "eehReconstructedParticles" by adding brem photons to electrons
+//
+// add(new AddBremPhotonsToElectrons());
+
// add(new MyJetFindingDriver());
- add(new PFOJetFindingDriver());
- add(new MyDBDvvHAnalysis());
+ add(new PFOJetFindingDriver());
+ add(new MyDBDvvHAnalysis());
+
+
};
+ protected String mcParticlesSkimmedName = "MCParticlesSkimmed";
+
public void process(EventHeader event) {
- super.process(event); // this takes care that the child Drivers are loaded and processed.
+ List<MCParticle> particles = event.get(MCParticle.class,mcParticlesSkimmedName);
+ int doit=0;
+ // a nu as the first particle is sufficient to identify this as nunuH
+ // DON'T FORGET TO REMOVE FOR BACKGROUNDS!
+ for (MCParticle particle : particles) {
+ if (Math.abs(particle.getPDGID())==12) doit=1;
+ break;
+ }
+ System.out.println("doit = " + doit);
+ if (doit==1)
+ super.process(event); // this takes care that the child Drivers are loaded and processed.
// List<List<MCParticle>> jetList = (List<List<MCParticle>>) event.get("jetList");
- List<List<ReconstructedParticle>> jetList = (List<List<ReconstructedParticle>>) event.get("jetList");
- aida.cloud1D("nJets").fill(jetList.size());
+
}
}
\ No newline at end of file
lcsim-contrib/src/main/java/org/lcsim/contrib/homer/vvHratio
diff -u -r1.1 -r1.2
--- MyDBDvvHAnalysis.java 11 Sep 2012 23:41:08 -0000 1.1
+++ MyDBDvvHAnalysis.java 27 Sep 2012 09:00:41 -0000 1.2
@@ -1,5 +1,6 @@
import org.lcsim.util.aida.AIDA;
import hep.physics.vec.VecOp;
+import hep.physics.vec.Hep3Vector;
import java.util.List;
import org.lcsim.util.Driver;
import org.lcsim.event.*;
@@ -13,13 +14,14 @@
import org.lcsim.spacegeom.SpaceVector;
+
public class MyDBDvvHAnalysis extends Driver
{
+
// collections
protected String inputPfoCollection = "PandoraPFOCollection";
protected String mcParticlesSkimmedName = "MCParticlesSkimmed";
-
// First get the reconstructed particles from the event
List<ReconstructedParticle> PFOs = new ArrayList<ReconstructedParticle>();
@@ -35,38 +37,103 @@
// print(HLEVEL_DEFAULT, e.getMessage(), true);
// }
+ // check for event with reconstruction errors
+ int skipEvent = 0;
+ for (ReconstructedParticle pfo : PFOs) {
+ boolean passPfoSelection = true;
+
+ int type = pfo.getType();
+// double cosTheta = momentum.cosTheta();
+ double energy = pfo.getEnergy();
+ // System.out.println("energy = "+energy + "type = " + type);
+ if (energy>1000.) {
+ skipEvent = 1;
+ continue;
+ }
+ }
+ if (skipEvent == 1) return; // skip this bad event
+
+
List<Track> tracklist = event.getTracks();
- int nTrks = tracklist.size();
+ // int nTrks = tracklist.size();
+
+ List<List<ReconstructedParticle>> jetList = (List<List<ReconstructedParticle>>) event.get("jetList");
+ int nJets = jetList.size();
+ aida.cloud1D("Trk #all Tracks").fill(tracklist.size());
+ aida.cloud1D("#Jets from PFO in ALL Events").fill(jetList.size());
- aida.cloud1D("Trk nTracks").fill(nTrks);
-
-// List<MCParticle> particles = event.get(MCParticle.class,event.MC_PARTICLES);
+ // List<MCParticle> particles = event.get(MCParticle.class,event.MC_PARTICLES);
List<MCParticle> particles = event.get(MCParticle.class,mcParticlesSkimmedName);
- aida.cloud1D("Number of MC particles").fill(particles.size());
- for (MCParticle particle : particles) {
- for (MCParticle mcp2 : particle.getParents()) {
- int parentid = mcp2.getPDGID();
- if (Math.abs(parentid)==25) {
- int id = particle.getPDGID();
- aida.cloud1D("All MC energy").fill(particle.getEnergy());
- aida.cloud1D("All MC cosTheta").fill(VecOp.cosTheta(particle.getMomentum()));
- aida.cloud1D("All MC phi").fill(VecOp.phi(particle.getMomentum()));
- aida.cloud1D("All MC H daughter ID").fill(id);
- System.out.println("Parent type:"+mcp2.getType().getName());
- System.out.println("Daughter type:"+particle.getType().getName());
-// System.out.println("Starting point: "+mcp2.getOrigin().toString());
-// System.out.println("Momentum: "+mcp2.getMomentum());
-// System.out.println("Endpoint: "+mcp2.getEndPoint());
-// System.out.println("# Daughters: "+mcp2.getDaughters().size());
- System.out.println("Daughter id:"+id);
- }
- }
- }
-
+ aida.cloud1D("Number of MC particles").fill(particles.size());
+ double MCM2vis = 0.0;
+ double EsumMC = 0.0;
+ double PsumMCX = 0.0;
+ double PsumMCY = 0.0;
+ double PsumMCZ = 0.0;
+
+ int nid=0;
+ for (MCParticle particle : particles) {
+ // System.out.println("Sim stat: " + nid + " " + (particle.getSimulatorStatus()).getValue());
+ if ((particle.getSimulatorStatus()).getValue()==67108864) {
+ // ||
+ // (particle.getSimulatorStatus()).getValue()==16777216 ) {
+ double[] momentum = (particle.getMomentum()).v();
+ double PX = momentum[0];
+ double PY = momentum[1];
+ double PZ = momentum[2];
+ double PT = Math.sqrt(PX*PX + PY*PY);
+ double P = Math.sqrt(PT*PT+PZ*PZ);
+ aida.cloud1D("All MC particle cosTheta").fill(VecOp.cosTheta(particle.getMomentum()));
+ // System.out.println("Sim stat: " + nid + " " + (particle.getSimulatorStatus()).getValue() + " " + PX + " " + PY + " " + PZ + " " + P);
+ if (Math.abs(VecOp.cosTheta(particle.getMomentum()))<0.95) {
+ EsumMC += P;
+ PsumMCX += PX;
+ PsumMCY += PY;
+ PsumMCZ += PZ;
+ // System.out.println("Sums: " + nid + " " + (particle.getSimulatorStatus()).getValue() + " " + PsumMCX + " " + PsumMCY + " " + PsumMCZ + " " + EsumMC);
+ }
+ }
+ //don't want
+ //33554432 left
+ //1107296256 left
+ //134217728 tracker
+ //want
+ //67108864 cal
+ //67108864 cal
+ //16777216 stopped
+
+
+ for (MCParticle mcp2 : particle.getParents()) {
+ int parentid = mcp2.getPDGID();
+ if (Math.abs(parentid)==25) {
+ int id = particle.getPDGID();
+ aida.cloud1D("All MC H daughters energy").fill(particle.getEnergy());
+ aida.cloud1D("All MC H daughters cosTheta").fill(VecOp.cosTheta(particle.getMomentum()));
+ aida.cloud1D("All MC H daughters phi").fill(VecOp.phi(particle.getMomentum()));
+ aida.cloud1D("All MC H daughter ID").fill(id);
+ System.out.println("Parent type:"+mcp2.getType().getName());
+ System.out.println("Daughter type:"+particle.getType().getName());
+ // System.out.println("Starting point: "+mcp2.getOrigin().toString());
+ // System.out.println("Momentum: "+mcp2.getMomentum());
+ // System.out.println("Endpoint: "+mcp2.getEndPoint());
+ // System.out.println("# Daughters: "+mcp2.getDaughters().size());
+ System.out.println("Daughter id:"+id);
+ }
+ }
+ nid++;
+ }
+ double EvisMC = EsumMC;
+ double PTvisMC = Math.sqrt(PsumMCX*PsumMCX+PsumMCY*PsumMCY);
+ double M2visMC = EsumMC*EsumMC - PTvisMC*PTvisMC - PsumMCZ*PsumMCZ;
+ aida.cloud1D("Evis from MC").fill(EvisMC);
+ aida.cloud1D("PTvis from MC").fill(PTvisMC);
+ aida.cloud1D("M2vis from MC").fill(M2visMC);
+ if (M2visMC > 0.) aida.cloud1D("Mvis from MC").fill(Math.sqrt(M2visMC));
// --------------------------------------------------
+ int nTrks = 0;
double EsumTrk = 0.0;
double PsumTrkX = 0.0;
double PsumTrkY = 0.0;
@@ -75,38 +142,42 @@
// aida.cloud1D("Trk cosTheta").fill(track.getMomentum());
TrackState trkstate = track.getTrackStates().get(0);
double[] momentum = trkstate.getMomentum();
- System.out.println("mom = " + momentum[0] + "," + momentum[1] + "," + momentum[2]);
+ // System.out.println("mom = " + momentum[0] + "," + momentum[1] + "," + momentum[2]);
double PX = momentum[0];
double PY = momentum[1];
double PZ = momentum[2];
double PT = Math.sqrt(PX*PX + PY*PY);
double P = Math.sqrt(PT*PT+PZ*PZ);
- EsumTrk += P;
- PsumTrkX += PX;
- PsumTrkY += PY;
- PsumTrkZ += PY;
// System.out.println("PT = "+PT);
-// System.out.println("P = "+P);
+// System.out.println("P = "+P+" EsumTrk = "+EsumTrk);
double q = trkstate.getOmega()/Math.abs(trkstate.getOmega());
// System.out.println("q = "+q);
double tlam = trkstate.getTanLambda();
// System.out.println("tlam = "+tlam);
aida.cloud1D("Trk cosTheta").fill(PZ/P);
-
+ if (Math.abs(PZ/P)<0.95) {
+ if (P>0.5) nTrks++;
+ EsumTrk += P;
+ PsumTrkX += PX;
+ PsumTrkY += PY;
+ PsumTrkZ += PZ;
+ }
}
double EvisTrk = EsumTrk;
double PTvisTrk = Math.sqrt(PsumTrkX*PsumTrkX+PsumTrkY*PsumTrkY);
double M2visTrk = EsumTrk*EsumTrk - PTvisTrk*PTvisTrk - PsumTrkZ*PsumTrkZ;
+ aida.cloud1D("#good Tracks").fill(nTrks);
aida.cloud1D("Evis from Tracks").fill(EvisTrk);
aida.cloud1D("PTvis from Tracks").fill(PTvisTrk);
- aida.cloud1D("M2vis from Tracks").fill(M2visTrk);
-
+ // aida.cloud1D("M2vis from Tracks").fill(M2visTrk);
+ if (M2visTrk>0.) aida.cloud1D("Mvis from Tracks").fill(Math.sqrt(M2visTrk));
- if ((PTvisTrk>20.0 && PTvisTrk<500.0) &&
- (EvisTrk>100.0 && EvisTrk<400.0) &&
- (nTrks>=6 && nTrks<=19)) {
+ if ((PTvisTrk>20.0 && PTvisTrk<250.0) &&
+ (EvisTrk>100.0 && EvisTrk<500.0) &&
+ (nTrks>=6 && nTrks<=45)) {
aida.cloud1D("Select M2vis from Tracks").fill(M2visTrk);
+ aida.cloud1D("#Jets from PFO in Track Selected Events").fill(jetList.size());
}
// -----------------------------------------------------------------------------
double EsumPFO = 0.0;
@@ -120,39 +191,50 @@
int type = pfo.getType();
SpaceVector momentum = new SpaceVector(pfo.getMomentum());
double pT = momentum.rxy();
- double p = momentum.rxyz();
+ double P = momentum.rxyz();
double PX = momentum.x();
double PY = momentum.y();
double PZ = momentum.z();
double cosTheta = momentum.cosTheta();
double energy = pfo.getEnergy();
- EsumPFO += energy;
- PsumPFOX += PX;
- PsumPFOY += PY;
- PsumPFOZ += PZ;
- eTotalInput += energy;
-// List<Cluster> pfoclusters = pfo.getClusters();
-// List<Track> pfotracks = pfo.getTracks();
+ // System.out.println("energy = "+energy + "type = " + type);
+ if (Math.abs(PZ/P)<0.95 && energy<1000.) {
+ EsumPFO += energy;
+ PsumPFOX += PX;
+ PsumPFOY += PY;
+ PsumPFOZ += PZ;
+ eTotalInput += energy;
+ // List<Cluster> pfoclusters = pfo.getClusters();
+ // List<Track> pfotracks = pfo.getTracks();
+ }
}
double EvisPFO = EsumPFO;
double PTvisPFO = Math.sqrt(PsumPFOX*PsumPFOX+PsumPFOY*PsumPFOY);
double M2visPFO = EsumPFO*EsumPFO - PTvisPFO*PTvisPFO - PsumPFOZ*PsumPFOZ;
aida.cloud1D("Evis from PFO").fill(EvisPFO);
aida.cloud1D("PTvis from PFO").fill(PTvisPFO);
- aida.cloud1D("M2vis from PFO").fill(M2visPFO);
+ if (M2visPFO>0.) aida.cloud1D("Mvis from PFO").fill(Math.sqrt(M2visPFO));
// ------------------MC parent info for PFO selected events ----------------------------------------
- if ((PTvisPFO>20.0 && PTvisPFO<500.0) &&
- (EvisPFO>100.0 && EvisPFO<400.0) &&
- (nTrks>=6 && nTrks<=19)) {
- aida.cloud1D("Select M2vis from PFO").fill(M2visPFO);
+ if ((PTvisPFO>20.0 && PTvisPFO<250.0) &&
+ (EvisPFO>100.0 && EvisPFO<500.0) &&
+ (nTrks>=6 && nTrks<=45)) {
+ aida.cloud1D("Select M2vis from PFO").fill(M2visPFO);
+ aida.cloud1D("#Jets from PFO in PFO Selected Events").fill(nJets);
+ if (nJets<4) aida.cloud1D("Select M2vis from PFO for #jets<4").fill(M2visPFO);
+ if (nJets<4 && M2visPFO>0.) {
+ aida.cloud1D("Select Mvisible from PFO for #jets<4").fill(Math.sqrt(M2visPFO));
+// aida.cloud1D("Select Mvisible from PFO for #jets<4 scaled using Evis").fill(Math.sqrt(M2visPFO)*500.0/EvisPFO);
+ }
for (MCParticle particle : particles) {
for (MCParticle mcp2 : particle.getParents()) {
int parentid = mcp2.getPDGID();
if (Math.abs(parentid)==25) {
int id = particle.getPDGID();
aida.cloud1D("PFO Selected MC H daughter ID").fill(id);
+ if (nJets<4) aida.cloud1D("PFO Selected MC H daughter ID for #jets<4").fill(id);
+ if (nJets<4 && M2visPFO>0. && Math.abs(id)==5) aida.cloud1D("Select Mvisible from PFO for #jets<4 for H-bb").fill(Math.sqrt(M2visPFO));
}
}
}