7 added files
lcsim/src/org/lcsim/contrib/Cassell/recon/analysis
diff -N CalClusterComposition.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CalClusterComposition.java 13 Sep 2007 15:10:11 -0000 1.1
@@ -0,0 +1,176 @@
+/*
+ * CalClusterComposition.java
+ *
+ * Created on September 8, 2007, 8:30 AM
+ *
+ * Class to decompose a cluster into contributions from
+ * photons, neutral hadrons, charged hadrons and electrons.
+ * Driven by a list of final state particles.
+ * Ron Cassell
+ */
+
+package org.lcsim.contrib.Cassell.recon.analysis;
+import org.lcsim.contrib.Cassell.recon.Cheat.TraceOrigin;
+import java.util.*;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimCalorimeterHit;
+import java.util.Map.Entry;
+
+/**
+ *
+ * @author cassell
+ */
+public class CalClusterComposition
+{
+ Cluster c;
+ double[] pfrac;
+ double[] pE;
+ int[] np;
+ MCParticle[] maxp;
+ double[] maxpE;
+ MCParticle mxp;
+ int mxptype;
+ double mxpE;
+ TraceOrigin tr;
+ int ntypes;
+ ParticleType pt;
+ Map<MCParticle,Double>[] Emap;
+
+ /** Creates a new instance of CalClusterComposition */
+ public CalClusterComposition(List<MCParticle> fsp)
+ {
+ tr = new TraceOrigin(fsp);
+ pt = new ParticleType();
+ ntypes = pt.getTypes().length;
+ pfrac = new double[ntypes];
+ np = new int[ntypes];
+ maxp = new MCParticle[ntypes];
+ maxpE = new double[ntypes];
+ pE = new double[ntypes];
+ }
+ public void decomposeCluster(Cluster cc)
+ {
+ c = cc;
+//
+// Map MCParticle to energy contribution for each type of particle separately
+//
+ Emap = new Map[ntypes];
+ for(int i=0;i<ntypes;i++)
+ {
+ pE[i] = 0.;
+ Emap[i] = new HashMap<MCParticle,Double>();
+ }
+//
+// Loop over the hits
+//
+ for(CalorimeterHit h:c.getCalorimeterHits())
+ {
+ SimCalorimeterHit hh = (SimCalorimeterHit) h;
+//
+// If only 1 contributor to the hit, no need to decompose the hit
+//
+ if(hh.getMCParticleCount() == 1)
+ {
+ MCParticle pp = tr.traceit(hh.getMCParticle(0));
+ if(pp == null)
+ {
+ pp = hh.getMCParticle(0);
+ }
+ int itype = pt.getType(pp);
+ double hitE = hh.getCorrectedEnergy();
+ pE[itype] += hitE;
+ if(Emap[itype].containsKey(pp))
+ {
+ double e = hitE + Emap[itype].get(pp).doubleValue();
+ Emap[itype].remove(pp);
+ Emap[itype].put(pp,new Double(e));
+ }
+ else
+ {
+ Emap[itype].put(pp,new Double(hitE));
+ }
+ }
+//
+// More than 1 contributor to the hit: decompose it
+//
+ else
+ {
+//
+// Only raw energy contributions stored, so corrected energy contribution is
+// is (fraction of raw energy)*(hit corrected energy)
+//
+ double totraw = 0.;
+ double[] cont = new double[hh.getMCParticleCount()];
+ MCParticle[] mccont = new MCParticle[hh.getMCParticleCount()];
+ for(int i=0;i<hh.getMCParticleCount();i++)
+ {
+ MCParticle pp = tr.traceit(hh.getMCParticle(i));
+ if(pp == null)pp = hh.getMCParticle(i);
+ mccont[i] = pp;
+ totraw += hh.getContributedEnergy(i);
+ cont[i] = hh.getContributedEnergy(i);
+ }
+ for(int i=0;i<hh.getMCParticleCount();i++)
+ {
+ MCParticle pp = mccont[i];
+ int itype = pt.getType(mccont[i]);
+ double hitE = hh.getCorrectedEnergy()*cont[i]/totraw;
+ pE[itype] += hitE;
+ if(Emap[itype].containsKey(pp))
+ {
+ double e = hitE + Emap[itype].get(pp).doubleValue();
+ Emap[itype].remove(pp);
+ Emap[itype].put(pp,new Double(e));
+ }
+ else
+ {
+ Emap[itype].put(pp,new Double(hitE));
+ }
+ }
+ }
+ }
+//
+// Have energy contributions, so calculate fractions, max contributor, etc.
+//
+ mxpE = 0.;
+ mxp = null;
+ mxptype = -1;
+ for(int i=0;i<ntypes;i++)
+ {
+ pfrac[i] = pE[i]/c.getEnergy();
+ np[i] = 0;
+ maxp[i] = null;
+ maxpE[i] = 0.;
+ for (Entry<MCParticle,Double> entry : Emap[i].entrySet())
+ {
+ MCParticle p = entry.getKey();
+ double e = entry.getValue().doubleValue();
+ if(e/c.getEnergy() > .1)np[i]++;
+ if(e > maxpE[i])
+ {
+ maxpE[i] = e;
+ maxp[i] = p;
+ }
+ if(e > mxpE)
+ {
+ mxpE = e;
+ mxp = p;
+ mxptype = i;
+ }
+ }
+ }
+ }
+ public Map<MCParticle,Double>[] getEmap(){return Emap;}
+ public Cluster getCluster(){return c;}
+ public MCParticle getMaxParticle(){return mxp;}
+ public double getMaxParticleE(){return mxpE;}
+ public int getMaxParticleType(){return mxptype;}
+ public double[] getEnergies(){return pE;}
+ public double[] getFractions(){return pfrac;}
+ public double[] getMaxEnergies(){return maxpE;}
+ public MCParticle[] getMaxParticles(){return maxp;}
+ public int[] getNParticles(){return np;}
+
+}
lcsim/src/org/lcsim/contrib/Cassell/recon/analysis
diff -N CalorimeterHitCollectionEnergies.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CalorimeterHitCollectionEnergies.java 13 Sep 2007 15:10:11 -0000 1.1
@@ -0,0 +1,124 @@
+/*
+ * CalorimeterHitCollectionEnergies.java
+ *
+ * Created on September 8, 2007, 5:26 PM
+ *
+ * Driver to store and optionally plot energies from a set of CalorimeterHit
+ * collectiosn as a function of particle type
+ */
+
+package org.lcsim.contrib.Cassell.recon.analysis;
+import java.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.contrib.Cassell.recon.Cheat.TraceOrigin;
+
+/**
+ *
+ * @author cassell
+ */
+public class CalorimeterHitCollectionEnergies extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ String[] coll;
+ String pre;
+ String FSname;
+ String[] type;
+ int nmax = 10000000;
+ TraceOrigin to;
+ double[] totE;
+ double[] runE;
+ boolean makeHists;
+ ParticleType pt;
+
+ /** Creates a new instance of CalorimeterHitCollectionEnergies */
+ public CalorimeterHitCollectionEnergies(String Prefix,String fsname,String[] collnames,boolean plot)
+ {
+ coll = collnames;
+ pre = Prefix+"/";
+ FSname = fsname;
+ makeHists = plot;
+ pt = new ParticleType();
+ type = pt.getTypes();
+ runE = new double[type.length];
+ }
+ public CalorimeterHitCollectionEnergies(String Prefix,String fsname,String[] collnames)
+ {
+ this(Prefix,fsname,collnames,false);
+ }
+ protected void process(EventHeader event)
+ {
+ List<MCParticle> fslist = event.get(MCParticle.class,FSname);
+ to = new TraceOrigin(fslist);
+ List<List<CalorimeterHit>> all = event.get(CalorimeterHit.class);
+ totE = new double[type.length];
+ for(List<CalorimeterHit> hl:all)
+ {
+ boolean useit = false;
+ for(int i=0;i<coll.length;i++)
+ {
+ if( (event.getMetaData(hl).getName().compareTo(coll[i]) == 0) )
+ {
+ useit = true;
+ break;
+ }
+ }
+ if(useit)
+ {
+ for(CalorimeterHit h:hl)
+ {
+ SimCalorimeterHit hh = (SimCalorimeterHit) h;
+ if(hh.getMCParticleCount() == 1)
+ {
+ MCParticle pp = to.traceit(hh.getMCParticle(0));
+ if(pp == null)
+ {
+ pp = hh.getMCParticle(0);
+ }
+ int itype = pt.getType(pp);
+ totE[itype] += h.getCorrectedEnergy();
+ runE[itype] += h.getCorrectedEnergy();
+ }
+ else
+ {
+ double totraw = 0.;
+ double[] cont = new double[type.length];
+ double fac = hh.getCorrectedEnergy()/hh.getRawEnergy();
+ for(int ii=0;ii<hh.getMCParticleCount();ii++)
+ {
+ MCParticle pp = to.traceit(hh.getMCParticle(ii));
+ if(pp == null)
+ {
+ pp = hh.getMCParticle(ii);
+ }
+ int itype = pt.getType(pp);
+ totraw += hh.getContributedEnergy(ii);
+ cont[itype] += hh.getContributedEnergy(ii);
+ }
+ for(int ii=0;ii<type.length;ii++)
+ {
+ totE[ii] += cont[ii]/totraw*h.getCorrectedEnergy();
+ runE[ii] += cont[ii]/totraw*h.getCorrectedEnergy();
+ }
+ }
+ }
+ }
+ }
+ if(makeHists)
+ {
+ for(int i=0;i<type.length;i++)
+ {
+ aida.cloud1D(pre+type[i]+" cal energy per event",nmax).fill(totE[i]);
+ }
+ }
+ }
+ public double[] getEnergies()
+ {return totE;}
+ public double[] getRunEnergies()
+ {return runE;}
+
+}
lcsim/src/org/lcsim/contrib/Cassell/recon/analysis
diff -N ClusterCollectionEnergies.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ClusterCollectionEnergies.java 13 Sep 2007 15:10:11 -0000 1.1
@@ -0,0 +1,100 @@
+/*
+ * ClusterCollectionEnergies.java
+ *
+ * Created on September 8, 2007, 12:09 PM
+ *
+ * Driver to store and optionally plot energies from a cluster collection as a
+ * function of particle type
+ */
+
+package org.lcsim.contrib.Cassell.recon.analysis;
+import java.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.aida.AIDA;
+import java.util.Map.Entry;
+
+/**
+ *
+ * @author cassell
+ */
+public class ClusterCollectionEnergies extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ String coll;
+ String pre;
+ String FSname;
+ String[] type;
+ int nmax = 10000000;
+ double[] totE;
+ double[] runE;
+ ParticleType pt;
+ CalClusterComposition ccc;
+ boolean makeHists;
+
+ /** Creates a new instance of ClusterCollectionEnergies */
+ public ClusterCollectionEnergies(String fsname,String collname,boolean plot)
+ {
+ pt = new ParticleType();
+ type = pt.getTypes();
+ coll = collname;
+ FSname = fsname;
+ makeHists = plot;
+ runE = new double[type.length];
+ pre = coll+"/";
+ }
+ public ClusterCollectionEnergies(String fsname,String collname)
+ {
+ this(fsname,collname,false);
+ }
+ protected void process(EventHeader event)
+ {
+ List<MCParticle> fslist = event.get(MCParticle.class,FSname);
+ ccc = new CalClusterComposition(fslist);
+ totE = new double[type.length];
+ List<Cluster> cl = event.get(Cluster.class,coll);
+ for(Cluster c:cl)
+ {
+ ccc.decomposeCluster(c);
+ Map<MCParticle,Double>[] Emap = ccc.getEmap();
+ double[] Etot = ccc.getEnergies();
+ for(int i=0;i<type.length;i++)
+ {
+ totE[i] += Etot[i];
+ runE[i] += Etot[i];
+ }
+ if(makeHists)
+ {
+ double[] fracs = ccc.getFractions();
+ aida.histogram1D(pre+"Cluster Energy", 200, 0., 200.).fill(c.getEnergy());
+ aida.histogram1D(
+ pre+"wted Fraction of cluster from max cont "+type[ccc.getMaxParticleType()],
+ 21, 0., 1.05).fill(ccc.getMaxParticleE()/c.getEnergy(), c.getEnergy());
+ for(int i=0;i<type.length;i++)
+ {
+ aida.histogram1D(pre+type[i]+" Energy per cluster", 200, 0.,200.);
+ aida.histogram1D(
+ pre+"wted Fraction of cluster from "+type[i],
+ 21, 0., 1.05).fill(fracs[i], c.getEnergy());
+ }
+ }
+ }
+ if(makeHists)
+ {
+ double evte = 0.;
+ for(int i=0;i<type.length;i++)
+ {
+ evte += totE[i];
+ aida.cloud1D(pre+"Energy per event from "+type[i], nmax).fill(totE[i]);
+ }
+ aida.cloud1D(pre+"Total Energy per event", nmax).fill(evte);
+ }
+ }
+ public double[] getEnergies(){return totE;}
+ public double[] getRunEnergies(){return runE;}
+
+}
lcsim/src/org/lcsim/contrib/Cassell/recon/analysis
diff -N ExampleAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ExampleAnalysisDriver.java 13 Sep 2007 15:10:11 -0000 1.1
@@ -0,0 +1,165 @@
+/*
+ * ExampleAnalysisDriver.java
+ *
+ * Created on September 8, 2007, 6:46 PM
+ *
+ * Example of some simple analysis using the analysis package
+ */
+
+package org.lcsim.contrib.Cassell.recon.analysis;
+import java.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.contrib.Cassell.recon.Cheat.*;
+import hep.physics.vec.*;
+import org.lcsim.util.aida.AIDA;
+import hep.aida.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class ExampleAnalysisDriver extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+//
+// Collection names
+//
+ String GenFSname = "GenFinalStateParticles";
+ String CheatReconRname = "ReconPerfectReconParticles";
+ String CheatReconFSname = "ReconFSParticles";
+ String PPRPflowRname = "PPRReconParticles";
+ String[] DigisimCollectionNames =
+ {"EcalBarrDigiHits", "EcalEndcapDigiHits","HcalBarrDigiHits", "HcalEndcapDigiHits"};
+ String[] ClustersFromReconnames;
+ String[] ClusterFromCheatReconRname;
+ String[] ClusterFromPPRPflowRname;
+//
+// Analysis Drivers
+//
+ FSParticleCollectionEnergies GenFSEnergies;
+ FSParticleCollectionEnergies ReconFSEnergies;
+ CalorimeterHitCollectionEnergies AllCalHitEnergies;
+ ClusterCollectionEnergies[] PerfectClusterEnergies;
+ ClusterCollectionEnergies[] PPRClusterEnergies;
+ ReconstructedParticleCollectionEnergies PerfectReconEnergies;
+ ReconstructedParticleCollectionEnergies PPRReconEnergies;
+
+ String[] type;
+ ParticleType pt;
+
+ int ievt;
+ int nmax = 1000000;
+ /** Creates a new instance of ExampleAnalysisDriver */
+ public ExampleAnalysisDriver()
+ {
+//
+// Do the cheating reconstruction (includes Digisim)
+//
+ CheatReconDriver crd = new CheatReconDriver();
+ crd.setCheatReconstructedParticleOutputName(CheatReconRname);
+ crd.setCheatFSParticleOutputName(CheatReconFSname);
+ add(crd);
+//
+// Make the perfect pattern recognition pflow reconstructed particles from the
+// Cheat Recon particles
+//
+ add(new PPRDriver(CheatReconRname,PPRPflowRname));
+//
+// Get the particle type strings
+//
+ pt = new ParticleType();
+ type = pt.getRTypes();
+//
+// Add the analysis classes
+//
+ GenFSEnergies = new FSParticleCollectionEnergies(GenFSname,true);
+ add(GenFSEnergies);
+ ReconFSEnergies = new FSParticleCollectionEnergies(CheatReconFSname,true);
+ add(ReconFSEnergies);
+ AllCalHitEnergies = new CalorimeterHitCollectionEnergies(
+ "AllCalHits",CheatReconFSname,DigisimCollectionNames,true);
+ add(AllCalHitEnergies);
+ PerfectReconEnergies = new ReconstructedParticleCollectionEnergies(
+ CheatReconFSname,CheatReconRname, true);
+ add(PerfectReconEnergies);
+ PPRReconEnergies = new ReconstructedParticleCollectionEnergies(
+ CheatReconFSname,PPRPflowRname, true);
+ add(PPRReconEnergies);
+
+//
+// Get the cluster collection names generated from ReconstructedParticles
+// and add the analysis classes
+//
+ ClusterFromCheatReconRname = new String[type.length];
+ ClusterFromPPRPflowRname = new String[type.length];
+ PerfectClusterEnergies = new ClusterCollectionEnergies[type.length];
+ PPRClusterEnergies = new ClusterCollectionEnergies[type.length];
+ for(int i=0;i<type.length;i++)
+ {
+ ClusterFromCheatReconRname[i] = CheatReconRname+type[i]+"Clusters";
+ ClusterFromPPRPflowRname[i] = PPRPflowRname+type[i]+"Clusters";
+ PerfectClusterEnergies[i] = new ClusterCollectionEnergies(
+ CheatReconFSname,ClusterFromCheatReconRname[i],true);
+ add(PerfectClusterEnergies[i]);
+ PPRClusterEnergies[i] = new ClusterCollectionEnergies(
+ CheatReconFSname,ClusterFromPPRPflowRname[i],true);
+ add(PPRClusterEnergies[i]);
+ }
+ ievt = 0;
+ }
+ protected void process(EventHeader event)
+ {
+//
+// Process events with both quarks costheta < .8
+//
+ List<MCParticle> mcl = event.get(MCParticle.class,"MCParticle");
+ boolean keepit = true;
+ double Zmass = 0.;
+ for(MCParticle p:mcl)
+ {
+ int id = Math.abs(p.getPDGID());
+ if( (id == 1 )||(id == 2 )||(id == 3 ) )
+ {
+ if(p.getParents().get(0).getPDGID() == 23 )
+ {
+ Zmass = p.getParents().get(0).getMass();
+ Hep3Vector pp = p.getMomentum();
+ double cost = Math.abs(pp.z()/pp.magnitude());
+ if(cost > .8)
+ {
+ keepit = false;
+ continue;
+ }
+ }
+ }
+ }
+ if(!keepit)
+ {
+ ievt++;
+ return;
+ }
+ System.out.println("Processing event "+ievt);
+ super.process(event);
+ aida.cloud1D(GenFSname+"/Delta Zmass",nmax).fill(GenFSEnergies.getNoNeutrinoMass()-Zmass);
+ aida.cloud1D(CheatReconFSname+"/Delta Zmass",nmax).fill(ReconFSEnergies.getNoNeutrinoMass()-Zmass);
+ aida.cloud1D(CheatReconRname+"/Delta Zmass",nmax).fill(PerfectReconEnergies.getMass()-Zmass);
+ aida.cloud1D(PPRPflowRname+"/Delta Zmass",nmax).fill(PPRReconEnergies.getMass()-Zmass);
+ String pre = PPRPflowRname+" to "+CheatReconRname+"/";
+ double[] Echeat = PerfectReconEnergies.getEnergies();
+ double[] Eppr = PPRReconEnergies.getEnergies();
+ for(int i=0;i<type.length;i++)
+ {
+ if(Echeat[i] > 1.)
+ {
+ double dE = Eppr[i] - Echeat[i];
+ aida.cloud1D(pre+"delta E "+type[i], nmax).fill(dE);
+ aida.cloud1D(pre+"delta E over E "+type[i], nmax).fill(dE/Echeat[i]);
+ aida.cloud1D(pre+"delta E over rootE "+type[i], nmax).fill(dE/Math.sqrt(Echeat[i]));
+ }
+ }
+ ievt++;
+ }
+
+}
lcsim/src/org/lcsim/contrib/Cassell/recon/analysis
diff -N FSParticleCollectionEnergies.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FSParticleCollectionEnergies.java 13 Sep 2007 15:10:11 -0000 1.1
@@ -0,0 +1,104 @@
+/*
+ * FSParticleCollectionEnergies.java
+ *
+ * Created on September 8, 2007, 7:38 AM
+ *
+ * Driver to accumulate and optionally histogram energies from a final state
+ * particle collection as a function of particle type
+ */
+
+package org.lcsim.contrib.Cassell.recon.analysis;
+import java.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author cassell
+ */
+public class FSParticleCollectionEnergies extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ String FSname;
+ String pre;
+ String[] type;
+ double[] totE;
+ double[] runE;
+ int nmax = 10000000;
+ ParticleType pt;
+ boolean makeHists;
+ double mass;
+ double nonmass;
+
+ /** Creates a new instance of FSParticleCollectionEnergies */
+ public FSParticleCollectionEnergies(String fsname, boolean hist)
+ {
+ FSname = fsname;
+ pt = new ParticleType();
+ type = pt.getTypes();
+ pre = FSname + "/";
+ makeHists = hist;
+ runE = new double[type.length];
+ }
+ public FSParticleCollectionEnergies(String fsname)
+ {
+ this(fsname, false);
+ }
+ protected void process(EventHeader event)
+ {
+ List<MCParticle> fslist = event.get(MCParticle.class,FSname);
+ totE = new double[type.length];
+ double px = 0.;
+ double py = 0.;
+ double pz = 0.;
+ double nonpx = 0.;
+ double nonpy = 0.;
+ double nonpz = 0.;
+ for(MCParticle p:fslist)
+ {
+ int itype = pt.getType(p);
+ totE[itype] += p.getEnergy();
+ runE[itype] += p.getEnergy();
+ Hep3Vector m = p.getMomentum();
+ px += m.x();
+ py += m.y();
+ pz += m.z();
+ if(type[itype].compareTo("neutrino") != 0)
+ {
+ nonpx += m.x();
+ nonpy += m.y();
+ nonpz += m.z();
+ }
+ if(makeHists)aida.histogram1D(pre+type[itype]+" Energy distibution of particles",200,0.,100.).fill(p.getEnergy());
+ }
+ double Etot = 0.;
+ double nonEtot = 0.;
+ for(int i=0;i<type.length;i++)
+ {
+ if(makeHists)aida.cloud1D(pre+type[i]+" Energy distibution per event",nmax).fill(totE[i]);
+ Etot += totE[i];
+ if(type[i].compareTo("neutrino") != 0)nonEtot += totE[i];
+ }
+ mass = Math.sqrt(Etot*Etot - px*px - py*py - pz*pz);
+ nonmass = Math.sqrt(nonEtot*nonEtot - nonpx*nonpx - nonpy*nonpy - nonpz*nonpz);
+ if(makeHists)
+ {
+ aida.cloud1D(pre+"Total energy per event",nmax).fill(Etot);
+ aida.cloud1D(pre+"Total noneutrino energy per event",nmax).fill(nonEtot);
+ aida.cloud1D(pre+"Event mass",nmax).fill(mass);
+ aida.cloud1D(pre+"Event noneutrino mass",nmax).fill(nonmass);
+ }
+ }
+ public double[] getEnergies()
+ {return totE;}
+ public double getMass()
+ {return mass;}
+ public double getNoNeutrinoMass()
+ {return nonmass;}
+ public double[] getRunEnergies()
+ {return runE;}
+
+}
lcsim/src/org/lcsim/contrib/Cassell/recon/analysis
diff -N ParticleType.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ParticleType.java 13 Sep 2007 15:10:11 -0000 1.1
@@ -0,0 +1,70 @@
+/*
+ * ParticleType.java
+ *
+ * Created on September 8, 2007, 7:14 AM
+ *
+ * Convenience class to return an integer "type" for a MCParticle
+ */
+
+package org.lcsim.contrib.Cassell.recon.analysis;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+
+/**
+ *
+ * @author cassell
+ */
+public class ParticleType
+{
+
+ String[] type = {"photon","ch hadron","n hadron","electron","neutrino"};
+ String[] rtype = {"photon","charged","n hadron"};
+ int[] ttrtype = {0,1,2,1,2};
+
+ /** Creates a new instance of ParticleType */
+ public ParticleType()
+ {
+ }
+//
+// Return photon=0, charged hadron=1, neutral hadron=2, electron=3, neutrino=4
+//
+ public int getType(MCParticle p)
+ {
+ int itype = 0;
+ int pdg = p.getPDGID();
+ if(pdg != 22)
+ {
+ if(p.getCharge() == 0)
+ {
+ itype = 2;
+ int apdg = Math.abs(pdg);
+ if( (apdg == 12)||(apdg == 14)||(apdg == 16) )itype = 4;
+ }
+ else
+ {
+ itype = 1;
+ if(Math.abs(pdg) == 11)itype = 3;
+ }
+ }
+ return itype;
+ }
+ public int getRType(ReconstructedParticle p)
+ {
+ if(p.getCharge() != 0.)return 1;
+ if(p.getMass() == 0)return 0;
+ return 2;
+ }
+ public String[] getTypes()
+ {
+ return type;
+ }
+ public String[] getRTypes()
+ {
+ return rtype;
+ }
+ public int mapTypeToRType(int itype)
+ {
+ if( (itype >= 0)&&(itype < type.length) )return ttrtype[itype];
+ return -1;
+ }
+}
lcsim/src/org/lcsim/contrib/Cassell/recon/analysis
diff -N ReconstructedParticleCollectionEnergies.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ReconstructedParticleCollectionEnergies.java 13 Sep 2007 15:10:11 -0000 1.1
@@ -0,0 +1,221 @@
+/*
+ * ReconstructedParticleCollectionEnergies.java
+ *
+ * Created on September 8, 2007, 3:08 PM
+ *
+ * Driver to store and optionally plot energies from a ReconstructedParticle
+ * collection as a function of particle type
+ */
+
+package org.lcsim.contrib.Cassell.recon.analysis;
+import java.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.contrib.Cassell.recon.Cheat.TraceOrigin;
+import org.lcsim.util.lcio.LCIOConstants;
+
+/**
+ *
+ * @author cassell
+ */
+public class ReconstructedParticleCollectionEnergies extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ String coll;
+ String pre;
+ String FSname;
+ String[] type;
+ String[] rtype;
+ int nmax = 10000000;
+ TraceOrigin to;
+ double[] totE;
+ double[] runE;
+ double[][] totCalE;
+ double[][] runCalE;
+ double mass;
+ double Etot;
+ boolean makeHists;
+ ParticleType pt;
+
+ /** Creates a new instance of ReconstructedParticleCollectionEnergies */
+ public ReconstructedParticleCollectionEnergies(String fsname,String collname,boolean plot)
+ {
+ coll = collname;
+ pre = coll+"/";
+ FSname = fsname;
+ pt = new ParticleType();
+ type = pt.getTypes();
+ rtype = pt.getRTypes();
+ runE = new double[rtype.length];
+ runCalE = new double[type.length][rtype.length];
+ makeHists = plot;
+ }
+ public ReconstructedParticleCollectionEnergies(String fsname,String collname)
+ {
+ this(fsname,collname,false);
+ }
+ protected void process(EventHeader event)
+ {
+ List<MCParticle> fslist = event.get(MCParticle.class,FSname);
+ to = new TraceOrigin(fslist);
+ List<ReconstructedParticle> rpl = event.get(ReconstructedParticle.class,coll);
+ totE = new double[rtype.length];
+ totCalE = new double[rtype.length][rtype.length];
+ Etot = 0.;
+ double pxtot = 0.;
+ double pytot = 0.;
+ double pztot = 0.;
+ List<Cluster>[] rpcl = new List[type.length];
+ for(int i=0;i<type.length;i++)
+ {
+ rpcl[i] = new ArrayList<Cluster>();
+ }
+ for(ReconstructedParticle rp:rpl)
+ {
+ double Epart = rp.getEnergy();
+ Hep3Vector rP = rp.getMomentum();
+ Etot += Epart;
+ pxtot += rP.x();
+ pytot += rP.y();
+ pztot += rP.z();
+ int irtype = pt.getRType(rp);
+ totE[irtype] += Epart;
+ runE[irtype] += Epart;
+ List<Cluster> cl = rp.getClusters();
+ double totEcal = 0.;
+ if(cl.size()>0)
+ {
+ BasicCluster rpcluster = new BasicCluster();
+ for(Cluster c:cl)
+ {
+ rpcluster.addCluster(c);
+ List<CalorimeterHit> hl = c.getCalorimeterHits();
+ for(CalorimeterHit h:hl)
+ {
+ totEcal += h.getCorrectedEnergy();
+ SimCalorimeterHit hh = (SimCalorimeterHit) h;
+ if(hh.getMCParticleCount() == 1)
+ {
+ MCParticle pp = to.traceit(hh.getMCParticle(0));
+ if(pp == null)
+ {
+ pp = hh.getMCParticle(0);
+ }
+ int itype = pt.getType(pp);
+ int ityper = pt.mapTypeToRType(itype);
+ totCalE[irtype][ityper] += h.getCorrectedEnergy();
+ runCalE[irtype][ityper] += h.getCorrectedEnergy();
+ }
+ else
+ {
+ double totraw = 0.;
+ double[] cont = new double[rtype.length];
+ double fac = hh.getCorrectedEnergy()/hh.getRawEnergy();
+ for(int ii=0;ii<hh.getMCParticleCount();ii++)
+ {
+ MCParticle pp = to.traceit(hh.getMCParticle(ii));
+ if(pp == null)
+ {
+ pp = hh.getMCParticle(ii);
+ }
+ int itype = pt.getType(pp);
+ int ityper = pt.mapTypeToRType(itype);
+ totraw += hh.getContributedEnergy(ii);
+ cont[ityper] += hh.getContributedEnergy(ii);
+ }
+ for(int ii=0;ii<rtype.length;ii++)
+ {
+ totCalE[irtype][ii] += cont[ii]/totraw*h.getCorrectedEnergy();
+ runCalE[irtype][ii] += cont[ii]/totraw*h.getCorrectedEnergy();
+ }
+ }
+ }
+ }
+ rpcl[irtype].add(rpcluster);
+ if(makeHists)
+ {
+ aida.histogram1D(pre+type[irtype]+" Energy distribution", 200, 0., 200.).fill(Epart);
+ aida.histogram1D(pre+type[irtype]+" Cal Energy distribution", 200, 0., 200.).fill(totEcal);
+ }
+ }
+ double dE = totEcal - Epart;
+ if(makeHists)
+ {
+ aida.histogram1D(pre+type[irtype]+" CalE-E", 200, -50., 50.).fill(dE);
+ aida.histogram1D(pre+type[irtype]+" CalE-E over E", 200, -1., 1.).fill(dE/Epart);
+ aida.histogram1D(pre+type[irtype]+" CalE-E over rootE", 200, -4., 4.).fill(dE/Math.sqrt(Epart));
+ }
+ }
+ for(int i=0;i<rtype.length;i++)
+ {
+ int flag = 1<<LCIOConstants.CLBIT_HITS;
+ event.put(coll+rtype[i]+"Clusters",rpcl[i],Cluster.class,flag);
+ }
+ mass = Math.sqrt(Etot*Etot - pxtot*pxtot - pytot*pytot - pztot*pztot);
+ if(makeHists)
+ {
+ for(int i=0;i<rtype.length;i++)
+ {
+ aida.cloud1D(pre+rtype[i]+" Energy distibution per event",nmax).fill(totE[i]);
+ double calE = 0.;;
+ for(int j=0;j<rtype.length;j++)
+ {
+ calE += totCalE[i][j];
+ }
+ if(calE > 0.)
+ {
+ for(int j=0;j<rtype.length;j++)
+ {
+ aida.histogram1D(
+ pre+rtype[i]+":wted fraction of calE from"+rtype[j],
+ 51, 0., 1.02).fill(totCalE[i][j]/calE, calE);
+ }
+ }
+ }
+ aida.cloud1D(pre+"Event Mass",nmax).fill(mass);
+ aida.cloud1D(pre+"Event Energy",nmax).fill(Etot);
+ }
+ }
+ public double[] getEnergies()
+ {return totE;}
+ public double[] getRunEnergies()
+ {return runE;}
+ public double[][] getCalEnergies()
+ {return totCalE;}
+ public double[][] getCalRunEnergies()
+ {return runCalE;}
+ public double getMass()
+ {return mass;}
+ public double getEnergy()
+ {return Etot;}
+ protected void endOfData()
+ {
+ System.out.println(coll+" run summary");
+ System.out.println(" Energy sums");
+ for(int i=0;i<rtype.length;i++)
+ {
+ System.out.println(rtype[i]+" = "+runE[i]);
+ }
+ System.out.println();
+ System.out.println(" Cal Energy sums");
+ for(int i=0;i<rtype.length;i++)
+ {
+ String str = rtype[i];
+ for(int j=0;j<rtype.length;j++)
+ {
+ str += " from "+rtype[j]+" = "+(int)runCalE[i][j]+";";
+ }
+ System.out.println(str);
+ }
+
+ }
+
+}
CVSspam 0.2.8