lcsim-contrib/src/main/java/org/lcsim/contrib/homer/vvHratio/sid/java
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);
}
}
}