lcsim/src/org/lcsim/contrib/uiowa
diff -N MassAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MassAnalysis.java 4 Apr 2007 18:16:58 -0000 1.1
@@ -0,0 +1,370 @@
+package org.lcsim.contrib.uiowa;
+
+// Overkill!
+import java.util.*;
+import hep.physics.vec.*;
+import hep.physics.jet.*;
+import hep.aida.*;
+
+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 MassAnalysis extends Driver
+{
+ // I think this is needed for Ron's analysis:
+ private AIDA aida = AIDA.defaultInstance();
+
+ 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_h1trueRes;
+ ICloud1D m_h1trueResAdHoc;
+ ICloud1D m_h1PassingCut;
+ ICloud1D m_h1trueResPassingCut;
+ ICloud1D m_h1trueResPassingCutAdHoc;
+ ICloud1D m_h1PassingCut2;
+ ICloud1D m_h1trueResPassingCut2;
+ ICloud1D m_h1trueResPassingCut2AdHoc;
+
+ public MassAnalysis() {
+ // PFA
+ m_pfa = new NonTrivialPFA(false);
+ 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");
+ jNonTrivialAdHoc.setOutputCollectionName("jetOutputAdHoc");
+ jNonTrivialAdHoc.setFinder(twojet);
+ add(jNonTrivialAdHoc);
+ 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_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_h1PassingCut = m_histoFactory.createCloud1D("eventMassInBarrel");
+ m_h1trueResPassingCut = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel");
+ m_h1trueResPassingCutAdHoc = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrelAdHoc");
+ m_h1PassingCut2 = m_histoFactory.createCloud1D("eventMassInBarrel2");
+ m_h1trueResPassingCut2 = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel2");
+ m_h1trueResPassingCut2AdHoc = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel2AdHoc");
+ } 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; }
+ }
+ }
+ }
+
+ // DEBUG: Cross-check event energy sum
+ if (false) {
+ List<ReconstructedParticle>recoParticles = event.get(ReconstructedParticle.class,"all particles (ron calib)");
+ List<ReconstructedParticle>recoParticlesAdHoc = event.get(ReconstructedParticle.class,"all particles");
+ List<ReconstructedParticle>perfParticles = event.get(ReconstructedParticle.class,"perfect particles");
+ double recoSum = 0.0;
+ double recoSumAdHoc = 0.0;
+ double perfSum = 0.0;
+ Hep3Vector perfSumVec = new BasicHep3Vector();
+ Hep3Vector recoSumVec = new BasicHep3Vector();
+ for (ReconstructedParticle part : recoParticles) {
+ recoSum += part.getEnergy();
+ recoSumVec = VecOp.add(recoSumVec, part.getMomentum());
+ // Debug
+ double partEnergy = part.getEnergy();
+ Hep3Vector partMomVec = part.getMomentum();
+ double partMomMag = partMomVec.magnitude();
+ double partMass = Math.sqrt(partEnergy*partEnergy - partMomMag*partMomMag);
+ System.out.println("DEBUG/RECO: Particle has Q="+part.getCharge()+", E="+partEnergy+", p="+partMomMag+" => m="+partMass+". SumE="+recoSum+", sumP="+recoSumVec.magnitude());
+
+ }
+ for (ReconstructedParticle part : recoParticlesAdHoc) { recoSumAdHoc += part.getEnergy(); }
+ for (ReconstructedParticle part : perfParticles) {
+ perfSum += part.getEnergy();
+ perfSumVec = VecOp.add(perfSumVec, part.getMomentum());
+ // Debug
+ double partEnergy = part.getEnergy();
+ Hep3Vector partMomVec = part.getMomentum();
+ double partMomMag = partMomVec.magnitude();
+ double partMass = Math.sqrt(partEnergy*partEnergy - partMomMag*partMomMag);
+ System.out.println("DEBUG/PERFECT: Particle has Q="+part.getCharge()+", E="+partEnergy+", p="+partMomMag+" => m="+partMass+". SumE="+perfSum+", sumP="+perfSumVec.magnitude());
+ }
+ }
+ // DEBUG: Check just photons:
+ if (false) {
+ List<ReconstructedParticle>recoParticles = event.get(ReconstructedParticle.class,"large photon particles");
+ for (ReconstructedParticle part : recoParticles) {
+ System.out.println("DEBUG/RECO: large photon particles: Particle has Q="+part.getCharge()+", E="+part.getEnergy()+", p="+part.getMomentum().magnitude());
+ }
+ recoParticles = event.get(ReconstructedParticle.class,"small photons (ron calib)");
+ for (ReconstructedParticle part : recoParticles) {
+ System.out.println("DEBUG/RECO: small photons: Particle has Q="+part.getCharge()+", E="+part.getEnergy()+", p="+part.getMomentum().magnitude());
+ }
+ recoParticles = event.get(ReconstructedParticle.class,"neutral hadron particles (ron calib)");
+ for (ReconstructedParticle part : recoParticles) {
+ System.out.println("DEBUG/RECO: neutral hadron particles (ron calib): Particle has Q="+part.getCharge()+", E="+part.getEnergy()+", p="+part.getMomentum().magnitude());
+ }
+ recoParticles = event.get(ReconstructedParticle.class,"charged hadron particles 2");
+ for (ReconstructedParticle part : recoParticles) {
+ System.out.println("DEBUG/RECO: charged particles: Particle has Q="+part.getCharge()+", E="+part.getEnergy()+", p="+part.getMomentum().magnitude());
+ }
+ }
+
+
+ List<ReconstructedParticle>jets = event.get(ReconstructedParticle.class,"jetOutput");
+ List<ReconstructedParticle>jetsAdHoc = event.get(ReconstructedParticle.class,"jetOutputAdHoc");
+ List<ReconstructedParticle>jetsPerfect = event.get(ReconstructedParticle.class,"jetOutputTrivial");
+ double esum = 0.;
+ double esumAdHoc = 0.;
+ double esumPerfect = 0.;
+ double pxsum = 0.;
+ double pxsumAdHoc = 0.;
+ double pxsumPerfect = 0.;
+ double pysum = 0.;
+ double pysumAdHoc = 0.;
+ double pysumPerfect = 0.;
+ double pzsum = 0.;
+ double pzsumAdHoc = 0.;
+ double pzsumPerfect = 0.;
+ double chargesum = 0.;
+ double chargesumAdHoc = 0.;
+ double chargesumPerfect = 0.;
+ double[] je = new double[jets.size()];
+ double[] jeAdHoc = new double[jets.size()];
+ double[] jePerfect = new double[jets.size()];
+ Hep3Vector[] jp = new Hep3Vector[jets.size()];
+ Hep3Vector[] jpAdHoc = new Hep3Vector[jets.size()];
+ Hep3Vector[] jpPerfect = new Hep3Vector[jets.size()];
+ int ind = 0;
+ int indAdHoc = 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: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 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 evtmassPerfect = Math.sqrt(esumPerfect*esumPerfect - psumPerfect*psumPerfect);
+ System.out.println("DEBUG: Energy sum = "+esum+" (reco-ron), "+esumAdHoc+" (reco-AH), "+esumPerfect+" (perf)");
+ System.out.println("DEBUG: Momentum sum = "+psum+" (reco-ron), "+psumAdHoc+" (reco-AH), "+psumPerfect+" (perf)");
+ System.out.println("DEBUG: => mass = "+evtmass+" (reco-ron), "+evtmassAdHoc+" (reco-AH), "+evtmassPerfect+" (perf)");
+ 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_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());
+ }
+
+ // 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());
+ }
+ }
+ }
+ }
+
+ 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_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();
+ }
+}
