3 added files
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon/analysis
diff -N AssociateReconJetsWithPerfectJets.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ AssociateReconJetsWithPerfectJets.java 28 Jun 2011 14:26:06 -0000 1.1
@@ -0,0 +1,147 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon.analysis;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import java.util.*;
+import hep.physics.vec.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class AssociateReconJetsWithPerfectJets extends Driver
+{
+ protected String jetListName;
+ List<ReconstructedParticle> pjlo;
+ List<ReconstructedParticle> reconjets;
+ double sigB = .005;
+ double sigEpercent = .075;
+ boolean[] arrInit;
+ int[][][] arr;
+ int maxj = 8;
+
+ public void setJetListName(String jetListName)
+ {
+ this.jetListName = jetListName;
+ }
+ public AssociateReconJetsWithPerfectJets()
+ {
+ arrInit = new boolean[maxj+1];
+ arr = new int[maxj+1][][];
+ }
+ public AssociateReconJetsWithPerfectJets(String jln)
+ {
+ jetListName = jln;
+ arrInit = new boolean[maxj+1];
+ arr = new int[maxj+1][][];
+ }
+ protected void process(EventHeader event)
+ {
+ pjlo = event.get(ReconstructedParticle.class,"PerfectJets");
+ reconjets = event.get(ReconstructedParticle.class,jetListName);
+ List<Double> chil = new ArrayList<Double>();
+ if((reconjets.size() != pjlo.size())||(pjlo.size() > maxj)||(pjlo.size() < 2))
+ {
+ chil.add(-1.);
+ event.put(jetListName+"AssChisq",chil);
+ return;
+ }
+ // Make the best association with the reconstructed jets
+ double[][] chisqcont = new double[pjlo.size()][pjlo.size()];
+ for(int i=0;i<pjlo.size();i++)
+ {
+ ReconstructedParticle pjet = pjlo.get(i);
+ for(int j=0;j<pjlo.size();j++)
+ {
+ ReconstructedParticle rjet = reconjets.get(j);
+ chisqcont[i][j] = getChisqCont(rjet,pjet);
+ }
+ }
+ int[][] assarray = getAssArray(pjlo.size());
+ int minio = -1;
+ double chisqmin = 9999999.;
+ for(int i=0;i<assarray.length;i++)
+ {
+ double chisq = 0.;
+ for(int j=0;j<pjlo.size();j++)
+ {
+ chisq += chisqcont[j][assarray[i][j]];
+ }
+ chisq /= pjlo.size();
+ if(chisq < chisqmin)
+ {
+ minio = i;
+ chisqmin = chisq;
+ }
+ }
+ chil.add(chisqmin);
+ event.put(jetListName+"AssChisq",chil);
+ List<ReconstructedParticle> treconjets = new ArrayList<ReconstructedParticle>(reconjets);
+ reconjets.removeAll(treconjets);
+ for(int i=0;i<pjlo.size();i++)
+ {
+ reconjets.add(treconjets.get(assarray[minio][i]));
+ }
+ }
+ public double getChisqCont(ReconstructedParticle p1,ReconstructedParticle p2)
+ {
+ double sigE = sigEpercent*p1.getEnergy();
+ double chisqc = 0.;
+ double E1 = p1.getEnergy();
+ double E2 = p2.getEnergy();
+ Hep3Vector P1 = p1.getMomentum();
+ Hep3Vector P2 = p2.getMomentum();
+ double dE = E1 - E2;
+ Hep3Vector dB = VecOp.sub(VecOp.mult(1./E1,P1),VecOp.mult(1./E2, P2));
+ chisqc += dE*dE/sigE;
+ chisqc += dB.x()*dB.x()/sigB;
+ chisqc += dB.y()*dB.y()/sigB;
+ chisqc += dB.z()*dB.z()/sigB;
+ chisqc /= 4.;
+ return chisqc;
+ }
+ public int[][] getAssArray(int n)
+ {
+ if(arrInit[n])return arr[n];
+ int ncomb = 1;
+ for(int i=n;i>1;i--)ncomb*=i;
+ arr[n] = new int[ncomb][n];
+ int repeat = ncomb;
+ for(int i=0;i<n;i++)
+ {
+ int ic = 0;
+ repeat /= (n-i);
+ int niter = ncomb/repeat/(n-i);
+ for(int iter=0;iter<niter;iter++)
+ {
+ for(int j=0;j<n;j++)
+ {
+ if(valid(arr[n],ic,i,j))
+ {
+ for(int k=0;k<repeat;k++)
+ {
+ arr[n][ic][i] = j;
+ ic++;
+ }
+ }
+ }
+ }
+ }
+ arrInit[n] = true;
+ return arr[n];
+ }
+ public boolean valid(int[][] a, int ic, int i, int j)
+ {
+ if(ic >= a.length)return false;
+ for(int k=0;k<i;k++)
+ {
+ if(a[ic][k] == j)return false;
+ }
+ return true;
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon/analysis
diff -N MakePerfectPDG94Jets.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MakePerfectPDG94Jets.java 28 Jun 2011 14:26:06 -0000 1.1
@@ -0,0 +1,168 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon.analysis;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import java.util.*;
+import hep.physics.vec.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class MakePerfectPDG94Jets extends Driver
+{
+ List<ReconstructedParticle> pjlo;
+ RPESort rpes;
+ public MakePerfectPDG94Jets()
+ {
+ rpes = new RPESort();
+ }
+ protected void process(EventHeader event)
+ {
+ List<MCParticle> all = event.get(MCParticle.class,"MCParticle");
+ List<MCParticle> fs = event.get(MCParticle.class,"GenFinalStateParticles");
+ // Ordered list of perfect jets
+ pjlo = new ArrayList<ReconstructedParticle>();
+ // Create the perfect pdg94 ReconstructedParticles
+ List<ReconstructedParticle> p94l = new ArrayList<ReconstructedParticle>();
+ for(MCParticle p:all)
+ {
+ if(p.getPDGID() == 94)
+ {
+ // Make a ReconstructedParticle from the pdg94
+ ReconstructedParticle r94 = makeRP94(p,fs);
+ // Add it to the list
+ p94l.add(r94);
+ // Sort the jets for this pdg94 by energy
+ Collections.sort(r94.getParticles(),rpes);
+ }
+ }
+ event.put("PerfectPDG94", p94l);
+ // Order the jets: jet 0 = maxE jet, jet1 = jet0 partner,
+ // jet2 = maxE jet from other pdg94, jet3 = jet2 partner
+ List<ReconstructedParticle> t94l = new ArrayList<ReconstructedParticle>(p94l);
+ while(t94l.size() > 0)
+ {
+ double maxE = 0.;
+ int n = t94l.size();
+ int i94 = -1;
+ for(int i=0;i<n;i++)
+ {
+ if(t94l.get(i).getParticles().get(0).getEnergy() > maxE)
+ {
+ maxE = t94l.get(i).getParticles().get(0).getEnergy();
+ i94 = i;
+ }
+ }
+ for(ReconstructedParticle jet:t94l.get(i94).getParticles())
+ {
+ pjlo.add(jet);
+ }
+ t94l.remove(i94);
+ }
+ event.put("PerfectJets",pjlo,ReconstructedParticle.class,0);
+ }
+ public ReconstructedParticle makeRP94(MCParticle pdg94,List<MCParticle> fsl)
+ {
+ // Make the original ReconstructedParticle
+ BaseReconstructedParticle pdg94rp = rpFromMC(pdg94);
+ // Construct and add the quark ReconstructedParticles
+ int type = 0;
+ int nused = 0;
+ for(MCParticle p:pdg94.getParents())
+ {
+ if(pdg94.getParents().size() > 2)
+ {
+ if(p.getDaughters().size() > 1)continue;
+ }
+ // Make the original quark reconstructedParticle
+ BaseReconstructedParticle rp = rpFromMC(p);
+ // Add it to the pdg94 particle
+ pdg94rp.addParticle(rp);
+ rp.setType(p.getPDGID());
+ type += p.getPDGID();
+ nused++;
+ }
+ // Adjust the pdg94 type
+ if(nused != 2)type -= 10;
+ else
+ {
+ if(pdg94.getParents().size() > 2)
+ {
+ type += 10;
+ }
+ }
+ pdg94rp.setType(type);
+ // Now add the fs particles traced to this pdg94 to the appropriate quark
+ for(MCParticle p:fsl)
+ {
+ if(trace(p,pdg94))
+ {
+ ReconstructedParticle rp = rpFromMC(p);
+ findNearestQ(rp,pdg94rp.getParticles(),pdg94rp).getParticles().add(rp);
+ }
+ }
+ // Now replace the quark 4-Vector with the actual sum of its components
+ for(ReconstructedParticle q:pdg94rp.getParticles())
+ {
+ HepLorentzVector v4 = new BasicHepLorentzVector();
+ for(ReconstructedParticle p:q.getParticles())
+ {
+ v4 = VecOp.add(v4,p.asFourVector());
+ }
+ BaseReconstructedParticle bq = (BaseReconstructedParticle) q;
+ bq.set4Vector(v4);
+ bq.setMass(v4.magnitude());
+ }
+ return pdg94rp;
+ }
+ public BaseReconstructedParticle rpFromMC(MCParticle p)
+ {
+ BaseReconstructedParticle rp = new BaseReconstructedParticle(p.getEnergy(),p.getMomentum());
+ rp.setMass(p.getMass());
+ rp.setType(p.getPDGID());
+ double charge = 0.;
+ try
+ {
+ charge = p.getCharge();
+ }
+ catch (hep.physics.particle.properties.UnknownParticleIDException e)
+ {
+ }
+ rp.setCharge(charge);
+ return rp;
+ }
+ public boolean trace(MCParticle fsp, MCParticle p94)
+ {
+ MCParticle pp = fsp;
+ while(pp.getParents().size() > 0)
+ {
+ if(pp.getParents().get(0) == p94)return true;
+ pp = pp.getParents().get(0);
+ }
+ return false;
+ }
+ public ReconstructedParticle findNearestQ(ReconstructedParticle p, List<ReconstructedParticle> ql, ReconstructedParticle bp)
+ {
+ ReconstructedParticle np = null;
+ double maxct = -2.;
+ HepLorentzVector pl = VecOp.boost(p.asFourVector(),bp.asFourVector());
+ for(ReconstructedParticle q:ql)
+ {
+ double ct = VecOp.dot(VecOp.boost(q.asFourVector(),bp.asFourVector()).v3(),pl.v3());
+ if(ct > maxct)
+ {
+ maxct = ct;
+ np = q;
+ }
+ }
+ return np;
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon/analysis
diff -N AssociateReconFullE4JetsWithPerfectJets.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ AssociateReconFullE4JetsWithPerfectJets.java 28 Jun 2011 14:26:06 -0000 1.1
@@ -0,0 +1,138 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon.analysis;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import java.util.*;
+import hep.physics.vec.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class AssociateReconFullE4JetsWithPerfectJets extends Driver
+{
+ protected String jetListName;
+ List<ReconstructedParticle> pjlo;
+ List<ReconstructedParticle> reconjets;
+ double sigB = .005;
+ double sigEpercent = .075;
+
+ public void setJetListName(String jetListName)
+ {
+ this.jetListName = jetListName;
+ }
+ public AssociateReconFullE4JetsWithPerfectJets()
+ {
+ }
+ public AssociateReconFullE4JetsWithPerfectJets(String jln)
+ {
+ jetListName = jln;
+ }
+ protected void process(EventHeader event)
+ {
+ pjlo = event.get(ReconstructedParticle.class,"PerfectJets");
+ reconjets = event.get(ReconstructedParticle.class,jetListName);
+ List<Double> chil = new ArrayList<Double>();
+ int[] order = {0,1,2,3};
+ if((reconjets.size() != 4)||(pjlo.size() != 4))
+ {
+ chil.add(-1.);
+ event.put(jetListName+"AssChisq",chil);
+ return;
+ }
+ // Make the best association with the reconstructed jets
+ // By construction, should already have jets paired
+ // So first associate dijets
+ HepLorentzVector[] rdijet = new HepLorentzVector[2];
+ HepLorentzVector[] pdijet = new HepLorentzVector[2];
+ rdijet[0] = new BasicHepLorentzVector(reconjets.get(0).getEnergy()+reconjets.get(1).getEnergy(),
+ VecOp.add(reconjets.get(0).getMomentum(),reconjets.get(1).getMomentum()));
+ rdijet[1] = new BasicHepLorentzVector(reconjets.get(2).getEnergy()+reconjets.get(3).getEnergy(),
+ VecOp.add(reconjets.get(2).getMomentum(),reconjets.get(3).getMomentum()));
+ pdijet[0] = new BasicHepLorentzVector(pjlo.get(0).getEnergy()+pjlo.get(1).getEnergy(),
+ VecOp.add(pjlo.get(0).getMomentum(),pjlo.get(1).getMomentum()));
+ pdijet[1] = new BasicHepLorentzVector(pjlo.get(2).getEnergy()+pjlo.get(3).getEnergy(),
+ VecOp.add(pjlo.get(2).getMomentum(),pjlo.get(3).getMomentum()));
+ boolean reversed = false;
+ if( (getChisqCont(rdijet[1],pdijet[0])+ getChisqCont(rdijet[0],pdijet[1])) <
+ (getChisqCont(rdijet[0],pdijet[0])+ getChisqCont(rdijet[1],pdijet[1])) )reversed = true;
+ if(reversed)
+ {
+ order[0] = 2;
+ order[1] = 3;
+ order[2] = 0;
+ order[3] = 1;
+ }
+ // Now get the pairings within the dijet
+ double chisq0 = getChisqCont(reconjets.get(order[0]),pjlo.get(0)) + getChisqCont(reconjets.get(order[1]),pjlo.get(1));
+ double chisq1 = getChisqCont(reconjets.get(order[1]),pjlo.get(0)) + getChisqCont(reconjets.get(order[0]),pjlo.get(1));
+ double chisq = chisq0;
+ if( chisq1 < chisq0)
+ {
+ int temp = order[0];
+ order[0] = order[1];
+ order[1] = temp;
+ chisq = chisq1;
+ }
+ double chis0 = getChisqCont(reconjets.get(order[2]),pjlo.get(2)) + getChisqCont(reconjets.get(order[3]),pjlo.get(3));
+ double chis1 = getChisqCont(reconjets.get(order[3]),pjlo.get(2)) + getChisqCont(reconjets.get(order[2]),pjlo.get(3));
+ double chis = chis0;
+ if(chis1 < chis0)
+ {
+ int temp = order[2];
+ order[2] = order[3];
+ order[3] = temp;
+ chis = chis1;
+ }
+ chisq += chis;
+ chisq /= 4.;
+ chil.add(chisq);
+ event.put(jetListName+"AssChisq",chil);
+ // Now rewrite the recon jets in the matched order
+ List<ReconstructedParticle> treconjets = new ArrayList<ReconstructedParticle>(reconjets);
+ reconjets.removeAll(treconjets);
+ for(int i=0;i<4;i++)
+ {
+ reconjets.add(treconjets.get(order[i]));
+ }
+ }
+ public double getChisqCont(ReconstructedParticle p1,ReconstructedParticle p2)
+ {
+ double sigE = sigEpercent*p1.getEnergy();
+ double chisqc = 0.;
+ double E1 = p1.getEnergy();
+ double E2 = p2.getEnergy();
+ Hep3Vector P1 = p1.getMomentum();
+ Hep3Vector P2 = p2.getMomentum();
+ double dE = E1 - E2;
+ Hep3Vector dB = VecOp.sub(VecOp.mult(1./E1,P1),VecOp.mult(1./E2, P2));
+ chisqc += dE*dE/sigE;
+ chisqc += dB.x()*dB.x()/sigB;
+ chisqc += dB.y()*dB.y()/sigB;
+ chisqc += dB.z()*dB.z()/sigB;
+ chisqc /= 4.;
+ return chisqc;
+ }
+ public double getChisqCont(HepLorentzVector v1,HepLorentzVector v2)
+ {
+ double sigE = sigEpercent*v1.t();
+ double chisqc = 0.;
+ double E1 = v1.t();
+ double E2 = v2.t();
+ Hep3Vector P1 = v1.v3();
+ Hep3Vector P2 = v2.v3();
+ double dE = E1 - E2;
+ Hep3Vector dB = VecOp.sub(VecOp.mult(1./E1,P1),VecOp.mult(1./E2, P2));
+ chisqc += dE*dE/sigE;
+ chisqc += dB.x()*dB.x()/sigB;
+ chisqc += dB.y()*dB.y()/sigB;
+ chisqc += dB.z()*dB.z()/sigB;
+ chisqc /= 4.;
+ return chisqc;
+ }
+}
CVSspam 0.2.8