 |
 |
| Commit in lcsim-contrib/src/main/java/org/lcsim/contrib/homer/vvHratio/sid on MAIN | |||
| java/MyDBDvvH.java | +56 | added 1.1 | |
| /MyDBDvvHAnalysis.java | +345 | added 1.1 | |
| /MyJetFindingDriver.java | +41 | added 1.1 | |
| /PFOJetFindingDriver.java | +57 | added 1.1 | |
| java-e/AddBremPhotonsToElectrons.java.hn2 | +128 | added 1.1 | |
| /IsolatedHighPElectronIdentifier.java.hn2 | +261 | added 1.1 | |
| /MyDBDvvHAnalysiseid.java | +433 | added 1.1 | |
| /MyDBDvvHeid.java | +56 | added 1.1 | |
| +1377 | |||
update of analysis and lepton ID testing code
diff -N MyDBDvvH.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MyDBDvvH.java 11 Oct 2012 16:17:22 -0000 1.1 @@ -0,0 +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());
+
+
+ };
+ protected String mcParticlesSkimmedName = "MCParticlesSkimmed";
+
+ public void process(EventHeader event) {
+ 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");
+
+ }
+}
\ No newline at end of file
diff -N MyDBDvvHAnalysis.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MyDBDvvHAnalysis.java 11 Oct 2012 16:17:22 -0000 1.1 @@ -0,0 +1,345 @@
+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.*;
+
+import java.util.ArrayList;
+import java.util.Collections;
+
+import org.lcsim.spacegeom.CartesianPoint;
+import org.lcsim.spacegeom.CartesianVector;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.spacegeom.SpaceVector;
+
+
+
+public class MyDBDvvHAnalysis extends Driver
+{
+
+ // collections
+ protected String inputPfoCollection = "PandoraPFOCollection";
+ protected String mcParticlesSkimmedName = "MCParticlesSkimmed";
+
+ protected double E_BEAM = 500.0;
+
+ // First get the reconstructed particles from the event
+ List<ReconstructedParticle> PFOs = new ArrayList<ReconstructedParticle>();
+
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ protected void process(EventHeader event)
+ {
+
+ // try {
+ PFOs = event.get(ReconstructedParticle.class, inputPfoCollection);
+ // } catch (IllegalArgumentException e) {
+ // 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();
+
+ // List<List<ReconstructedParticle>> jetList = (List<List<ReconstructedParticle>>) event.get("JetOut");
+ List<ReconstructedParticle> jetList = (List<ReconstructedParticle>)event.get("JetOut");
+ int nJets = jetList.size();
+
+ 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());
+ 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.getGeneratorStatus() == particle.FINAL_STATE) {
+ int pdgid = particle.getPDGID();
+ int apdgid = Math.abs(pdgid);
+ if (particle.getGeneratorStatus() == particle.FINAL_STATE && !((apdgid == 12 || apdgid == 14 || apdgid == 16) && (particle.getParents().get(0).getPDGID()==pdgid))) {
+ // 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;
+ double PsumTrkZ = 0.0;
+ for (Track track : tracklist) {
+ // 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]);
+
+ 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);
+// System.out.println("PT = "+PT);
+// 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);
+ if (M2visTrk>0.) aida.cloud1D("Mvis from Tracks").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());
+ }
+// ---------------------PFO-----------------------------------------------------
+ double EsumPFO = 0.0;
+ double PsumPFOX = 0.0;
+ double PsumPFOY = 0.0;
+ double PsumPFOZ = 0.0;
+ double eTotalInput = 0.0;
+ int nelec = 0;
+ for (ReconstructedParticle pfo : PFOs) {
+ boolean passPfoSelection = true;
+
+ int type = pfo.getType();
+ SpaceVector momentum = new SpaceVector(pfo.getMomentum());
+ double pT = momentum.rxy();
+ double P = momentum.rxyz();
+ double PX = momentum.x();
+ double PY = momentum.y();
+ double PZ = momentum.z();
+ double cosTheta = momentum.cosTheta();
+ double energy = pfo.getEnergy();
+
+ if (pfo.getParticleIDUsed() != null) {
+ int pidu = (pfo.getParticleIDUsed()).getPDG();
+ if (Math.abs(pidu)==11) nelec++;
+ }
+
+ // 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);
+ if (M2visPFO>0.) aida.cloud1D("Mvis from PFO").fill(Math.sqrt(M2visPFO));
+ aida.cloud1D("Number of electrons per event").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);
+ }
+ 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));
+ }
+ }
+ }
+ }
+
+
+// -------------JETS------------------------------------------------------
+ double EsumJETS = 0.0;
+ double PsumJETSX = 0.0;
+ double PsumJETSY = 0.0;
+ double PsumJETSZ = 0.0;
+ double eJETTotalInput = 0.0;
+ // int nelec = 0;
+ boolean passjetSelection = true;
+ for (ReconstructedParticle jet : jetList) {
+
+ int ntrksjet = (jet.getParticles()).size();
+ aida.cloud1D("#tracks in jet").fill(ntrksjet);
+
+ // int type = pfo.getType();
+ SpaceVector momentum = new SpaceVector(jet.getMomentum());
+ double pT = momentum.rxy();
+ double P = momentum.rxyz();
+ double PX = momentum.x();
+ double PY = momentum.y();
+ double PZ = momentum.z();
+ double cosTheta = momentum.cosTheta();
+ double energy = jet.getEnergy();
+
+ // if (pfo.getParticleIDUsed() != null) {
+ // int pidu = (pfo.getParticleIDUsed()).getPDG();
+ // if (Math.abs(pidu)==11) nelec++;
+ // }
+
+ // System.out.println("energy = "+energy + "type = " + type);
+ if (Math.abs(PZ/P)>0.90 || energy>(2.0*E_BEAM)) {
+ passjetSelection = false;
+ } else {
+ aida.cloud1D("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);
+ }
+ 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 ----------------------------------------
+
+// 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);
+ }
+ 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));
+ }
+ }
+ }
+ }
+
+
+
+
+ }
+}
+
+
diff -N MyJetFindingDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MyJetFindingDriver.java 11 Oct 2012 16:17:22 -0000 1.1 @@ -0,0 +1,41 @@
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+import hep.physics.jet.JadeEJetFinder;
+import hep.physics.jet.JetFinder;
+import hep.physics.vec.HepLorentzVector;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.Driver;
+
+public class MyJetFindingDriver extends Driver {
+ private double yCut = 0.001;
+
+ public MyJetFindingDriver() {
+ }
+
+ public void process(EventHeader event) {
+ List<MCParticle> particles = event.getMCParticles();
+ Map<HepLorentzVector, MCParticle> map = new HashMap<HepLorentzVector, MCParticle>();
+ for (MCParticle particle : particles) {
+ if (particle.getCharge() != 0) {
+ map.put(particle.asFourVector(), particle);
+ }
+ }
+ JetFinder jetFind = new JadeEJetFinder(yCut);
+ List<List<MCParticle>> jetList = new ArrayList<List<MCParticle>>();
+ jetFind.setEvent(map.keySet());
+ for (int i = 0; i < jetFind.njets(); ++i) {
+ List<HepLorentzVector> jetVecList = jetFind.particlesInJet(i);
+ List<MCParticle> jetPartList = new ArrayList<MCParticle>();
+ for (HepLorentzVector iVec : jetVecList) {
+ jetPartList.add(map.get(iVec));
+ }
+ jetList.add(jetPartList);
+ }
+ event.put("jetList", jetList);
+ }
+}
diff -N PFOJetFindingDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ PFOJetFindingDriver.java 11 Oct 2012 16:17:22 -0000 1.1 @@ -0,0 +1,57 @@
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+import hep.physics.jet.JadeEJetFinder;
+import hep.physics.jet.JetFinder;
+import hep.physics.vec.HepLorentzVector;
+import hep.physics.vec.VecOp;
+
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.Driver;
+import org.lcsim.event.*;
+
+import org.lcsim.spacegeom.CartesianPoint;
+import org.lcsim.spacegeom.CartesianVector;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.spacegeom.SpaceVector;
+
+
+public class PFOJetFindingDriver extends Driver {
+ private double yCut = 0.02;
+ protected String inputPfoCollection = "PandoraPFOCollection";
+
+
+ // First get the reconstructed particles from the event
+ List<ReconstructedParticle> PFOs = new ArrayList<ReconstructedParticle>();
+
+ public PFOJetFindingDriver() {
+ }
+
+ public void process(EventHeader event) {
+ PFOs = event.get(ReconstructedParticle.class, inputPfoCollection);
+
+ Map<HepLorentzVector, ReconstructedParticle> map = new HashMap<HepLorentzVector, ReconstructedParticle>();
+ for (ReconstructedParticle pfo : PFOs) {
+// if (pfo.getCharge() != 0) {
+ if (Math.abs(VecOp.cosTheta(pfo.getMomentum()))<0.95) {
+ map.put(pfo.asFourVector(), pfo);
+ }
+// }
+ }
+ JetFinder jetFind = new JadeEJetFinder(yCut);
+ List<List<ReconstructedParticle>> jetList = new ArrayList<List<ReconstructedParticle>>();
+ jetFind.setEvent(map.keySet());
+ for (int i = 0; i < jetFind.njets(); ++i) {
+ List<HepLorentzVector> jetVecList = jetFind.particlesInJet(i);
+ List<ReconstructedParticle> jetPartList = new ArrayList<ReconstructedParticle>();
+ for (HepLorentzVector iVec : jetVecList) {
+ jetPartList.add(map.get(iVec));
+ }
+ jetList.add(jetPartList);
+ }
+ event.put("jetList", jetList);
+ }
+}
\ No newline at end of file
diff -N AddBremPhotonsToElectrons.java.hn2 --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ AddBremPhotonsToElectrons.java.hn2 11 Oct 2012 16:17:22 -0000 1.1 @@ -0,0 +1,128 @@
+package org.lcsim.recon.postrecon.leptonID.electron;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Track;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import hep.physics.vec.*;
+import java.util.*;
+import java.io.IOException;
+public class AddBremPhotonsToElectrons extends Driver
+{
+ String rpl = "PandoraPFOCollection";
+ String outrpn = "eehReconstructedParticles";
+//
+// Mass squared cut for associating photon with electron
+//
+ double msqcut = 0.49;
+//
+// # sigma cut on electron cluster match to ignore photon
+// correction
+//
+ double nsigematch = 4.;
+//
+// # sigma cut on combined cluster energy to ignore photon
+// correction
+//
+ double nsigtotignore = -999.;
+ public AddBremPhotonsToElectrons()
+ {
+ }
+ public void setMassSquaredCut(double x){msqcut = x;}
+ public void setOutputRPName(String x){outrpn = x;}
+ protected void process(EventHeader event)
+ {
+ List<ReconstructedParticle> el = new ArrayList<ReconstructedParticle>();
+ List<ReconstructedParticle> phl = new ArrayList<ReconstructedParticle>();
+ List<ReconstructedParticle> nrp = new ArrayList<ReconstructedParticle>(event.get(ReconstructedParticle.class,rpl));
+ List<ReconstructedParticle> rp = event.get(ReconstructedParticle.class,rpl);
+
+//
+// Find the reconstructed electrons and photons
+//
+ for(ReconstructedParticle p:rp)
+ {
+ if(p.getCharge() == 0.)
+ {
+ if(p.getMass() == 0.)phl.add(p);
+ }
+ else
+ {
+ if(p.getMass() < .001)
+ {
+ el.add(p);
+ }
+ }
+ }
+//
+// For each electron, loop over photons and use inv mass to
+// identify brem photons
+//
+ for(ReconstructedParticle pp:el)
+ {
+ Hep3Vector ep = pp.getMomentum();
+ double eE = pp.getEnergy();
+ List<ReconstructedParticle> addphl = new ArrayList<ReconstructedParticle>();
+ double addE = 0.;
+ for(ReconstructedParticle p:phl)
+ {
+ Hep3Vector php = p.getMomentum();
+ double phE = p.getEnergy();
+ double masssq = (eE+phE)*(eE+phE) - (ep.x()+php.x())*(ep.x()+php.x()) -
+ (ep.y()+php.y())*(ep.y()+php.y()) - (ep.z()+php.z())*(ep.z()+php.z());
+ if(masssq < msqcut)
+ {
+ addphl.add(p);
+ addE += phE;
+ }
+ }
+ double addf = 1. + addE/eE;
+ if(addE > 0.)
+ {
+//
+// Check for electron ID without an e/p match. In that case, assume combined brem and electron
+// shower, and just use the Track measurement.
+//
+ double clE = 0.;
+ for(Cluster c:pp.getClusters())
+ {
+ clE += c.getEnergy();
+ }
+ if(clE/eE < (1. - nsigematch*.18/Math.sqrt(eE)))addf = 1.;
+//
+// Also check that adding the brem photons exceed e/p match. If not,
+// assume the photons are electron fragments.
+//
+ if( (clE+addE)/eE < (1. + nsigtotignore*.18/Math.sqrt(eE)) )addf = 1.;
+//
+// Now modify the output ReconstructedParticle list
+//
+ BaseReconstructedParticle newe = new BaseReconstructedParticle(addf*eE,addf*ep.x(),addf*ep.y(),addf*ep.z());
+ newe.setMass(pp.getMass());
+ newe.setCharge(pp.getCharge());
+ newe.setReferencePoint(pp.getReferencePoint());
+ newe.setParticleIdUsed(pp.getParticleIDUsed());
+ newe.addTrack(pp.getTracks().get(0));
+ nrp.remove(pp);
+ nrp.add(newe);
+ for(Cluster c:pp.getClusters())
+ {
+ newe.addCluster(c);
+ }
+ newe.addParticle(pp);
+ for(ReconstructedParticle ph:addphl)
+ {
+ nrp.remove(ph);
+ newe.addParticle(ph);
+ for(Cluster c:ph.getClusters())
+ {
+ newe.addCluster(c);
+ }
+ }
+ }
+ }
+ event.put(outrpn,nrp,ReconstructedParticle.class,0);
+ }
+}
diff -N IsolatedHighPElectronIdentifier.java.hn2 --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ IsolatedHighPElectronIdentifier.java.hn2 11 Oct 2012 16:17:22 -0000 1.1 @@ -0,0 +1,261 @@
+package org.lcsim.recon.postrecon.leptonID.electron;
+import org.lcsim.util.Driver;
+import hep.physics.particle.properties.ParticlePropertyManager;
+import hep.physics.particle.properties.ParticlePropertyProvider;
+import org.lcsim.event.ParticleID;
+import org.lcsim.event.base.BaseParticleID;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackState;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import org.lcsim.util.lcio.LCIOConstants;
+import org.lcsim.geometry.*;
+import hep.physics.vec.*;
+import java.util.*;
+import java.io.IOException;
+public class IsolatedHighPElectronIdentifier extends Driver
+{
+ String rpl = "ReconstructedParticles";
+ protected String inputPfoCollection = "PandoraPFOCollection";
+
+ double conect;
+ double maxE;
+ double minE;
+ final ParticlePropertyProvider dPPP;
+ ParticleID eppid;
+ ParticleID empid;
+ double emass;
+ Hep3Vector origin;
+ public IsolatedHighPElectronIdentifier()
+ {
+ conect = 0.90;
+ maxE = 125.;
+ minE = 15.;
+ dPPP = ParticlePropertyManager.getParticlePropertyProvider();
+ eppid = new BaseParticleID(dPPP.get(-11));
+ empid = new BaseParticleID(dPPP.get(11));
+ emass = dPPP.get(11).getMass();
+ origin = new BasicHep3Vector(0.,0.,0.);
+ }
+ public void setConeCosTheta(double x){conect = x;}
+ public void setMaxE(double x){maxE = x;}
+ public void setMinE(double x){minE = x;}
+
+
+ // List<ReconstructedParticle> rp = new ArrayList<ReconstructedParticle>();
+
+ protected void process(EventHeader event)
+ {
+//
+// Match Tracks with Recon Particles
+//
+ System.out.println("This is the one!");
+
+ // if (1==1) return;
+
+ List<Track> allTracks = event.get(Track.class,"Tracks");
+ List<Track> unusedTracks = new ArrayList<Track>(allTracks);
+ // List<ReconstructedParticle> rp = event.get(ReconstructedParticle.class,rpl);
+ List<ReconstructedParticle> rp = new ArrayList<ReconstructedParticle>();
+ rp = event.get(ReconstructedParticle.class, inputPfoCollection);
+
+ Map<Track,ReconstructedParticle> ttorp = new HashMap<Track,ReconstructedParticle>();
+ Map<Cluster,ReconstructedParticle> ctorp = new HashMap<Cluster,ReconstructedParticle>();
+ for(ReconstructedParticle p:rp)
+ {
+ if(p.getCharge() != 0.)
+ {
+ Track t = p.getTracks().get(0);
+ ttorp.put(t,p);
+ unusedTracks.remove(t);
+ }
+ else
+ {
+ if(p.getClusters().size() > 0)ctorp.put(p.getClusters().get(0),p);
+ }
+ }
+//
+// Map non-muon tracks > minE GeV to clusters within a 25 degree cone
+//
+ List<Track> xtrae = new ArrayList<Track>();
+ // List<Cluster> cl = event.get(Cluster.class,"Clusters");
+ List<Cluster> cl = event.get(Cluster.class,"ReconClusters");
+ Map<Track,List<Cluster>> ttocll = new HashMap<Track,List<Cluster>>();
+ for(Track trck:allTracks)
+ {
+ TrackState t = trck.getTrackStates().get(0);
+ boolean eid = false;
+ Hep3Vector tp = new BasicHep3Vector(t.getMomentum());
+ if(tp.magnitude() < minE)continue;
+ if(tp.magnitude() > maxE)continue;
+ if(ttorp.containsKey(trck))
+ {
+ double mass = ttorp.get(trck).getMass();
+ if( Math.abs(mass - .105) < .005)continue;
+ if( Math.abs(mass - .0005) < .0001)eid = true;
+ }
+ List<Cluster> tcl = new ArrayList<Cluster>();
+ double HcalE = 0.;
+ for(Cluster c:cl)
+ {
+ Hep3Vector cp = new BasicHep3Vector(c.getPosition());
+ double ct = (tp.x()*cp.x()+tp.y()*cp.y()+tp.z()*cp.z())/tp.magnitude()/cp.magnitude();
+ // double hec = 0.0;
+ // double hec = c.getSubdetectorEnergies()[3] + c.getSubdetectorEnergies()[7];
+ // System.out.println(" ....... energy = " + hec);
+ double hec = 0.0;
+
+ List<CalorimeterHit> calhitsList = c.getCalorimeterHits();
+ for (CalorimeterHit pfochit : calhitsList) {
+ IDDecoder _decoder = pfochit.getIDDecoder();
+ String dname = _decoder.getSubdetector().getName();
+ // System.out.println(" name of subdetector for this hit " + dname);
+ if (dname.equals("HcalBarrel") || dname.equals("HcalEndcap")) {
+ hec += pfochit.getCorrectedEnergy();
+ // System.out.println(" HCAL energy contribution = "+pfochit.getCorrectedEnergy());
+ }
+ }
+ System.out.println(" hec = "+hec);
+
+ if(ct > conect)
+ {
+ tcl.add(c);
+ HcalE += hec;
+ }
+ }
+ double f = HcalE/tp.magnitude();
+ if(!eid)
+ {
+ System.out.println(" f = "+f);
+ if(f < 1000.)xtrae.add(trck); // .04
+ }
+ ttocll.put(trck,tcl);
+ }
+ if(xtrae.size() < 1)return;
+//
+// Have Tracks IDed as electrons. If we already have a ReconstructedParticle
+// for this Track, just replace it with clone with electron ID. Otherwise, we
+// need to make a new ReconstructedParticle and do some association
+// to determine which neutrals to remove.
+//
+ List<ReconstructedParticle> addlist = new ArrayList<ReconstructedParticle>();
+ List<ReconstructedParticle> removelist = new ArrayList<ReconstructedParticle>();
+ for(Track t:xtrae)
+ {
+ if(ttorp.containsKey(t))
+ {
+ ReconstructedParticle erp = makeElectronRP(ttorp.get(t));
+ addlist.add(erp);
+ removelist.add(ttorp.get(t));
+ }
+ else
+ {
+ BaseReconstructedParticle erp = makeElectronBRP(t);
+ if(ttocll.get(t).size() < 1)
+ {
+ addlist.add(erp);
+ }
+ else
+ {
+ List<Cluster> ocl = new ArrayList<Cluster>();
+ List<Cluster> mcl = new ArrayList<Cluster>(ttocll.get(t));
+ Hep3Vector eP = erp.getMomentum();
+ double eE = erp.getEnergy();
+ double cE = 0.;
+ for(int i=0;i<ttocll.get(t).size();i++)
+ {
+ double maxct = -1.;
+ Cluster c = null;
+ for(int j=0;j<mcl.size();j++)
+ {
+ Cluster tc = mcl.get(j);
+ Hep3Vector cp = new BasicHep3Vector(tc.getPosition());
+ double ct = (eP.x()*cp.x() + eP.y()*cp.y() + eP.z()*cp.z())/eP.magnitude()/cp.magnitude();
+ if(ct > maxct)
+ {
+ maxct = ct;
+ c = tc;
+ }
+ }
+ ocl.add(i,c);
+ mcl.remove(c);
+ }
+ List<Cluster> acl = new ArrayList<Cluster>();
+ for(Cluster c:ocl)
+ {
+ if( (cE + c.getEnergy()) < (eE + 3.*.2*Math.sqrt(eE)) )
+ {
+ cE += c.getEnergy();
+ acl.add(c);
+ }
+ }
+ addlist.add(erp);
+ double Erem = 0.;
+ int nrem = 0;
+ for(Cluster c:acl)
+ {
+ if(ctorp.containsKey(c))
+ {
+ erp.addCluster(c);
+ removelist.add(ctorp.get(c));
+ nrem++;
+ Erem += ctorp.get(c).getEnergy();
+ }
+ }
+ }
+ }
+ }
+ System.out.println("Removing and Adding Reconstructed Particles");
+
+ for(ReconstructedParticle p:removelist)
+ {
+ System.out.println("removing particle; mass = " + p.getMass());
+ rp.remove(p);
+ }
+ for(ReconstructedParticle p:addlist)
+ {
+ System.out.println("adding particle; mass = " + p.getMass());
+ rp.add(p);
+ }
+ }
+ public ReconstructedParticle makeElectronRP(ReconstructedParticle pi)
+ {
+ Hep3Vector m = pi.getMomentum();
+ double E = Math.sqrt(emass*emass + m.magnitudeSquared());
+ BaseReconstructedParticle p = new BaseReconstructedParticle(E,m);
+ p.setMass(emass);
+ p.setCharge(pi.getCharge());
+ p.setReferencePoint(pi.getReferencePoint());
+ if(pi.getCharge() > 0.)p.setParticleIdUsed(eppid);
+ else p.setParticleIdUsed(empid);
+ for(Cluster c:pi.getClusters())
+ {
+ p.addCluster(c);
+ }
+ for(Track c:pi.getTracks())
+ {
+ p.addTrack(c);
+ }
+ return p;
+ }
+ public BaseReconstructedParticle makeElectronBRP(Track trck)
+ {
+ TrackState t = trck.getTrackStates().get(0);
+ Hep3Vector m = new BasicHep3Vector(t.getMomentum());
+ double E = Math.sqrt(emass*emass + m.magnitudeSquared());
+ BaseReconstructedParticle p = new BaseReconstructedParticle(E,m);
+ p.setMass(emass);
+ double tcharge = t.getOmega()/Math.abs(t.getOmega());
+ // p.setCharge((double) (t.getOmega()/Math.abs(t.getOmega)));
+ p.setCharge(tcharge);
+ p.setReferencePoint(origin);
+ if(tcharge > 0.)p.setParticleIdUsed(eppid);
+ else p.setParticleIdUsed(empid);
+ p.addTrack(trck);
+ return p;
+ }
+}
diff -N MyDBDvvHAnalysiseid.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MyDBDvvHAnalysiseid.java 11 Oct 2012 16:17:22 -0000 1.1 @@ -0,0 +1,433 @@
+//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 java.util.List;
+import org.lcsim.util.Driver;
+import org.lcsim.event.*;
+import org.lcsim.util.lcio.LCIOConstants;
+import org.lcsim.geometry.*;
+import org.lcsim.recon.util.McTruthLinker;
+
+import java.util.ArrayList;
+import java.util.Collections;
+
+import org.lcsim.spacegeom.CartesianPoint;
+import org.lcsim.spacegeom.CartesianVector;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.spacegeom.SpaceVector;
+
+
+
+public class MyDBDvvHAnalysiseid extends Driver
+{
+
+ // collections
+ protected String inputPfoCollection = "PandoraPFOCollection";
+ protected String mcParticlesSkimmedName = "MCParticle";
+ protected String recoMCTruthName = "RecoMCTruthLink";
+
+ boolean first = true;
+
+ IProfile1D effprof1;
+ IProfile1D fakeprof1;
+
+ protected double E_BEAM = 500.0;
+
+ // First get the reconstructed particles from the event
+ List<ReconstructedParticle> PFOs = new ArrayList<ReconstructedParticle>();
+
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ protected void process(EventHeader event)
+ {
+
+ // try {
+ PFOs = event.get(ReconstructedParticle.class, inputPfoCollection);
+ // } catch (IllegalArgumentException e) {
+ // 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();
+
+ // List<List<ReconstructedParticle>> jetList = (List<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);
+ return;
+ }
+ int nJets = jetList.size();
+
+ 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());
+ 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.getGeneratorStatus() == particle.FINAL_STATE) {
+ int pdgid = particle.getPDGID();
+ int apdgid = Math.abs(pdgid);
+ if (particle.getGeneratorStatus() == particle.FINAL_STATE && !((apdgid == 12 || apdgid == 14 || apdgid == 16) && (particle.getParents().get(0).getPDGID()==pdgid))) {
+ // 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("final state MC particle PT").fill(PT);
+ aida.cloud1D("final state MC particle P").fill(P);
+ aida.cloud1D("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);
+ // }
+ }
+ //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;
+ double PsumTrkZ = 0.0;
+ for (Track track : tracklist) {
+ // 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]);
+
+ 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);
+// System.out.println("PT = "+PT);
+// 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);
+ if (M2visTrk>0.) aida.cloud1D("Mvis from Tracks").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());
+ }
+// ---------------------PFO-----------------------------------------------------
+ double EsumPFO = 0.0;
+ double PsumPFOX = 0.0;
+ double PsumPFOY = 0.0;
+ double PsumPFOZ = 0.0;
+ double eTotalInput = 0.0;
+ int nelec = 0;
+ for (ReconstructedParticle pfo : PFOs) {
+ boolean passPfoSelection = true;
+
+ MCParticle mcpart = null;
+ int pfoPDGID = 0;
+ List<LCRelation> MCTs = new ArrayList<LCRelation>();
+ MCTs = event.get(LCRelation.class, recoMCTruthName);
+ for(LCRelation mct:MCTs) {
+ if( (ReconstructedParticle) (mct.getFrom()) == pfo ) {
+ mcpart = (MCParticle) (mct.getTo());
+ pfoPDGID = mcpart.getPDGID();
+ break;
+ }
+ }
+ if (mcpart != null) {
+ double[] momentum = (mcpart.getMomentum()).v();
+ double mPX = momentum[0];
+ double mPY = momentum[1];
+ double mPZ = momentum[2];
+ SpaceVector pP = new SpaceVector(pfo.getMomentum());
+ double pPX = pP.x();
+ double pPY = pP.y();
+ double pPZ = pP.z();
+ System.out.println(" mPX,Y,Z = "+mPX+","+mPY+","+mPZ);
+ System.out.println(" pPX,Y,Z = "+pPX+","+pPY+","+pPZ);
+ System.out.println(" -------------------------------");
+ } else {
+ System.out.println(" ----assoc. mcpart no found-----");
+ }
+
+ /* -------------- just testing for the electron identifier -------------
+ double hec = 0.0;
+ List<Cluster> pfoclusList = pfo.getClusters();
+ for (Cluster pfoclus : pfoclusList) {
+ List<CalorimeterHit> calhitsList = pfoclus.getCalorimeterHits();
+ for (CalorimeterHit pfochit : calhitsList) {
+ IDDecoder _decoder = pfochit.getIDDecoder();
+ String dname = _decoder.getSubdetector().getName();
+ System.out.println(" name of subdetector for this hit " + dname);
+ if (dname.equals("HcalBarrel") || dname.equals("HcalEndcap")) {
+ hec += pfochit.getCorrectedEnergy();
+ System.out.println(" HCAL energy contribution = "+pfochit.getCorrectedEnergy());
+ }
+ }
+ System.out.println(" hec = "+hec);
+ }
+ ---------------------------------------------------------------------*/
+
+ int type = pfo.getType();
+ SpaceVector momentum = new SpaceVector(pfo.getMomentum());
+ double pT = momentum.rxy();
+ double P = momentum.rxyz();
+ double PX = momentum.x();
+ double PY = momentum.y();
+ double PZ = momentum.z();
+ double cosTheta = momentum.cosTheta();
+ double energy = pfo.getEnergy();
+
+ int foundelec = 0;
+ if (pfo.getParticleIDUsed() != null) {
+ int pidu = (pfo.getParticleIDUsed()).getPDG();
+ if (Math.abs(pidu)==11) {
+ nelec++;
+ foundelec = 1;
+ }
+ }
+
+ double CorrecteFound = 0.0;
+ if (Math.abs(pfoPDGID) == 11 && foundelec != 0 ) {
+ CorrecteFound = 1.0;
+ }
+ // aida.cloud2D("electron finding efficiency vs. PFO momentum").fill(P,CorrecteFound);
+ effprof1 = aida.profile1D("electron finding efficiency vs. PFO momentum", 10, 0., 100.);
+ effprof1.fill(P,CorrecteFound);
+
+
+ double FakeeFound = 0.0;
+ if (Math.abs(pfoPDGID) != 11 && foundelec != 0 ) {
+ FakeeFound = 1.0;
+ }
+ // aida.cloud2D("fake electron finding fraction vs. PFO momentum").fill(P,FakeeFound);
+ fakeprof1 = aida.profile1D("fake electron finding fraction vs. PFO momentum", 10, 0., 100.);
+ fakeprof1.fill(P,FakeeFound);
+
+
+ // 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);
+ if (M2visPFO>0.) aida.cloud1D("Mvis from PFO").fill(Math.sqrt(M2visPFO));
+ aida.cloud1D("Number of electrons per event").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);
+ }
+ 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));
+ }
+ }
+ }
+ }
+
+
+// -------------JETS------------------------------------------------------
+ double EsumJETS = 0.0;
+ double PsumJETSX = 0.0;
+ double PsumJETSY = 0.0;
+ double PsumJETSZ = 0.0;
+ double eJETTotalInput = 0.0;
+ // int nelec = 0;
+ boolean passjetSelection = true;
+ for (ReconstructedParticle jet : jetList) {
+
+ int ntrksjet = (jet.getParticles()).size();
+ aida.cloud1D("#tracks in jet").fill(ntrksjet);
+
+ // int type = pfo.getType();
+ SpaceVector momentum = new SpaceVector(jet.getMomentum());
+ double pT = momentum.rxy();
+ double P = momentum.rxyz();
+ double PX = momentum.x();
+ double PY = momentum.y();
+ double PZ = momentum.z();
+ double cosTheta = momentum.cosTheta();
+ double energy = jet.getEnergy();
+
+ // if (pfo.getParticleIDUsed() != null) {
+ // int pidu = (pfo.getParticleIDUsed()).getPDG();
+ // if (Math.abs(pidu)==11) nelec++;
+ // }
+
+ // System.out.println("energy = "+energy + "type = " + type);
+ if (Math.abs(PZ/P)>0.75 || energy>(2.0*E_BEAM)) {
+ passjetSelection = false;
+ } else {
+ aida.cloud1D("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);
+ }
+ 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 ----------------------------------------
+
+// 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);
+ }
+ 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));
+ }
+ }
+ }
+ }
+
+
+
+
+ }
+}
+
+
diff -N MyDBDvvHeid.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MyDBDvvHeid.java 11 Oct 2012 16:17:22 -0000 1.1 @@ -0,0 +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 MyDBDvvHeid extends Driver {
+
+ String rplnew = "eehReconstructedParticles";
+ String rpl = "ReconstructedParticles";
+
+
+ private AIDA aida = AIDA.defaultInstance();
+ public MyDBDvvHeid(){
+//
+// 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 MyDBDvvHAnalysiseid());
+
+
+ };
+ protected String mcParticlesSkimmedName = "MCParticle";
+
+ public void process(EventHeader event) {
+ 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");
+
+ }
+}
\ No newline at end of file
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