lcsim/src/org/lcsim/contrib/uiowa
diff -N MassPlots.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MassPlots.java 27 Aug 2007 21:27:39 -0000 1.1
@@ -0,0 +1,217 @@
+package org.lcsim.contrib.uiowa;
+
+// Overkill!
+import java.util.*;
+import hep.physics.vec.*;
+import hep.physics.jet.*;
+import hep.aida.*;
+
+import org.lcsim.event.EventHeader.LCMetaData;
+
+import org.lcsim.util.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.util.hitmap.*;
+import org.lcsim.event.*;
+import org.lcsim.event.util.*;
+import org.lcsim.digisim.*;
+import org.lcsim.recon.cluster.mst.*;
+import org.lcsim.recon.cluster.mipfinder.*;
+import org.lcsim.recon.cluster.clumpfinder.*;
+import org.lcsim.recon.cluster.structural.*;
+import org.lcsim.recon.cluster.structural.likelihood.*;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.recon.cluster.analysis.*;
+import org.lcsim.mc.fast.tracking.*;
+import org.lcsim.recon.pfa.identifier.*;
+import org.lcsim.recon.pfa.output.*;
+import org.lcsim.recon.cluster.cheat.*;
+import org.lcsim.contrib.uiowa.NonTrivialPFA;
+import org.lcsim.contrib.uiowa.TrivialPFA;
+import java.io.IOException;
+
+import org.lcsim.util.aida.AIDA;
+
+/**
+ * Analysis: Compute the dijet invariant mass.
+ * Code mostly stolen from Ron.
+ * Assumes event topology is e+e- -> ZZ, Z_1 -> nu nubar, Z_2 -> q qbar
+ */
+
+public class MassPlots extends Driver
+{
+
+ ITree m_tree = null;
+ IHistogramFactory m_histoFactory = null;
+ ICloud1D m_h1;
+ ICloud1D m_h2;
+ ICloud1D m_h3;
+ ICloud1D m_h1trueRes;
+ ICloud1D m_h1PassingCut;
+ ICloud1D m_h1trueResPassingCut;
+ ICloud1D m_h1PassingCut2;
+ ICloud1D m_h1trueResPassingCut2;
+ String m_jetListName;
+
+ public MassPlots(String inputList, String outputFilename) {
+ m_jetListName = new String();
+ m_jetListName += "jetOutput__";
+ m_jetListName += inputList;
+ JetFinder twojet = new FixNumberOfJetsFinder(2);
+ JetDriver jNonTrivial = new JetDriver();
+ jNonTrivial.setInputCollectionName(inputList);
+ jNonTrivial.setOutputCollectionName(m_jetListName);
+ jNonTrivial.setFinder(twojet);
+ add(jNonTrivial);
+
+ m_eventCount = 0;
+
+ // Init plots
+ IAnalysisFactory af = IAnalysisFactory.create();
+ try {
+ m_tree = af.createTreeFactory().create(outputFilename, "xml", false, true);
+ m_histoFactory = af.createHistogramFactory(m_tree);
+ m_h1 = m_histoFactory.createCloud1D("eventMass");
+ m_h2 = m_histoFactory.createCloud1D("eventEnergy");
+ m_h3 = m_histoFactory.createCloud1D("eventMomentum");
+ m_h1trueRes = m_histoFactory.createCloud1D("eventMassResidualsToTruth");
+ m_h1PassingCut = m_histoFactory.createCloud1D("eventMassInBarrel");
+ m_h1trueResPassingCut = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel");
+ m_h1PassingCut2 = m_histoFactory.createCloud1D("eventMassInBarrel2");
+ m_h1trueResPassingCut2 = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel2");
+ } catch (IOException ioe1) {
+ ioe1.printStackTrace();
+ }
+ }
+
+ int m_eventCount;
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+
+ // Event cut: Both reconstructed jets have cos(theta) < 0.8
+ // Event cut 2: Both quarks have cos(theta) < 0.8
+ boolean passesTruthAcceptanceCut = false;
+ MCParticle partZTovv = null;
+ MCParticle partZToqq = null;
+ List<MCParticle> partZToqq_daughters = null;
+ {
+ List<MCParticle> mcParticles = event.getMCParticles();
+ for (MCParticle part : mcParticles) {
+ if (part.getPDGID() == 23) {
+ Vector<MCParticle> foundLightQuarkDaughters = new Vector<MCParticle>();
+ Vector<MCParticle> foundNeutrinoDaughters = new Vector<MCParticle>();
+ List<MCParticle> daughters = part.getDaughters();
+ for (MCParticle dau : daughters) {
+ int pdg = dau.getPDGID();
+ if (Math.abs(pdg)>0 && Math.abs(pdg)<4) {
+ foundLightQuarkDaughters.add(dau);
+ } else if (Math.abs(pdg)==12 || Math.abs(pdg)==14 || Math.abs(pdg)==16) {
+ foundNeutrinoDaughters.add(dau);
+ }
+ }
+ if (foundLightQuarkDaughters.size() == 2) {
+ partZToqq = part;
+ partZToqq_daughters = foundLightQuarkDaughters;
+ }
+ if (foundNeutrinoDaughters.size() == 2) {
+ partZTovv = part;
+ }
+ }
+ }
+ if (partZTovv == null) {
+ System.out.println("ERROR: no Z -> vv cand");
+ }
+ if (partZToqq == null) {
+ System.out.println("ERROR: no Z -> qq cand");
+ } else {
+ System.out.println("DEBUG: Mass of Z -> qq = "+partZToqq.getMass());
+ passesTruthAcceptanceCut = true;
+ for (MCParticle dau : partZToqq_daughters) {
+ double cosTheta = Math.abs(dau.getMomentum().z() / dau.getMomentum().magnitude());
+ if (cosTheta>0.8) { passesTruthAcceptanceCut = false; }
+ }
+ }
+ }
+
+ List<ReconstructedParticle>jets = event.get(ReconstructedParticle.class, m_jetListName);
+ double esum = 0.;
+ double pxsum = 0.;
+ double pysum = 0.;
+ double pzsum = 0.;
+ double chargesum = 0.;
+ double[] je = new double[jets.size()];
+ Hep3Vector[] jp = new Hep3Vector[jets.size()];
+ int ind = 0;
+
+ // Loop over jets
+ for(ReconstructedParticle jet:jets) {
+ esum += jet.getEnergy();
+ je[ind] = jet.getEnergy();
+ Hep3Vector v = jet.getMomentum();
+ jp[ind] = v;
+ pxsum += v.x();
+ pysum += v.y();
+ pzsum += v.z();
+ chargesum += jet.getCharge();
+ ind++;
+ }
+
+ // Compute total momentum and hence event mass
+ double psum = Math.sqrt(pxsum*pxsum+pysum*pysum+pzsum*pzsum);
+ double evtmass = Math.sqrt(esum*esum - psum*psum);
+ if (partZToqq != null) {
+ System.out.println(" => Event residuals: mass-truth="+(evtmass-partZToqq.getMass()));
+ }
+
+ m_h1.fill(evtmass);
+ m_h2.fill(esum);
+ m_h3.fill(psum);
+
+ if (partZToqq != null) {
+ m_h1trueRes.fill(evtmass-partZToqq.getMass());
+ }
+
+ // Look at individual jets
+ for(int i=0;i<jets.size();i++) {
+ double ct = jp[i].z()/Math.sqrt(jp[i].x()*jp[i].x() + jp[i].y()*jp[i].y() + jp[i].z()*jp[i].z());
+ double jm = Math.sqrt(jp[i].x()*jp[i].x() + jp[i].y()*jp[i].y() + jp[i].z()*jp[i].z());
+ for(int j=i+1;j<jets.size();j++) {
+ double ct1 = jp[j].z()/Math.sqrt(jp[j].x()*jp[j].x() + jp[j].y()*jp[j].y() + jp[j].z()*jp[j].z());
+ double masssq = (je[i]+je[j])*(je[i]+je[j]) - (jp[i].x()+jp[j].x())*(jp[i].x()+jp[j].x()) - (jp[i].y()+jp[j].y())*(jp[i].y()+jp[j].y()) - (jp[i].z()+jp[j].z())*(jp[i].z()+jp[j].z());
+ double sign = 1.;
+ if(masssq < 0.)sign = -1.;
+ double signedMass = sign*Math.sqrt(sign*masssq);
+ if( (Math.abs(ct) < .8)&&(Math.abs(ct1) < .8) ) {
+ System.out.println("EVENT PASSES CUT, since cos(theta) = "+ct+", "+ct1+". Jet-jet mass = "+signedMass);
+ m_h1PassingCut.fill(evtmass);
+ if (partZToqq != null) {
+ m_h1trueResPassingCut.fill(evtmass-partZToqq.getMass());
+ }
+ }
+ }
+ }
+
+ if (passesTruthAcceptanceCut) {
+ System.out.println("EVENT PASSES CUT2. Jet-jet mass = "+evtmass);
+ m_h1PassingCut2.fill(evtmass);
+ m_h1trueResPassingCut2.fill(evtmass-partZToqq.getMass());
+ }
+
+
+ m_eventCount++;
+ if (m_eventCount % 50 == 0) {
+ // Checkpoint
+ System.out.println("DEBUG: Checkpoint at "+m_eventCount);
+ try { m_tree.commit(); } catch(IOException ioe1) { ioe1.printStackTrace(); }
+ }
+ }
+
+ public void suspend() {
+ try {
+ m_tree.commit();
+ } catch(IOException ioe1) {
+ ioe1.printStackTrace();
+ }
+ super.suspend();
+ }
+}
lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.4 -r1.5
--- MassAnalysis.java 23 Apr 2007 21:46:50 -0000 1.4
+++ MassAnalysis.java 27 Aug 2007 21:27:39 -0000 1.5
@@ -2,40 +2,17 @@
// Overkill!
import java.util.*;
-import hep.physics.vec.*;
-import hep.physics.jet.*;
import hep.aida.*;
-import org.lcsim.event.EventHeader.LCMetaData;
-
import org.lcsim.util.*;
-import org.lcsim.util.decision.*;
-import org.lcsim.util.hitmap.*;
import org.lcsim.event.*;
import org.lcsim.event.util.*;
-import org.lcsim.digisim.*;
-import org.lcsim.recon.cluster.mst.*;
-import org.lcsim.recon.cluster.mipfinder.*;
-import org.lcsim.recon.cluster.clumpfinder.*;
-import org.lcsim.recon.cluster.structural.*;
-import org.lcsim.recon.cluster.structural.likelihood.*;
-import org.lcsim.recon.cluster.util.*;
-import org.lcsim.recon.cluster.analysis.*;
-import org.lcsim.mc.fast.tracking.*;
-import org.lcsim.recon.pfa.identifier.*;
-import org.lcsim.recon.pfa.output.*;
-import org.lcsim.recon.cluster.cheat.*;
import org.lcsim.contrib.uiowa.NonTrivialPFA;
import org.lcsim.contrib.uiowa.TrivialPFA;
-import java.io.IOException;
+import org.lcsim.contrib.uiowa.MassPlots;
-import org.lcsim.util.aida.AIDA;
-/**
- * Analysis: Compute the dijet invariant mass.
- * Code mostly stolen from Ron.
- * Assumes event topology is e+e- -> ZZ, Z_1 -> nu nubar, Z_2 -> q qbar
- */
+import org.lcsim.util.aida.AIDA;
public class MassAnalysis extends Driver
{
@@ -44,34 +21,6 @@
NonTrivialPFA m_pfa = null;
TrivialPFA m_pfa2 = null;
- ITree m_tree = null;
- IHistogramFactory m_histoFactory = null;
- ICloud1D m_h1;
- ICloud1D m_h2;
- ICloud1D m_h3;
- ICloud1D m_h1res;
- ICloud1D m_h2res;
- ICloud1D m_h3res;
- ICloud1D m_h1perf;
- ICloud1D m_h2perf;
- ICloud1D m_h3perf;
- ICloud1D m_h1resAdHoc;
- ICloud1D m_h2resAdHoc;
- ICloud1D m_h3resAdHoc;
- ICloud1D m_h1resNoMerge;
- ICloud1D m_h2resNoMerge;
- ICloud1D m_h3resNoMerge;
- ICloud1D m_h1trueRes;
- ICloud1D m_h1trueResAdHoc;
- ICloud1D m_h1trueResNoMerge;
- ICloud1D m_h1PassingCut;
- ICloud1D m_h1trueResPassingCut;
- ICloud1D m_h1trueResPassingCutAdHoc;
- ICloud1D m_h1trueResPassingCutNoMerge;
- ICloud1D m_h1PassingCut2;
- ICloud1D m_h1trueResPassingCut2;
- ICloud1D m_h1trueResPassingCut2AdHoc;
- ICloud1D m_h1trueResPassingCut2NoMerge;
public MassAnalysis() {
// PFA
@@ -79,285 +28,7 @@
add(m_pfa);
m_pfa2 = new TrivialPFA(true);
add(m_pfa2);
- // Jets
- JetFinder twojet = new FixNumberOfJetsFinder(2);
- JetDriver jNonTrivial = new JetDriver();
- jNonTrivial.setInputCollectionName("all particles (ron calib)");
- jNonTrivial.setOutputCollectionName("jetOutput");
- jNonTrivial.setFinder(twojet);
- add(jNonTrivial);
- JetDriver jNonTrivialAdHoc = new JetDriver();
- jNonTrivialAdHoc.setInputCollectionName("all particles (adhoc calib)");
- jNonTrivialAdHoc.setOutputCollectionName("jetOutputAdHoc");
- jNonTrivialAdHoc.setFinder(twojet);
- add(jNonTrivialAdHoc);
- JetDriver jNonTrivialNoMerge = new JetDriver();
- jNonTrivialNoMerge.setInputCollectionName("all particles [no frag merge] (ron calib)");
- jNonTrivialNoMerge.setOutputCollectionName("jetOutputNoMerge");
- jNonTrivialNoMerge.setFinder(twojet);
- add(jNonTrivialNoMerge);
- JetDriver jTrivial = new JetDriver();
- jTrivial.setInputCollectionName("perfect particles");
- jTrivial.setOutputCollectionName("jetOutputTrivial");
- jTrivial.setFinder(twojet);
- add(jTrivial);
-
- m_eventCount = 0;
-
- // Init plots
- IAnalysisFactory af = IAnalysisFactory.create();
- try {
- m_tree = af.createTreeFactory().create("output.aida","xml",false,true);
- m_histoFactory = af.createHistogramFactory(m_tree);
- m_h1 = m_histoFactory.createCloud1D("eventMass");
- m_h2 = m_histoFactory.createCloud1D("eventEnergy");
- m_h3 = m_histoFactory.createCloud1D("eventMomentum");
- m_h1res = m_histoFactory.createCloud1D("eventMassResiduals");
- m_h2res = m_histoFactory.createCloud1D("eventEnergyResiduals");
- m_h3res = m_histoFactory.createCloud1D("eventMomentumResiduals");
- m_h1resAdHoc = m_histoFactory.createCloud1D("eventMassResidualsAdHoc");
- m_h2resAdHoc = m_histoFactory.createCloud1D("eventEnergyResidualsAdHoc");
- m_h3resAdHoc = m_histoFactory.createCloud1D("eventMomentumResidualsAdHoc");
- m_h1resNoMerge = m_histoFactory.createCloud1D("eventMassResidualsNoMerge");
- m_h2resNoMerge = m_histoFactory.createCloud1D("eventEnergyResidualsNoMerge");
- m_h3resNoMerge = m_histoFactory.createCloud1D("eventMomentumResidualsNoMerge");
- m_h1perf = m_histoFactory.createCloud1D("eventMassPerfect");
- m_h2perf = m_histoFactory.createCloud1D("eventEnergyPerfect");
- m_h3perf = m_histoFactory.createCloud1D("eventMomentumPerfect");
- m_h1trueRes = m_histoFactory.createCloud1D("eventMassResidualsToTruth");
- m_h1trueResAdHoc = m_histoFactory.createCloud1D("eventMassResidualsToTruthAdHoc");
- m_h1trueResNoMerge = m_histoFactory.createCloud1D("eventMassResidualsToTruthNoMerge");
- m_h1PassingCut = m_histoFactory.createCloud1D("eventMassInBarrel");
- m_h1trueResPassingCut = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel");
- m_h1trueResPassingCutAdHoc = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrelAdHoc");
- m_h1trueResPassingCutNoMerge = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrelNoMerge");
- m_h1PassingCut2 = m_histoFactory.createCloud1D("eventMassInBarrel2");
- m_h1trueResPassingCut2 = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel2");
- m_h1trueResPassingCut2AdHoc = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel2AdHoc");
- m_h1trueResPassingCut2NoMerge = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel2NoMerge");
- } catch (IOException ioe1) {
- ioe1.printStackTrace();
- }
- }
-
- int m_eventCount;
- protected void process(EventHeader event)
- {
- super.process(event);
-
- // Event cut: Both reconstructed jets have cos(theta) < 0.8
- // Event cut 2: Both quarks have cos(theta) < 0.8
- boolean passesTruthAcceptanceCut = false;
-
- // DEBUG: Print out MC lists
- MCParticle partZTovv = null;
- MCParticle partZToqq = null;
- List<MCParticle> partZToqq_daughters = null;
- {
- List<MCParticle> mcParticles = event.getMCParticles();
- for (MCParticle part : mcParticles) {
- if (part.getPDGID() == 23) {
- Vector<MCParticle> foundLightQuarkDaughters = new Vector<MCParticle>();
- Vector<MCParticle> foundNeutrinoDaughters = new Vector<MCParticle>();
- List<MCParticle> daughters = part.getDaughters();
- for (MCParticle dau : daughters) {
- int pdg = dau.getPDGID();
- if (Math.abs(pdg)>0 && Math.abs(pdg)<4) {
- foundLightQuarkDaughters.add(dau);
- } else if (Math.abs(pdg)==12 || Math.abs(pdg)==14 || Math.abs(pdg)==16) {
- foundNeutrinoDaughters.add(dau);
- }
- }
- if (foundLightQuarkDaughters.size() == 2) {
- partZToqq = part;
- partZToqq_daughters = foundLightQuarkDaughters;
- }
- if (foundNeutrinoDaughters.size() == 2) {
- partZTovv = part;
- }
- }
- }
- if (partZTovv == null) {
- System.out.println("ERROR: no Z -> vv cand");
- }
- if (partZToqq == null) {
- System.out.println("ERROR: no Z -> qq cand");
- } else {
- System.out.println("DEBUG: Mass of Z -> qq = "+partZToqq.getMass());
- passesTruthAcceptanceCut = true;
- for (MCParticle dau : partZToqq_daughters) {
- double cosTheta = Math.abs(dau.getMomentum().z() / dau.getMomentum().magnitude());
- if (cosTheta>0.8) { passesTruthAcceptanceCut = false; }
- }
- }
- }
-
- List<ReconstructedParticle>jets = event.get(ReconstructedParticle.class,"jetOutput");
- List<ReconstructedParticle>jetsAdHoc = event.get(ReconstructedParticle.class,"jetOutputAdHoc");
- List<ReconstructedParticle>jetsNoMerge = event.get(ReconstructedParticle.class,"jetOutputNoMerge");
- List<ReconstructedParticle>jetsPerfect = event.get(ReconstructedParticle.class,"jetOutputTrivial");
- double esum = 0.;
- double esumAdHoc = 0.;
- double esumNoMerge = 0.;
- double esumPerfect = 0.;
- double pxsum = 0.;
- double pxsumAdHoc = 0.;
- double pxsumNoMerge = 0.;
- double pxsumPerfect = 0.;
- double pysum = 0.;
- double pysumAdHoc = 0.;
- double pysumNoMerge = 0.;
- double pysumPerfect = 0.;
- double pzsum = 0.;
- double pzsumAdHoc = 0.;
- double pzsumNoMerge = 0.;
- double pzsumPerfect = 0.;
- double chargesum = 0.;
- double chargesumAdHoc = 0.;
- double chargesumNoMerge = 0.;
- double chargesumPerfect = 0.;
- double[] je = new double[jets.size()];
- double[] jeAdHoc = new double[jets.size()];
- double[] jeNoMerge = new double[jets.size()];
- double[] jePerfect = new double[jets.size()];
- Hep3Vector[] jp = new Hep3Vector[jets.size()];
- Hep3Vector[] jpAdHoc = new Hep3Vector[jets.size()];
- Hep3Vector[] jpNoMerge = new Hep3Vector[jets.size()];
- Hep3Vector[] jpPerfect = new Hep3Vector[jets.size()];
- int ind = 0;
- int indAdHoc = 0;
- int indNoMerge = 0;
- int indPerfect = 0;
-
- // Loop over jets
- for(ReconstructedParticle jet:jets) {
- esum += jet.getEnergy();
- je[ind] = jet.getEnergy();
- Hep3Vector v = jet.getMomentum();
- jp[ind] = v;
- pxsum += v.x();
- pysum += v.y();
- pzsum += v.z();
- chargesum += jet.getCharge();
- ind++;
- }
- for(ReconstructedParticle jet:jetsAdHoc) {
- esumAdHoc += jet.getEnergy();
- jeAdHoc[indAdHoc] = jet.getEnergy();
- Hep3Vector v = jet.getMomentum();
- jpAdHoc[indAdHoc] = v;
- pxsumAdHoc += v.x();
- pysumAdHoc += v.y();
- pzsumAdHoc += v.z();
- chargesumAdHoc += jet.getCharge();
- indAdHoc++;
- }
- for(ReconstructedParticle jet:jetsNoMerge) {
- esumNoMerge += jet.getEnergy();
- jeNoMerge[indNoMerge] = jet.getEnergy();
- Hep3Vector v = jet.getMomentum();
- jpNoMerge[indNoMerge] = v;
- pxsumNoMerge += v.x();
- pysumNoMerge += v.y();
- pzsumNoMerge += v.z();
- chargesumNoMerge += jet.getCharge();
- indNoMerge++;
- }
- for(ReconstructedParticle jet:jetsPerfect) {
- esumPerfect += jet.getEnergy();
- jePerfect[indPerfect] = jet.getEnergy();
- Hep3Vector v = jet.getMomentum();
- jpPerfect[indPerfect] = v;
- pxsumPerfect += v.x();
- pysumPerfect += v.y();
- pzsumPerfect += v.z();
- chargesumPerfect += jet.getCharge();
- indPerfect++;
- }
-
- // Compute total momentum and hence event mass
- double psum = Math.sqrt(pxsum*pxsum+pysum*pysum+pzsum*pzsum);
- double psumAdHoc = Math.sqrt(pxsumAdHoc*pxsumAdHoc+pysumAdHoc*pysumAdHoc+pzsumAdHoc*pzsumAdHoc);
- double psumNoMerge = Math.sqrt(pxsumNoMerge*pxsumNoMerge+pysumNoMerge*pysumNoMerge+pzsumNoMerge*pzsumNoMerge);
- double psumPerfect = Math.sqrt(pxsumPerfect*pxsumPerfect+pysumPerfect*pysumPerfect+pzsumPerfect*pzsumPerfect);
- double evtmass = Math.sqrt(esum*esum - psum*psum);
- double evtmassAdHoc = Math.sqrt(esumAdHoc*esumAdHoc - psumAdHoc*psumAdHoc);
- double evtmassNoMerge = Math.sqrt(esumNoMerge*esumNoMerge - psumNoMerge*psumNoMerge);
- double evtmassPerfect = Math.sqrt(esumPerfect*esumPerfect - psumPerfect*psumPerfect);
- System.out.println("DEBUG: Energy sum = "+esum+" (reco-ron), "+esumAdHoc+" (reco-AH), "+esumPerfect+" (perf), "+esumNoMerge+" (reco-ron,nomerge)");
- System.out.println("DEBUG: Momentum sum = "+psum+" (reco-ron), "+psumAdHoc+" (reco-AH), "+psumPerfect+" (perf), "+psumNoMerge+" (reco-ron,nomerge)");
- System.out.println("DEBUG: => mass = "+evtmass+" (reco-ron), "+evtmassAdHoc+" (reco-AH), "+evtmassPerfect+" (perf), "+evtmassNoMerge+" (reco-ron,nomerge)");
- if (partZToqq != null) {
- System.out.println(" => Event residuals: mass = "+(evtmass-evtmassPerfect)+", event energy = "+(esum-esumPerfect)+", mom = "+(psum-psumPerfect)+", mass-truth="+(evtmass-partZToqq.getMass()));
- } else {
- System.out.println(" => Event residuals: mass = "+(evtmass-evtmassPerfect)+", event energy = "+(esum-esumPerfect)+", mom = "+(psum-psumPerfect));
- }
-
- m_h1.fill(evtmass);
- m_h2.fill(esum);
- m_h3.fill(psum);
- m_h1res.fill(evtmass-evtmassPerfect);
- m_h2res.fill(esum-esumPerfect);
- m_h3res.fill(psum-psumPerfect);
- m_h1resAdHoc.fill(evtmassAdHoc-evtmassPerfect);
- m_h2resAdHoc.fill(esumAdHoc-esumPerfect);
- m_h3resAdHoc.fill(psumAdHoc-psumPerfect);
- m_h1resNoMerge.fill(evtmassNoMerge-evtmassPerfect);
- m_h2resNoMerge.fill(esumNoMerge-esumPerfect);
- m_h3resNoMerge.fill(psumNoMerge-psumPerfect);
- m_h1perf.fill(evtmassPerfect);
- m_h2perf.fill(esumPerfect);
- m_h3perf.fill(psumPerfect);
- if (partZToqq != null) {
- m_h1trueRes.fill(evtmass-partZToqq.getMass());
- m_h1trueResAdHoc.fill(evtmassAdHoc-partZToqq.getMass());
- m_h1trueResNoMerge.fill(evtmassNoMerge-partZToqq.getMass());
- }
-
- // Look at individual jets
- for(int i=0;i<jets.size();i++) {
- double ct = jp[i].z()/Math.sqrt(jp[i].x()*jp[i].x() + jp[i].y()*jp[i].y() + jp[i].z()*jp[i].z());
- double jm = Math.sqrt(jp[i].x()*jp[i].x() + jp[i].y()*jp[i].y() + jp[i].z()*jp[i].z());
- for(int j=i+1;j<jets.size();j++) {
- double ct1 = jp[j].z()/Math.sqrt(jp[j].x()*jp[j].x() + jp[j].y()*jp[j].y() + jp[j].z()*jp[j].z());
- double masssq = (je[i]+je[j])*(je[i]+je[j]) - (jp[i].x()+jp[j].x())*(jp[i].x()+jp[j].x()) - (jp[i].y()+jp[j].y())*(jp[i].y()+jp[j].y()) - (jp[i].z()+jp[j].z())*(jp[i].z()+jp[j].z());
- double sign = 1.;
- if(masssq < 0.)sign = -1.;
- double signedMass = sign*Math.sqrt(sign*masssq);
- if( (Math.abs(ct) < .8)&&(Math.abs(ct1) < .8) ) {
- System.out.println("EVENT PASSES CUT, since cos(theta) = "+ct+", "+ct1+". Jet-jet mass = "+signedMass);
- m_h1PassingCut.fill(evtmass);
- if (partZToqq != null) {
- m_h1trueResPassingCut.fill(evtmass-partZToqq.getMass());
- m_h1trueResPassingCutAdHoc.fill(evtmassAdHoc-partZToqq.getMass());
- m_h1trueResPassingCutNoMerge.fill(evtmassNoMerge-partZToqq.getMass());
- }
- }
- }
- }
-
- if (passesTruthAcceptanceCut) {
- System.out.println("EVENT PASSES CUT2. Jet-jet mass = "+evtmass);
- m_h1PassingCut2.fill(evtmass);
- m_h1trueResPassingCut2.fill(evtmass-partZToqq.getMass());
- m_h1trueResPassingCut2AdHoc.fill(evtmassAdHoc-partZToqq.getMass());
- m_h1trueResPassingCut2NoMerge.fill(evtmassNoMerge-partZToqq.getMass());
- }
-
- m_eventCount++;
- if (m_eventCount % 50 == 0) {
- // Checkpoint
- System.out.println("DEBUG: Checkpoint at "+m_eventCount);
- try { m_tree.commit(); } catch(IOException ioe1) { ioe1.printStackTrace(); }
- }
- }
-
- public void suspend() {
- try {
- m_tree.commit();
- } catch(IOException ioe1) {
- ioe1.printStackTrace();
- }
- super.suspend();
+ // Plots
+ add(new MassPlots("PFAReconstructedParticles", "massAnalysis.aida"));
}
}