5 added files
lcsim-contrib/src/main/java/org/lcsim/contrib/LGilbert
diff -N ClusterToReconDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ClusterToReconDriver.java 17 Aug 2009 16:43:42 -0000 1.1
@@ -0,0 +1,57 @@
+package org.lcsim.contrib.LGilbert;
+import java.util.*;
+import org.lcsim.util.*;
+import org.lcsim.event.util.*;
+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 org.lcsim.event.base.BaseParticleID;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.particle.properties.ParticlePropertyManager;
+import hep.physics.particle.properties.ParticlePropertyProvider;
+import java.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.util.*;
+import org.lcsim.event.*;
+import org.lcsim.event.base.*;
+import org.lcsim.recon.cluster.directedtree.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.recon.pfa.identifier.*;
+import org.lcsim.recon.cluster.util.CalorimeterHitTimeCutDecision;
+import org.lcsim.recon.cluster.util.UpperSubLayerDecision;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.vec.*;
+import Jama.Matrix;
+import org.lcsim.digisim.*;
+import hep.aida.*;
+
+public class ClusterToReconDriver extends Driver
+{
+ String inname;
+ String outname;
+ final ParticlePropertyProvider dPPP;
+ MakeReconstructedParticleFromCluster mr;
+
+ public ClusterToReconDriver(String in, String out)
+ {
+ inname = in;
+ outname = out;
+ dPPP = ParticlePropertyManager.getParticlePropertyProvider();
+ }
+
+ protected void process(EventHeader event)
+ {
+ mr = new MakeReconstructedParticleFromCluster();
+ List<Cluster> inlist = event.get(Cluster.class,inname);
+ List<ReconstructedParticle> outlist = new ArrayList<ReconstructedParticle>();
+ for(Cluster p:inlist)
+ {
+ outlist.add(mr.makeZeroMassRPFromCluster(p));
+ }
+ event.put(outname,outlist);
+ }
+}
\ No newline at end of file
lcsim-contrib/src/main/java/org/lcsim/contrib/LGilbert
diff -N MakeClustersForCalOnlyJets.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MakeClustersForCalOnlyJets.java 17 Aug 2009 16:43:42 -0000 1.1
@@ -0,0 +1,155 @@
+package org.lcsim.contrib.LGilbert;
+import java.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.util.*;
+import org.lcsim.event.*;
+import org.lcsim.event.base.*;
+import org.lcsim.recon.cluster.directedtree.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.recon.pfa.identifier.*;
+import org.lcsim.recon.cluster.util.CalorimeterHitTimeCutDecision;
+import org.lcsim.recon.cluster.util.UpperSubLayerDecision;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.vec.*;
+import Jama.Matrix;
+import org.lcsim.digisim.*;
+import hep.aida.*;
+
+/**
+ * Use DT clusters in the Ecal and Hcal > 10 hits to define particle objects,
+ * then add remaining hits to the nearest cluster
+ * Ron Cassell
+ */
+public class MakeClustersForCalOnlyJets extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ String outname;
+ String[] cn = {"EcalBarrDigiHitsDTreeClusters","EcalEndcapDigiHitsDTreeClusters",
+ "HcalBarrDigiHitsDTreeClusters","HcalEndcapDigiHitsDTreeClusters"};
+ String[] hn = {"CorrMuonBarrDigiHits","CorrMuonEndcapDigiHits"};
+ double alpha = -.23;
+ public MakeClustersForCalOnlyJets()
+ {
+ this("CalOnlyClusters");
+ }
+ public MakeClustersForCalOnlyJets(String out)
+ {
+ outname = out;
+//
+// Always use digisim
+//
+ add(new DigiPackageDriver());
+ DecisionMakerSingle<CalorimeterHit> upperLayer = new UpperSubLayerDecision();
+ DecisionMakerSingle<CalorimeterHit> lowerLayer = new NotDecisionMakerSingle(upperLayer);
+ DecisionMakerSingle<CalorimeterHit> timeCut = new CalorimeterHitTimeCutDecision(100);
+ AndDecisionMakerSingle<CalorimeterHit> lowerLayerAndTimeCut = new AndDecisionMakerSingle<CalorimeterHit>();
+ lowerLayerAndTimeCut.addDecisionMaker(lowerLayer);
+ lowerLayerAndTimeCut.addDecisionMaker(timeCut);
+ add(new ListFilterDriver(lowerLayerAndTimeCut, "MuonBarrHits", "CorrMuonBarrDigiHits", CalorimeterHit.class));
+ add(new ListFilterDriver(lowerLayerAndTimeCut, "MuonEndcapHits", "CorrMuonEndcapDigiHits", CalorimeterHit.class));
+//
+// Cluster with the DT clusterer
+//
+ add(new DirectedTreeDriver());
+ }
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+ List<Cluster> coc = new ArrayList<Cluster>();
+ List<Cluster> small = new ArrayList<Cluster>();
+ for(int i=0;i<cn.length;i++)
+ {
+ List<Cluster> cl = event.get(Cluster.class,cn[i]);
+ for(Cluster c:cl)
+ {
+ double ac = 1.;
+ if(i == 2)
+ {
+ double[] pos = c.getPosition();
+ double ct = Math.abs(pos[2])/Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
+ ac = 1./(1. + alpha*(1./Math.sqrt(1.- ct*ct) - 1));
+ }
+ if(i == 3)
+ {
+ double[] pos = c.getPosition();
+ double ct = Math.abs(pos[2])/Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
+ ac = 1./(1. + alpha*(1./ct - 1));
+ }
+ BasicCluster cc = (BasicCluster) c;
+ cc.setEnergy(ac*c.getEnergy());
+ if(c.getCalorimeterHits().size() < 10)
+ {
+ small.add(c);
+ }
+ else
+ {
+ coc.add(c);
+ }
+ }
+ }
+ double[] Ec = new double[coc.size()];
+ for(int i=0;i<coc.size();i++)
+ {
+ Ec[i] = coc.get(i).getEnergy();
+ }
+ for(int i=0;i<small.size();i++)
+ {
+ Hep3Vector p = new BasicHep3Vector(small.get(i).getPosition());
+ double ct = -1.;
+ int ic = -1;
+ for(int j=0;j<coc.size();j++)
+ {
+ Hep3Vector pc = new BasicHep3Vector(coc.get(j).getPosition());
+ double ctc = (p.x()*pc.x()+p.y()*pc.y()+p.z()*pc.z())/p.magnitude()/pc.magnitude();
+ if(ctc > ct)
+ {
+ ct = ctc;
+ ic = j;
+ }
+ }
+ Ec[ic] += small.get(i).getEnergy();
+ }
+ for(int j=0;j<hn.length;j++)
+ {
+ List<CalorimeterHit> hc = event.get(CalorimeterHit.class,hn[j]);
+ for(CalorimeterHit h:hc)
+ {
+ double ac = 1.;
+ if(j == 0)
+ {
+ double[] pos = h.getPosition();
+ double ct = Math.abs(pos[2])/Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
+ ac = 1./(1. + alpha*(1./Math.sqrt(1.- ct*ct) - 1));
+ }
+ else
+ {
+ double[] pos = h.getPosition();
+ double ct = Math.abs(pos[2])/Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
+ ac = 1./(1. + alpha*(1./ct - 1));
+ }
+ Hep3Vector p = new BasicHep3Vector(h.getPosition());
+ double ct = -1.;
+ int ic = -1;
+ for(int i=0;i<coc.size();i++)
+ {
+ Hep3Vector pc = new BasicHep3Vector(coc.get(i).getPosition());
+ double ctc = (p.x()*pc.x()+p.y()*pc.y()+p.z()*pc.z())/p.magnitude()/pc.magnitude();
+ if(ctc > ct)
+ {
+ ct = ctc;
+ ic = i;
+ }
+ }
+ Ec[ic] += ac*h.getCorrectedEnergy();
+ }
+ }
+ for(int i=0;i<coc.size();i++)
+ {
+ BasicCluster bc = (BasicCluster) coc.get(i);
+ bc.setEnergy(Ec[i]);
+ }
+ event.put(outname,coc);
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/LGilbert
diff -N MakePerfectReconParticles.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MakePerfectReconParticles.java 17 Aug 2009 16:43:42 -0000 1.1
@@ -0,0 +1,47 @@
+package org.lcsim.contrib.LGilbert;
+import java.util.*;
+import org.lcsim.util.*;
+import org.lcsim.event.util.*;
+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 org.lcsim.event.base.BaseParticleID;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.particle.properties.ParticlePropertyManager;
+import hep.physics.particle.properties.ParticlePropertyProvider;
+public class MakePerfectReconParticles extends Driver
+{
+ String inname;
+ String outname;
+ final ParticlePropertyProvider dPPP;
+ public MakePerfectReconParticles(String in, String out)
+ {
+ inname = in;
+ outname = out;
+ dPPP = ParticlePropertyManager.getParticlePropertyProvider();
+ }
+ protected void process(EventHeader event)
+ {
+ List<MCParticle> inlist = event.get(MCParticle.class,inname);
+ List<ReconstructedParticle> outlist = new ArrayList<ReconstructedParticle>();
+ for(MCParticle p:inlist)
+ {
+ if (!ParticleTypeClassifier.isNeutrino(p.getPDGID()))
+ {
+ outlist.add(makeRPFromMC(p));
+ }
+ }
+ event.put(outname,outlist);
+ }
+ public ReconstructedParticle makeRPFromMC(MCParticle p)
+ {
+ BaseReconstructedParticle rp = new BaseReconstructedParticle(p.getEnergy(),p.getMomentum());
+ rp.setMass(p.getMass());
+ rp.setCharge(p.getCharge());
+ rp.setReferencePoint(p.getOrigin());
+ rp.setParticleIdUsed(new BaseParticleID(dPPP.get(p.getType().getPDGID())));
+ return rp;
+ }
+}
\ No newline at end of file
lcsim-contrib/src/main/java/org/lcsim/contrib/LGilbert
diff -N MakeReconstructedParticleFromCluster.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MakeReconstructedParticleFromCluster.java 17 Aug 2009 16:43:42 -0000 1.1
@@ -0,0 +1,39 @@
+package org.lcsim.contrib.LGilbert;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import org.lcsim.event.Cluster;
+import hep.physics.vec.*;
+import hep.physics.particle.properties.ParticlePropertyManager;
+import hep.physics.particle.properties.ParticlePropertyProvider;
+import org.lcsim.event.ParticleID;
+import org.lcsim.event.base.BaseParticleID;
+
+/**
+ * Make 0 mass reconstructed particle from a cluster
+ */
+public class MakeReconstructedParticleFromCluster
+{
+ final ParticlePropertyProvider dPPP;
+ Hep3Vector origin;
+ ParticleID phpid;
+ public MakeReconstructedParticleFromCluster()
+ {
+ dPPP = ParticlePropertyManager.getParticlePropertyProvider();
+ origin = new BasicHep3Vector(0.,0.,0.);
+ phpid = new BaseParticleID(dPPP.get(22));
+ }
+ public ReconstructedParticle makeZeroMassRPFromCluster(Cluster c)
+ {
+ double E = c.getEnergy();
+ Hep3Vector x = new BasicHep3Vector(c.getPosition());
+ double r = E/x.magnitude();
+ Hep3Vector m = new BasicHep3Vector(x.x()*r,x.y()*r,x.z()*r);
+ BaseReconstructedParticle p = new BaseReconstructedParticle(E,m);
+ p.setMass(0.);
+ p.setCharge(0.);
+ p.setReferencePoint(origin);
+ p.setParticleIdUsed(phpid);
+ p.addCluster(c);
+ return p;
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/LGilbert
diff -N Reconstructed.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ Reconstructed.java 17 Aug 2009 16:43:42 -0000 1.1
@@ -0,0 +1,228 @@
+package org.lcsim.contrib.LGilbert;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.*;
+
+import hep.physics.jet.*;
+import hep.physics.vec.*;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import org.lcsim.util.*;
+import org.lcsim.event.util.*;
+import org.lcsim.util.aida.*;
+import hep.aida.*;
+
+public class Reconstructed extends Driver {
+ private double yCut = 0.001;
+ private AIDA aida = AIDA.defaultInstance();
+ String FS = "GenFSParticles";
+ String[] input = {"Durham/Real Reconstructed Particles/", "Durham/Perfect Reconstructed Particles/", "Durham/Cluster/", "Fix/Real Reconstructed Particles/", "Fix/Perfect Reconstructed Particles/", "Fix/Cluster/"};
+ String[] printfriendly = {"Real Reconstructed Particles w/ Durham Jet Finder", "Perfect Reconstructed Particles w/ Durham Jet Finder", "Clustering w/ Durham Jet Finder", "Real Reconstructed Particles w/ Fixed Number of Jets", "Perfect Reconstructed Particles w/ Fixed Number of Jets", "Clustering w/ Fixed Number of Jets"};
+ String inputcollectionname[] = {"ReconstructedParticles", "ReconParticleList", "ClusterReconParticles"};
+ String outputcollectionname[] = {"DurhamReconstructed", "DurhamPerfectRecon", "DurhamCluster", "FixReconstructed", "FixPerfect", "FixCluster"};
+ double rms90;
+ double mean90;
+ List <ICloud1D> mass;
+ List <ICloud1D> e;
+ MakeReconstructedParticleFromCluster mr;
+ boolean first;
+
+ public Reconstructed() {
+ rms90 = 0;
+ mean90 = 0;
+ CreateFinalStateMCParticleList cmc = new CreateFinalStateMCParticleList("Gen");
+ cmc.setCollectionName(FS);
+ add(cmc);
+
+ mass = new ArrayList<ICloud1D>();
+ e = new ArrayList<ICloud1D>();
+
+ first = true;
+
+ add(new MakePerfectReconParticles(FS, "ReconParticleList"));
+ add(new MakeClustersForCalOnlyJets("ClusterList"));
+ add(new ClusterToReconDriver("ClusterList", "ClusterReconParticles"));
+
+ JetDriver Durham1 = new JetDriver();
+ Durham1.setInputCollectionName(inputcollectionname[0]);
+ Durham1.setOutputCollectionName(outputcollectionname[0]);
+ Durham1.setFinder(new DurhamJetFinder(0.0075));
+ add(Durham1);
+
+ JetDriver Durham2 = new JetDriver();
+ Durham2.setInputCollectionName(inputcollectionname[1]);
+ Durham2.setOutputCollectionName(outputcollectionname[1]);
+ Durham2.setFinder(new DurhamJetFinder(0.0075));
+ add(Durham2);
+
+ JetDriver Durham3 = new JetDriver();
+ Durham3.setInputCollectionName(inputcollectionname[2]);
+ Durham3.setOutputCollectionName(outputcollectionname[2]);
+ Durham3.setFinder(new DurhamJetFinder(0.0075));
+ add(Durham3);
+
+ JetDriver Fix1 = new JetDriver();
+ Fix1.setInputCollectionName(inputcollectionname[0]);
+ Fix1.setOutputCollectionName(outputcollectionname[3]);
+ Fix1.setFinder(new FixNumberOfJetsFinder(2));
+ add(Fix1);
+
+ JetDriver Fix2 = new JetDriver();
+ Fix2.setInputCollectionName(inputcollectionname[1]);
+ Fix2.setOutputCollectionName(outputcollectionname[4]);
+ Fix2.setFinder(new FixNumberOfJetsFinder(2));
+ add(Fix2);
+
+ JetDriver Fix3 = new JetDriver();
+ Fix3.setInputCollectionName(inputcollectionname[2]);
+ Fix3.setOutputCollectionName(outputcollectionname[5]);
+ Fix3.setFinder(new FixNumberOfJetsFinder(2));
+ add(Fix3);
+ }
+
+ public double calculaterms90(ICloud1D cloud)
+ {
+ int entries = cloud.entries();
+ List<Double> elist = new ArrayList<Double>();
+
+ for(int i=0;i<entries;i++)
+ {
+ elist.add(cloud.value(i));
+ }
+
+ rms90 = 999999.;
+ int ntail = (int) (.1*entries);
+ int ncore = entries - ntail;
+ Collections.sort(elist);
+
+ for(int k=0;k<ntail;k++)
+ {
+ double sv = 0.;
+ double sv2 = 0.;
+
+ for(int i=k;i<k+ncore;i++)
+ {
+ sv += elist.get(i);
+ sv2 += elist.get(i)*elist.get(i);
+ }
+
+ double mean = sv/ncore;
+ double rms = Math.sqrt(sv2/ncore - mean*mean);
+
+ if(rms < rms90)
+ {
+ rms90 = rms;
+ mean90 = mean;
+ }
+ }
+ return rms90;
+ }
+
+ protected void suspend()
+ {
+ for (int i=0; i<mass.size(); i++)
+ {
+ System.out.println("\n" + printfriendly[i]);
+ System.out.println("RMS 90 for Mass: " + calculaterms90(mass.get(i)));
+ System.out.println("Mean 90 for Mass: " + mean90);
+ System.out.println("RMS 90 for Energy: " + calculaterms90(e.get(i)));
+ System.out.println("Mean 90 for Energy: " + mean90);
+ }
+
+ System.out.println("Number of mass entries: " + mass.size());
+ System.out.println("Number of energy entries: " + e.size());
+ }
+
+ public double getmean90(){return mean90;}
+
+ public void process(EventHeader event)
+ {
+ List<MCParticle> mcl = event.get(MCParticle.class,"MCParticle");
+ double Zmass = 0.;
+ double ZE = 0.;
+ double maxCosTheta = 0.;
+ //
+ // Loop over all MCParticles
+ //
+ for(MCParticle p:mcl)
+ {
+ int id = Math.abs(p.getPDGID());
+ //
+ // Look for light quarks
+ //
+ if( (id == 1 )||(id == 2 )||(id == 3 ) )
+ {
+ //
+ // If the light quark parent is a Z, use it
+ //
+ if(p.getParents().get(0).getPDGID() == 23 )
+ {
+ Zmass = p.getParents().get(0).getMass();
+ ZE = p.getParents().get(0).getEnergy();
+ double CosTheta = (Math.abs(p.getMomentum().z())/p.getMomentum().magnitude());
+ if (CosTheta > maxCosTheta)
+ {
+ maxCosTheta = CosTheta;
+ }
+ }
+ }
+ }
+ if (maxCosTheta > 0.9)
+ {
+ return;
+ }
+
+ super.process(event);
+
+ aida.cloud1D("Zmass").fill(Zmass);
+
+ if (first)
+ {
+ first = false;
+
+ for (int j = 0; j <6; j++)
+ {
+ mass.add(aida.cloud1D(input[j] + "difference in mass"));
+ e.add(aida.cloud1D(input[j] + "difference in energy"));
+ }
+ }
+
+ for (int lcv = 0; lcv < 6; lcv++)
+ {
+ List<ReconstructedParticle> jets = event.get(ReconstructedParticle.class, outputcollectionname[lcv]);
+
+ double px = 0;
+ double py = 0;
+ double pz = 0;
+ double energy = 0;
+ double totalmass = 0;
+ double difference = 0;
+ double deltae = 0;
+
+ for (ReconstructedParticle jet : jets)
+ {
+ Hep3Vector v = jet.getMomentum();
+ px += v.x();
+ py += v.y();
+ pz += v.z();
+ energy += jet.getEnergy();
+ }
+
+ totalmass = Math.sqrt(energy * energy - px * px - py * py - pz * pz);
+ aida.cloud1D(input[lcv] + "Total event mass").fill(totalmass);
+ aida.cloud1D(input[lcv] + "Total event energy").fill(energy);
+ aida.histogram1D(input[lcv] + "Number of jets", 10, -0.5, 9.5).fill(jets.size());
+
+ difference = totalmass - Zmass;
+ deltae = energy - ZE;
+
+ mass.get(lcv).fill(difference);
+ e.get(lcv).fill(deltae);
+ }
+ }
+}
CVSspam 0.2.8
