Commit in lcsim-contrib/src/main/java/org/lcsim/contrib/homer/vvHratio/sid/java on MAIN | |||
MyDBDvvH.java | +20 | -7 | 1.1 -> 1.2 |
MyDBDvvHAnalysis.java | +437 | -98 | 1.1 -> 1.2 |
PFOJetFindingDriver.java | +2 | 1.1 -> 1.2 | |
+459 | -105 |
update of the ffH n-tuple production code
diff -u -r1.1 -r1.2 --- MyDBDvvH.java 11 Oct 2012 16:17:22 -0000 1.1 +++ MyDBDvvH.java 23 Apr 2013 23:33:58 -0000 1.2 @@ -1,3 +1,4 @@
+package org.lcsim.contrib.homer;
import java.util.List; //import org.lcsim.event.*;
@@ -7,6 +8,8 @@
import org.lcsim.util.Driver; import org.lcsim.util.aida.AIDA; import org.lcsim.recon.postrecon.leptonID.electron.*;
+import org.lcsim.contrib.homer.PFOJetFindingDriver; +import org.lcsim.contrib.homer.MyDBDvvHAnalysis;
import java.util.*; import java.io.IOException;
@@ -24,7 +27,9 @@
// // Overwrite "ReconstructedParticles" with extra electron ID //
+
add(new IsolatedHighPElectronIdentifier());
+
// // Make "eehReconstructedParticles" by adding brem photons to electrons //
@@ -32,6 +37,7 @@
// add(new MyJetFindingDriver()); add(new PFOJetFindingDriver());
+
add(new MyDBDvvHAnalysis());
@@ -39,18 +45,25 @@
protected String mcParticlesSkimmedName = "MCParticlesSkimmed"; public void process(EventHeader event) {
+// System.out.println("RUN = " + event.getRunNumber());
List<MCParticle> particles = event.get(MCParticle.class,mcParticlesSkimmedName);
- int doit=0;
+ int isanu=0; + int isahiggs=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;
+ if (Math.abs(particle.getPDGID())==12) isanu=1; + if (Math.abs(particle.getPDGID())==25) isahiggs=1; + if (isanu==1 && isahiggs==1) break;
}
- System.out.println("doit = " + doit); - if (doit==1) - super.process(event); // this takes care that the child Drivers are loaded and processed.
+ System.out.println("isahiggs = " + isahiggs + " , isanu = " + isanu); +// if (isahiggs==0 || (isahiggs==1 && isanu==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"); }
-}
\ No newline at end of file
+ + protected void endOfData() { + super.endOfData(); + } +}
diff -u -r1.1 -r1.2 --- MyDBDvvHAnalysis.java 11 Oct 2012 16:17:22 -0000 1.1 +++ MyDBDvvHAnalysis.java 23 Apr 2013 23:33:58 -0000 1.2 @@ -1,6 +1,12 @@
-import org.lcsim.util.aida.AIDA;
+package org.lcsim.contrib.homer; +//import org.lcsim.util.aida.AIDA; +import hep.aida.*; +import org.lcsim.util.aida.*;
import hep.physics.vec.VecOp; import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHepLorentzVector; +import hep.physics.vec.HepLorentzVector; +import hep.physics.jet.EventShape;
import java.util.List; import org.lcsim.util.Driver; import org.lcsim.event.*;
@@ -13,16 +19,34 @@
import org.lcsim.spacegeom.SpacePoint; import org.lcsim.spacegeom.SpaceVector;
-
+import java.nio.*; +import java.io.*;
public class MyDBDvvHAnalysis extends Driver { // collections
+ protected String RefinedJetsrelName = "RefinedJets_rel"; + protected String RefinedJetsName = "RefinedJets";
protected String inputPfoCollection = "PandoraPFOCollection"; protected String mcParticlesSkimmedName = "MCParticlesSkimmed";
+ IHistogram1D evisjets; + IHistogram1D ptvisjets; + + IHistogram1D mvisjets; + IHistogram1D mvispfo; + IHistogram1D mvistrks; + IHistogram1D mvismc; +
protected double E_BEAM = 500.0;
+ protected int ikt = 0; + protected int ievt = 0; +// protected double ktvals[] = {0.70,0.70,1.00,1.50}; + protected double ktvals[] = {1.50,1.50,1.50,1.50}; + protected int iflav[] = {0,0,0,0}; + protected int badidrup[]={35344,35342,35356,35354,35338,35400,35398,35436,35434,35394,35372,35370,35384,35382,35368,35516,35514,35552,35550,35512}; +
// First get the reconstructed particles from the event List<ReconstructedParticle> PFOs = new ArrayList<ReconstructedParticle>();
@@ -30,9 +54,34 @@
private AIDA aida = AIDA.defaultInstance();
+ protected FileWriter fstream = null; + protected BufferedWriter out = null; + + protected void endOfData() + { + if (fstream != null) { + try { + out.close(); + } catch (IOException e) { + e.printStackTrace(); + } + } + } +
protected void process(EventHeader event) {
-
+ ievt++; + if (fstream == null) { + try { + fstream = new FileWriter("TMVA-vars.txt"); + out = new BufferedWriter(fstream); + } catch (IOException e) { + e.printStackTrace(); + } + } + + double evwgt = event.getWeight(); +
// try { PFOs = event.get(ReconstructedParticle.class, inputPfoCollection); // } catch (IllegalArgumentException e) {
@@ -47,29 +96,125 @@
int type = pfo.getType(); // double cosTheta = momentum.cosTheta(); double energy = pfo.getEnergy();
- // System.out.println("energy = "+energy + "type = " + type);
+// System.out.println("energy = "+energy + "type = " + type);
if (energy>1000.) { skipEvent = 1; continue; } }
- // if (skipEvent == 1) return; // skip this bad event
+ // if (skipEvent == 1) return; // skip this bad event + if (event.getRunNumber()==17267118 && event.getEventNumber()==0) ikt++; + double ktval = ktvals[ikt];
+ int idrup_value=0; + if(event.getIntegerParameters().containsKey("idrup")) + idrup_value=event.getIntegerParameters().get("idrup")[0]; + + skipEvent = 0; + for (int irup=0;irup<badidrup.length;irup++) { + if (idrup_value==badidrup[irup]) { + skipEvent = 1; + continue; + } + } + if (skipEvent == 1) { + System.out.println("skipping: bad event, idrup="+idrup_value); + return; + }
List<Track> tracklist = event.getTracks();
- // int nTrks = tracklist.size();
+ // int nTrks = tracklist.size();
// List<List<ReconstructedParticle>> jetList = (List<List<ReconstructedParticle>>) event.get("JetOut");
- List<ReconstructedParticle> jetList = (List<ReconstructedParticle>)event.get("JetOut");
+ + List<ReconstructedParticle> jetList = null; + try { + jetList = (List<ReconstructedParticle>)event.get("JetOut"); + } catch (IllegalArgumentException e) { + // print(HLEVEL_DEFAULT, e.getMessage(), true); + System.out.println("skipping: no JetOut"); + return; + } +
int nJets = jetList.size();
-
+ + List<ReconstructedParticle> jetPartsList = null; + try { + jetPartsList = (List<ReconstructedParticle>)event.get("JetParts2Jet"); + } catch (IllegalArgumentException e) { + // print(HLEVEL_DEFAULT, e.getMessage(), true); + return; + } + + int nJetsParts= jetPartsList.size(); + + List<ReconstructedParticle> colJet = null; + try { + colJet = (List<ReconstructedParticle>)event.get(RefinedJetsName); + } catch (IllegalArgumentException e) { + // print(HLEVEL_DEFAULT, e.getMessage(), true); + System.out.println("skipping: no refined jets"); + return; + } + + + List<LCRelation>colJetrel = new ArrayList<LCRelation>(); + try { + colJetrel = event.get(LCRelation.class, RefinedJetsrelName); + } catch (IllegalArgumentException e) { + // print(HLEVEL_DEFAULT, e.getMessage(), true); + return; + } + +
aida.cloud1D("Trk #all Tracks").fill(tracklist.size()); aida.cloud1D("#Jets from PFO in ALL Events").fill(jetList.size()); // 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());
+ aida.cloud1D("MC:/Number of MC particles").fill(particles.size());
+
+ int did = 0;
+ int hasanu = 0;
+ // classify the event by flavor
+ for (MCParticle particle : particles) {
+
+ if (Math.abs(particle.getPDGID())==12) hasanu=1;
+
+ for (MCParticle mcp2 : particle.getParents()) {
+ int parentid = mcp2.getPDGID();
+
+ if (Math.abs(parentid)==25) {
+ int id = particle.getPDGID();
+ if (Math.abs(id)!=25) {
+ aida.cloud1D("MC:/All MC H daughters energy").fill(particle.getEnergy());
+ aida.cloud1D("MC:/All MC H daughters cosTheta").fill(VecOp.cosTheta(particle.getMomentum()));
+ aida.cloud1D("MC:/All MC H daughters phi").fill(VecOp.phi(particle.getMomentum()));
+ aida.cloud1D("MC:/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);
+
+ did = Math.abs(id);
+ }
+ }
+ // System.out.println("H decay type identified as:"+" ( Flavor = "+did+" )");
+ }
+
+ }
+ // ---------------------------------------------
+ if (did==5) iflav[0]++;
+ if (did==4) iflav[1]++;
+ if (did==24) iflav[2]++;
+ if (did==21) iflav[3]++;
+ System.out.println("bb= " + iflav[0] + " cc= " + iflav[1] + " ww= " + iflav[2] + " gg= " + iflav[3]);
+
+
double MCM2vis = 0.0; double EsumMC = 0.0; double PsumMCX = 0.0;
@@ -87,18 +232,27 @@
// || // (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;
+ + if (Math.abs(VecOp.cosTheta(particle.getMomentum()))<0.96) { + + 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("MC:/final state MC particle PT").fill(PT); + aida.cloud1D("MC:/final state MC particle P").fill(P); + aida.cloud1D("MC:/final state 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); // } }
@@ -112,35 +266,20 @@
//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));
+ aida.cloud1D("MC:/Evis from MC"+" ( Flavor = "+did+" )").fill(EvisMC); + aida.cloud1D("MC:/PTvis from MC"+" ( Flavor = "+did+" )").fill(PTvisMC); + aida.cloud1D("MC:/M2vis from MC"+" ( Flavor = "+did+" )").fill(M2visMC); + mvismc = aida.histogram1D("MC:/Mvis from MC ( Flavor = "+did+" )", 100, 1., 500.); + if (M2visMC > 0.) mvismc.fill(Math.sqrt(M2visMC));
// -------------------------------------------------- int nTrks = 0;
+ int nChTrks = 0; // # of charged tracks
double EsumTrk = 0.0; double PsumTrkX = 0.0; double PsumTrkY = 0.0;
@@ -162,29 +301,34 @@
// System.out.println("q = "+q); double tlam = trkstate.getTanLambda(); // System.out.println("tlam = "+tlam);
- aida.cloud1D("Trk cosTheta").fill(PZ/P);
+ aida.cloud1D("TRKS:/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;
- }
+ } + if (q!=0) { + aida.cloud1D("TRKS:/ChTrk cosTheta").fill(PZ/P); + nChTrks++; + } +
} 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("TRKS:/#good Tracks ( Flavor = "+did+" )").fill(nTrks); + aida.cloud1D("TRKS:/Evis from Tracks").fill(EvisTrk); + aida.cloud1D("TRKS:/PTvis from Tracks").fill(PTvisTrk);
// aida.cloud1D("M2vis from Tracks").fill(M2visTrk);
- if (M2visTrk>0.) aida.cloud1D("Mvis from Tracks").fill(Math.sqrt(M2visTrk));
+ mvistrks = aida.histogram1D("TRKS:/Mvis from Tracks ( Flavor = "+did+" )", 100, 1., 500.); + if (M2visTrk>0.) mvistrks.fill(Math.sqrt(M2visTrk));
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());
+ aida.cloud1D("TRKS:/Select M2vis from Tracks").fill(M2visTrk);
} // ---------------------PFO----------------------------------------------------- double EsumPFO = 0.0;
@@ -211,6 +355,16 @@
if (Math.abs(pidu)==11) nelec++; }
+ + /* + if (P>50.0) { // only high P electrons + if (pfo.getParticleIDUsed() != null) { + int pidu = (pfo.getParticleIDUsed()).getPDG(); + if (Math.abs(pidu)==11) nelec++; + } else if (Math.abs(type)==11) nelec++; + } + */ +
// System.out.println("energy = "+energy + "type = " + type); if (Math.abs(PZ/P)<0.95 && energy<1000.) { EsumPFO += energy;
@@ -225,52 +379,105 @@
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); - if (M2visPFO>0.) aida.cloud1D("Mvis from PFO").fill(Math.sqrt(M2visPFO)); - aida.cloud1D("Number of electrons per event").fill(nelec);
+ aida.cloud1D("PFO:/Evis from PFO ( Flavor = "+did+" )").fill(EvisPFO); + aida.cloud1D("PFO:/PTvis from PFO ( Flavor = "+did+" )").fill(PTvisPFO); + mvispfo = aida.histogram1D("PFO:/Mvis from PFO ( Flavor = "+did+" )", 100, 1., 500.); + double MvisPFO = 0.; + if (M2visPFO>0.) MvisPFO = Math.sqrt(M2visPFO); + mvispfo.fill(MvisPFO); + aida.cloud1D("Number of electrons per event ( Flavor = "+did+" )").fill(nelec);
// ------------------MC parent info for PFO selected events ---------------------------------------- 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);
+ aida.cloud1D("PFO:/Select M2vis from PFO").fill(M2visPFO); + aida.cloud1D("PFO:/#Jets from PFO in PFO Selected Events").fill(nJets); + if (M2visPFO>0.) { + aida.cloud1D("PFO:/Mvis from PFO w/selection").fill(Math.sqrt(M2visPFO));
} 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<6) aida.cloud1D("PFO Selected MC H daughter ID for #jets<6").fill(id); - if (nJets<6 && nelec==0) aida.cloud1D("PFO Selected MC H daughter ID for #jets<6 and nelec=0").fill(id); - if (nJets<6 && nelec>=3) aida.cloud1D("PFO Selected MC H daughter ID for #jets<6 and nelec>=3").fill(id); - if (nJets<6 && M2visPFO>0. && Math.abs(id)==5) aida.cloud1D("Select Mvisible from PFO for #jets<6 for H-bb").fill(Math.sqrt(M2visPFO));
+ aida.cloud1D("PFO:/PFO Selected MC H daughter ID").fill(id); + if (nelec==0) aida.cloud1D("PFO:/PFO Selected MC H daughter ID nelec=0").fill(id); + if (nelec>=3) aida.cloud1D("PFO:/PFO Selected MC H daughter ID nelec>=3").fill(id);
} } } }
+ // ----- flavor tag info ----- + double bprob[] = {0.0,0.0}; + double cprob[] = {0.0,0.0}; + int nvtxall[] = {0,0}; + int nvtx[] = {0,0}; + double btag[]=null; + int ijet=-1; + for (ReconstructedParticle rjet : colJet) { + ijet++; + // rjet.getAlgorithmID("lcfiplus"); + // double btag = rjet.get("BTag"); + // System.out.println("rjet ParticleIDs being checked"); + + for (LCRelation rjetrel : colJetrel) { + ReconstructedParticle rjetcmp = (ReconstructedParticle) rjetrel.getFrom(); + Vertex vtx = (Vertex) rjetrel.getTo(); + if (rjetcmp == rjet) { + nvtxall[ijet]++; + if (vtx.getChi2()>0.0) nvtx[ijet]++; + } + } + // System.out.println("ijet = "+ijet+" nvtxall = "+nvtxall[ijet]+" nvtx = "+nvtx[ijet]); +
+ if (rjet.getParticleIDs() != null) { + // System.out.println("rjet ParticleIDs !NULL"); + List<ParticleID> parts = rjet.getParticleIDs(); + // List<String> parmns = rjet.getParameterName(); + for (ParticleID pidh : parts) { + btag=null; + btag = pidh.getParameters(); + // btag = pidh.getDoubleParameters(); + // System.out.println("btag size = "+btag.length); + if (btag != null) { + int algo = pidh.getAlgorithmType(); + // System.out.println("algo type = "+algo); + if (btag.length==3) { + bprob[ijet] = btag[0]; + cprob[ijet] = btag[1]; + } + } + } + } + aida.cloud1D("FLAV:/BPROB ( Flavor = "+did+" )").fill(bprob[ijet]); + aida.cloud1D("FLAV:/CPROB ( Flavor = "+did+" )").fill(cprob[ijet]);
+ }
// -------------JETS------------------------------------------------------ double EsumJETS = 0.0; double PsumJETSX = 0.0; double PsumJETSY = 0.0; double PsumJETSZ = 0.0; double eJETTotalInput = 0.0;
- // int nelec = 0;
+ + double m2jet = 0.0; + double mjet[] = {0.0, 0.0}; + double phijet[] = {0.0, 0.0}; + double cjet[] = {0.0, 0.0}; + double ejet[] = {0.0, 0.0}; + double qjet[] = {0.0, 0.0}; + + int iijet = -1; +
boolean passjetSelection = true; for (ReconstructedParticle jet : jetList) {
+ iijet++;
int ntrksjet = (jet.getParticles()).size();
- aida.cloud1D("#tracks in jet").fill(ntrksjet);
+ aida.cloud1D("JETS:/#tracks in jet"+" ( Flavor = "+did+" )").fill(ntrksjet);
- // int type = pfo.getType();
SpaceVector momentum = new SpaceVector(jet.getMomentum()); double pT = momentum.rxy(); double P = momentum.rxyz();
@@ -278,59 +485,191 @@
double PY = momentum.y(); double PZ = momentum.z(); double cosTheta = momentum.cosTheta();
+ double phi =momentum.phi();
double energy = jet.getEnergy();
+ ejet[iijet] = energy;
- // if (pfo.getParticleIDUsed() != null) { - // int pidu = (pfo.getParticleIDUsed()).getPDG(); - // if (Math.abs(pidu)==11) nelec++; - // }
+ qjet[iijet] = jet.getCharge();
// System.out.println("energy = "+energy + "type = " + type);
- if (Math.abs(PZ/P)>0.90 || energy>(2.0*E_BEAM)) {
+ if (Math.abs(PZ/P)>0.75 || energy>(2.0*E_BEAM)) {
passjetSelection = false; } else {
- aida.cloud1D("good JETS cos theta").fill(cosTheta);
+ aida.cloud1D("JETS:/good JETS cos theta").fill(cosTheta);
} EsumJETS += energy; PsumJETSX += PX; PsumJETSY += PY; PsumJETSZ += PZ; eJETTotalInput += energy;
- // List<Cluster> pfoclusters = pfo.getClusters(); - // List<Track> pfotracks = pfo.getTracks(); - // } - aida.cloud1D("all JETS cos theta").fill(cosTheta);
+ + m2jet = energy*energy-(PX*PX+PY*PY+PZ*PZ); + mjet[iijet] = 0.0; + if (m2jet>0.0) mjet[iijet] = Math.sqrt(m2jet); + + phijet[iijet] = phi; + cjet[iijet] = cosTheta; + + + aida.cloud1D("JETS:/all JETS cos theta").fill(cosTheta);
} double EvisJETS = EsumJETS; double PTvisJETS = Math.sqrt(PsumJETSX*PsumJETSX+PsumJETSY*PsumJETSY); double M2visJETS = EsumJETS*EsumJETS - PTvisJETS*PTvisJETS - PsumJETSZ*PsumJETSZ;
- aida.cloud1D("Evis from JETS").fill(EvisJETS); - aida.cloud1D("PTvis from JETS").fill(PTvisJETS); - if (M2visJETS>0.) aida.cloud1D("Mvis from JETS").fill(Math.sqrt(M2visJETS)); - // aida.cloud1D("Number of electrons per event").fill(nelec); -// ------------------MC parent info for JETS selected events ----------------------------------------
+ double MvisJETS = 0.0; + if (M2visJETS>0.0) MvisJETS = Math.sqrt(M2visJETS); + + evisjets = aida.histogram1D("JETS:/Evis from JETS ( Flavor = "+did+" )", 200, 1., 1000.); + evisjets.fill(EvisJETS); + + ptvisjets = aida.histogram1D("JETS:/PTvis from JETS ( Flavor = "+did+" )", 100, 1., 500.); + ptvisjets.fill(PTvisJETS); + + mvisjets = aida.histogram1D("JETS:/Mvis from JETS (Kt "+ktval +") ( Flavor = "+did+" )", 100, 1., 500.); + mvisjets.fill(MvisJETS); + + mvisjets = aida.histogram1D("JETS:/Mvis from JETS (Kt "+ktval +") - ALL ", 100, 1., 500.); + mvisjets.fill(MvisJETS); + + if (did>0) { + mvisjets = aida.histogram1D("JETS:/Mvis from JETS (Kt "+ktval +") - signal ", 100, 1., 500.); + mvisjets.fill(MvisJETS); + } else { + mvisjets = aida.histogram1D("JETS:/Mvis from JETS (Kt "+ktval +") - background ", 100, 1., 500.); + mvisjets.fill(MvisJETS); + } + // ------ for mumu info + double ejpart[] = { 0.0 , 0.0 }; + int tjpart[] = { 0 , 0 }; + int njparts = 0; + for (ReconstructedParticle jetPart : jetPartsList) { + int type = jetPart.getType(); + double energy = jetPart.getEnergy(); + if (energy>ejpart[0]) { + ejpart[1] = ejpart[0]; + tjpart[1] = tjpart[0]; + ejpart[0] = energy; + tjpart[0] = type; + } else if (energy>ejpart[1]) { + ejpart[1] = energy; + tjpart[1] = type; + } + if (energy>20.) njparts++; + } + List<HepLorentzVector> recothrustlist =new ArrayList<HepLorentzVector>(); + + for (ReconstructedParticle jetPart : jetPartsList) { + recothrustlist.add(jetPart.asFourVector()); + } + EventShape myrecoshape = new EventShape(); + myrecoshape.setEvent(recothrustlist); + double jetthrust = myrecoshape.thrust().magnitude(); + +// ------------------ h decay mode selections ----------------------------------------
-// if ((PTvisJETS>20.0 && PTvisJETS<250.0) && -// (EvisJETS>100.0 && EvisJETS<(2.0*E_BEAM)) && -// (nTrks>=6 && nTrks<=45)) { - if (passjetSelection) { - // aida.cloud1D("Select M2vis from JETS").fill(M2visJETS); - aida.cloud1D("#Jets from JETS in JETS Selected Events").fill(nJets); - // aida.cloud1D("Select M2vis from JETS f").fill(M2visJETS); - if (nJets<4 && M2visJETS>0.) { - aida.cloud1D("Mvisible for events with good Jets").fill(Math.sqrt(M2visJETS)); - // aida.cloud1D("Select Mvisible from JETS for #jets<4 scaled using Evis").fill(Math.sqrt(M2visJETS)*500.0/EvisJETS);
+ boolean btagged = false; + if (bprob[0]>0.85 && bprob[1]>0.85) btagged = true; + boolean ctagged = false; + if (cprob[0]>0.6 && cprob[1]>0.6) ctagged = true; + + String line = null; + + int hdecmode = 0; + + + line = String.valueOf(did)+" " + + String.valueOf(nTrks)+" " + + String.valueOf(bprob[0])+" " + + String.valueOf(bprob[1])+" " + + String.valueOf(cprob[0])+" " + + String.valueOf(cprob[1])+" " + + String.valueOf(nelec)+" " + + String.valueOf(MvisJETS)+" " + + String.valueOf(MvisPFO)+" " + + String.valueOf(mjet[0])+" " + + String.valueOf(mjet[1])+" " + + String.valueOf(PTvisJETS)+" " + + String.valueOf(EvisJETS)+" " + + String.valueOf(event.getRunNumber())+" " + + String.valueOf(evwgt)+" " + + String.valueOf(idrup_value)+" " + + String.valueOf(cjet[0])+" " + + String.valueOf(cjet[1])+" " + + String.valueOf(ejet[0])+" " + + String.valueOf(ejet[1])+" " + + String.valueOf(qjet[0])+" " + + String.valueOf(qjet[1])+" " + + String.valueOf(phijet[1]-phijet[0]+" " + + String.valueOf(ejpart[0])+" " + + String.valueOf(ejpart[1])+" " + + String.valueOf(tjpart[0])+" " + + String.valueOf(tjpart[1])+" " + + String.valueOf(njparts)+" " + + String.valueOf(nJetsParts)+" " + + String.valueOf(nChTrks)+" " + + String.valueOf(jetthrust)+" " + + String.valueOf(nvtxall[0])+" " + + String.valueOf(nvtxall[1])+" " + + String.valueOf(nvtx[0])+" " + + String.valueOf(nvtx[1])+" " + + String.valueOf(hasanu)) + ; + + try { + out.write(line); + out.newLine(); + } catch (IOException e) { + e.printStackTrace();
}
+ + System.out.println("writing out variables: ievt="+ievt); + + // if ( (PTvisJETS>20.0 && PTvisJETS<250.0) && + // (EvisJETS>100.0 && EvisJETS<400.0) && (nelec==0)){ + if ( (PTvisJETS>20.0 && PTvisJETS<250.0) && + (EvisJETS>100.0 && EvisJETS<400.0)){ + /* ---------------------------- + line = String.valueOf(did)+" " + + String.valueOf(nTrks)+" " + + String.valueOf(btagged)+" " + + String.valueOf(ctagged)+" " + + String.valueOf(nelec)+" " + + String.valueOf(MvisJETS)+" " + + String.valueOf(MvisPFO)+" " + + String.valueOf(mjet[0])+" " + + String.valueOf(mjet[1])+" " + + String.valueOf(evwgt) + ; + // line = nTrks+btagged+ctagged+nelec+MvisJets+mjet[0]+mjet[1]; + try { + out.write(line); + out.newLine(); + } catch (IOException e) { + e.printStackTrace(); + } + --------------------------- */ + + if ((nTrks>=6 && nTrks<=45) && btagged && !ctagged) + hdecmode =1; // bb (5) + else if ((nTrks>=5 && nTrks<=30) && !btagged && ctagged) + hdecmode =2; // cc (4) + else if ((nTrks>=16 && nTrks<=60) && !btagged && !ctagged) + hdecmode =3; // WW (24) + else if ((nTrks>=11 && nTrks<=70) && !btagged && !ctagged) + hdecmode =4; // gg (21) + + +// if (passjetSelection) { + aida.cloud1D("JETS:/Mvis from JETS w/preselection"+" (Kt "+ktval +") ( Flavor = "+did+" )").fill(MvisJETS); + aida.cloud1D("JETS:/selection#"+hdecmode+"/Mvis from JETS passing selection#"+hdecmode+" (Kt "+ktval +") ( Flavor = "+did+" )").fill(MvisJETS);
for (MCParticle particle : particles) { for (MCParticle mcp2 : particle.getParents()) { int parentid = mcp2.getPDGID(); if (Math.abs(parentid)==25) { int id = particle.getPDGID();
- aida.cloud1D("JETS Selected MC H daughter ID").fill(id); - if (nJets<6) aida.cloud1D("JETS Selected MC H daughter ID for #jets<6").fill(id); - if (nJets<6 && nelec==0) aida.cloud1D("JETS Selected MC H daughter ID for #jets<6 and nelec=0").fill(id); - if (nJets<6 && nelec>=3) aida.cloud1D("JETS Selected MC H daughter ID for #jets<6 and nelec>=3").fill(id); - if (nJets<6 && M2visJETS>0. && Math.abs(id)==5) aida.cloud1D("Select Mvisible from JETS for #jets<6 for H-bb").fill(Math.sqrt(M2visJETS));
+ aida.cloud1D("JETS:/selection#"+hdecmode+"/MC H daughter ID for events passing selection#"+hdecmode).fill(id); + if (nelec==0) aida.cloud1D("JETS:/JETS Selected MC H daughter ID for nelec=0").fill(id); + if (nelec>=3) aida.cloud1D("JETS:/JETS Selected MC H daughter ID for nelec>=3").fill(id);
} } }
diff -u -r1.1 -r1.2 --- PFOJetFindingDriver.java 11 Oct 2012 16:17:22 -0000 1.1 +++ PFOJetFindingDriver.java 23 Apr 2013 23:33:58 -0000 1.2 @@ -1,3 +1,5 @@
+package org.lcsim.contrib.homer; +
import java.util.ArrayList; import java.util.Collections; import java.util.HashMap;
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