3 added + 1 removed + 3 modified, total 7 files
lcsim/src/org/lcsim/recon/cluster/analysis
diff -N ClusterAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ClusterAnalysis.java 17 Feb 2006 17:11:36 -0000 1.1
@@ -0,0 +1,11 @@
+package org.lcsim.recon.cluster.analysis;
+import org.lcsim.event.EventHeader;
+
+/**
+ * Interface for doing cluster analysis
+ * @author Ron Cassell
+ */
+public interface ClusterAnalysis
+{
+ public void analyzeEvent(EventHeader event);
+}
lcsim/src/org/lcsim/recon/cluster/analysis
diff -N ClusterAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ClusterAnalysisDriver.java 17 Feb 2006 17:11:36 -0000 1.1
@@ -0,0 +1,114 @@
+package org.lcsim.recon.cluster.analysis;
+import java.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.util.*;
+
+/**
+ * Cluster analysis driver
+ */
+public class ClusterAnalysisDriver extends Driver
+{
+ String fsParticleListName;
+ String[] hitcollnames;
+ String[] clustercollnames;
+ String MCPCname;
+ String CMCPname;
+ String foldername;
+ boolean first;
+ CreateClusterAnalysisLists crcllists;
+ ClusterAnalysis analyzer;
+ static String[] defhitcolls = {"EcalBarrHits","EcalEndcapHits","HcalBarrHits",
+ "HcalEndcapHits"};
+ static String defFSPlist = "SimFinalStateParticles";
+ CreateFinalStateMCParticleList fsDriver;
+ boolean gotit;
+ public ClusterAnalysisDriver(String[] clnames)
+ {
+ this(clnames,defhitcolls,defFSPlist, clnames[0]);
+ }
+ public ClusterAnalysisDriver(String[] clnames, String[] hcnames)
+ {
+ this(clnames,hcnames,defFSPlist, clnames[0]);
+ }
+ public ClusterAnalysisDriver(String[] clnames, String[] hcnames,
+ String fsPname)
+ {
+ this(clnames,hcnames,fsPname, clnames[0]);
+ }
+ public ClusterAnalysisDriver(String[] clnames, String[] hcnames,
+ String fsPname, String fname)
+ {
+ first = true;
+ clustercollnames = clnames;
+ hitcollnames = hcnames;
+ MCPCname = clnames[0]+"MCPC";
+ CMCPname = clnames[0]+"CMCP";
+ foldername = fname;
+ fsParticleListName = fsPname;
+ String hclist = hitcollnames[0];
+ for(int i=1;i<hitcollnames.length;i++)
+ {
+ hclist += ", "+hitcollnames[i];
+ }
+ String cllist = clustercollnames[0];
+ for(int i=1;i<clustercollnames.length;i++)
+ {
+ cllist += ", "+clustercollnames[i];
+ }
+ System.out.println("ClusterAnalysisDriver initialized");
+ System.out.println("Hit lists = "+hclist);
+ System.out.println("Cluster lists = "+cllist);
+ System.out.println("FSParticle list = "+fsParticleListName);
+ System.out.println("Folder name = "+foldername);
+ DefaultClusterAnalysis defanal = new DefaultClusterAnalysis(MCPCname,
+ CMCPname, foldername);
+ analyzer = defanal;
+ gotit = false;
+ crcllists = new CreateClusterAnalysisLists(hitcollnames,clustercollnames,
+ fsParticleListName,MCPCname,CMCPname);
+ }
+ public void setAnalyzer(ClusterAnalysis anal)
+ {
+ analyzer = anal;
+ }
+ protected void process(EventHeader event)
+ {
+ if(first)
+ {
+ gotit = false;
+ first = false;
+ List<List<MCParticle>> mclists = event.get(MCParticle.class);
+ for(List<MCParticle> mclist:mclists)
+ {
+ if(event.getMetaData(mclist).getName().compareTo(fsParticleListName) == 0)
+ gotit = true;
+ }
+ if(!gotit)
+ {
+ fsDriver = new CreateFinalStateMCParticleList("Sim");
+ fsDriver.setCollectionName(fsParticleListName);
+ fsDriver.setKeepContinuousElectrons();
+ fsDriver.setKeepContinuousPhotons();
+ fsDriver.setKeepContinuousHadrons();
+ }
+ }
+
+//
+// Create the Final state particle list, if needed
+//
+ if(!gotit)
+ {
+ fsDriver.process(event);
+ }
+//
+// Create the analysis lists
+//
+ crcllists.CreateLists(event);
+//
+// Do the analysis
+//
+ analyzer.analyzeEvent(event);
+ }
+}
lcsim/src/org/lcsim/recon/cluster/analysis
diff -N DefaultClusterAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DefaultClusterAnalysis.java 17 Feb 2006 17:11:36 -0000 1.1
@@ -0,0 +1,553 @@
+package org.lcsim.recon.cluster.analysis;
+
+import java.util.*;
+import hep.aida.*;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.event.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.recon.cluster.cheat.*;
+import org.lcsim.event.Cluster;
+import hep.physics.vec.*;
+
+/**
+ * Default class to generate a multitude of plots analyzing
+ * clusters.
+ * @author Ron Cassell
+ */
+public class DefaultClusterAnalysis implements ClusterAnalysis
+{
+ int ievt = 0;
+ String MCPCname;
+ String CMCPname;
+ String Treename;
+ int ntypes = 6;
+ String[] typeNames = {"neutrino","muon","electron","gamma","nHad","cHad","All"};
+ // Event sums
+ double[] EvtGenEall;
+ double[] EvtGenEvis;
+ double[] EvtGenEfound;
+ double[] EvtVisEvis;
+ double[] EvtVisEfound;
+ double[] EvtMeasEprim;
+ double[] EvtMeasEfrag;
+ double[] EvtMeasEuncl;
+ double[] EvtNhitsAll;
+ double[] EvtNhitsClustered;
+ double[] EvtNhitsPrim;
+ double[] EvtNhitsFrag;
+ double[] EvtNhitsUncl;
+ double[] EvtNClusters;
+ double[] EvtNFragClusters;
+ //
+ //Event plots
+ ICloud1D[] cEvtGenFSE;
+ ICloud1D[] cEvtGenFSEvis;
+ ICloud1D[] cEvtGenFSEfound;
+ ICloud1D[] cEvtVisEvis;
+ ICloud1D[] cEvtVisEfound;
+ ICloud1D[] cEvtMeasEprim;
+ ICloud1D[] cEvtMeasEfrag;
+ ICloud1D[] cEvtMeasEuncl;
+ ICloud1D[] cEvtNhitsAll;
+ ICloud1D[] cEvtNhitsClustered;
+ ICloud1D[] cEvtNhitsPrim;
+ ICloud1D[] cEvtNhitsFrag;
+ ICloud1D[] cEvtNhitsUncl;
+ ICloud1D[] cEvtFracFragCl;
+ ICloud1D[] cEvtFracHitsUncl;
+ ICloud1D[] cEvtFracEPrim;
+ ICloud1D[] cEvtFracEFrag;
+ ICloud1D[] cEvtFracEUncl;
+ ICloud1D[] cEvtFracGenEvis;
+ ICloud1D[] cEvtFracGenEfound;
+ ICloud1D[] cEvtNClusters;
+ //
+ // Cluster plots
+ ICloud1D[] cClPrimFracPartE;
+ ICloud1D[] cClPrimFracPartN;
+ ICloud1D[] cClPrimFracClE;
+ ICloud1D[] cClPrimFracClN;
+ ICloud1D[] cClFragFracPartE;
+ ICloud1D[] cClFragFracPartN;
+ ICloud1D[] cClFragFracClE;
+ ICloud1D[] cClFragFracClN;
+ IProfile1D[] pClPrimFracPartE;
+ IProfile1D[] pClPrimFracPartN;
+ IProfile1D[] pClPrimFracClE;
+ IProfile1D[] pClPrimFracClN;
+ IProfile1D[] pClFragFracPartE;
+ IProfile1D[] pClFragFracPartN;
+ IProfile1D[] pClFragFracClE;
+ IProfile1D[] pClFragFracClN;
+ //
+ //Particle plots
+ IHistogram1D[] hPartGenEall;
+ IHistogram1D[] hPartGenEvis;
+ IHistogram1D[] hPartGenEfound;
+ ICloud1D[] cPartGenEall;
+ ICloud1D[] cPartGenEvis;
+ ICloud1D[] cPartGenEfound;
+ IHistogram1D[] hPartCosall;
+ IHistogram1D[] hPartCosvis;
+ IHistogram1D[] hPartCosfound;
+ IHistogram1D[] hPartGenEvisoverall;
+ IHistogram1D[] hPartGenEfoundovervis;
+ ITree tree;
+ IHistogramFactory histogramFactory;
+ String TreeName;
+ int plotlevel;
+ //
+ /**
+ *
+ */
+ public DefaultClusterAnalysis(String mcpcname,String cmcpname,String treename)
+ {
+ this(mcpcname,cmcpname,treename,100);
+ }
+ public DefaultClusterAnalysis(String mcpcname,String cmcpname,String treename,
+ int level)
+ {
+ plotlevel = level;
+ MCPCname = mcpcname;
+ CMCPname = cmcpname;
+ Treename = treename;
+ IAnalysisFactory analysisFactory = IAnalysisFactory.create();
+ tree = analysisFactory.createTreeFactory().create();
+ histogramFactory = analysisFactory.createHistogramFactory(tree);
+ TreeName = "/"+Treename;
+ tree.mkdir(TreeName);
+ tree.mkdir(TreeName+"/GeneratedParticleProperties");
+ tree.mkdir(TreeName+"/EventProperties");
+ tree.mkdir(TreeName+"/ClusterProperties");
+
+ //Event sums
+ EvtGenEall = new double[ntypes+1];
+ EvtGenEvis = new double[ntypes+1];
+ EvtGenEfound = new double[ntypes+1];
+ EvtVisEvis = new double[ntypes+1];
+ EvtVisEfound = new double[ntypes+1];
+ EvtMeasEprim = new double[ntypes+1];
+ EvtMeasEfrag = new double[ntypes+1];
+ EvtMeasEuncl = new double[ntypes+1];
+ EvtNhitsAll = new double[ntypes+1];
+ EvtNhitsClustered = new double[ntypes+1];
+ EvtNhitsPrim = new double[ntypes+1];
+ EvtNhitsFrag = new double[ntypes+1];
+ EvtNhitsUncl = new double[ntypes+1];
+ EvtNClusters = new double[ntypes+1];
+ EvtNFragClusters = new double[ntypes+1];
+ //
+ //Event plots
+ cEvtGenFSE = new ICloud1D[ntypes+1];
+ cEvtGenFSEvis = new ICloud1D[ntypes+1];
+ cEvtGenFSEfound = new ICloud1D[ntypes+1];
+ cEvtVisEvis = new ICloud1D[ntypes+1];
+ cEvtVisEfound = new ICloud1D[ntypes+1];
+ cEvtMeasEprim = new ICloud1D[ntypes+1];
+ cEvtMeasEfrag = new ICloud1D[ntypes+1];
+ cEvtMeasEuncl = new ICloud1D[ntypes+1];
+ cEvtNhitsAll = new ICloud1D[ntypes+1];
+ cEvtNhitsClustered = new ICloud1D[ntypes+1];
+ cEvtNhitsPrim = new ICloud1D[ntypes+1];
+ cEvtNhitsFrag = new ICloud1D[ntypes+1];
+ cEvtNhitsUncl = new ICloud1D[ntypes+1];
+ cEvtFracFragCl = new ICloud1D[ntypes+1];
+ cEvtFracHitsUncl = new ICloud1D[ntypes+1];
+ cEvtFracEPrim = new ICloud1D[ntypes+1];
+ cEvtFracEFrag = new ICloud1D[ntypes+1];
+ cEvtFracEUncl = new ICloud1D[ntypes+1];
+ cEvtFracGenEvis = new ICloud1D[ntypes+1];
+ cEvtFracGenEfound = new ICloud1D[ntypes+1];
+ cEvtNClusters = new ICloud1D[ntypes+1];
+ //
+ // Cluster plots
+ cClPrimFracPartE = new ICloud1D[ntypes+1];
+ cClPrimFracPartN = new ICloud1D[ntypes+1];
+ cClPrimFracClE = new ICloud1D[ntypes+1];
+ cClPrimFracClN = new ICloud1D[ntypes+1];
+ cClFragFracPartE = new ICloud1D[ntypes+1];
+ cClFragFracPartN = new ICloud1D[ntypes+1];
+ cClFragFracClE = new ICloud1D[ntypes+1];
+ cClFragFracClN = new ICloud1D[ntypes+1];
+ pClPrimFracPartE = new IProfile1D[ntypes+1];
+ pClPrimFracPartN = new IProfile1D[ntypes+1];
+ pClPrimFracClE = new IProfile1D[ntypes+1];
+ pClPrimFracClN = new IProfile1D[ntypes+1];
+ pClFragFracPartE = new IProfile1D[ntypes+1];
+ pClFragFracPartN = new IProfile1D[ntypes+1];
+ pClFragFracClE = new IProfile1D[ntypes+1];
+ pClFragFracClN = new IProfile1D[ntypes+1];
+ //
+ //Particle plots
+ hPartGenEall = new IHistogram1D[ntypes+1];
+ hPartGenEvis = new IHistogram1D[ntypes+1];
+ hPartGenEfound = new IHistogram1D[ntypes+1];
+ cPartGenEall = new ICloud1D[ntypes+1];
+ cPartGenEvis = new ICloud1D[ntypes+1];
+ cPartGenEfound = new ICloud1D[ntypes+1];
+ hPartCosall = new IHistogram1D[ntypes+1];
+ hPartCosvis = new IHistogram1D[ntypes+1];
+ hPartCosfound = new IHistogram1D[ntypes+1];
+ hPartGenEvisoverall = new IHistogram1D[ntypes+1];
+ hPartGenEfoundovervis = new IHistogram1D[ntypes+1];
+ //
+
+ for(int i=0;i<ntypes+1;i++)
+ {
+ tree.mkdir(TreeName+"/EventProperties/"+typeNames[i]);
+ tree.cd(TreeName+"/EventProperties/"+typeNames[i]);
+ cEvtGenFSE[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Final State Generated Energy per event");
+ cEvtGenFSEvis[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Visible Generated Energy per event");
+ cEvtGenFSEfound[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Generated Energy of found particles per event");
+ cEvtVisEvis[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Visible Energy measuered per event");
+ cEvtVisEfound[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Found Particles Measured Energy per event");
+ cEvtMeasEprim[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Primary Clusters Measured Energy per event");
+ cEvtMeasEfrag[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fragment Clusters Measured Energy per event");
+ cEvtMeasEuncl[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Unclustered Hits Measured Energy per event");
+ cEvtNhitsAll[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":# Calorimeter Hits per event");
+ cEvtNhitsClustered[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":# Clustered Hits per event");
+ cEvtNhitsPrim[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":# Primary Clustered Hits per event");
+ cEvtNhitsFrag[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":# Fragment Clustered Hits per event");
+ cEvtNhitsUncl[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":# Unclustered Hits per event");
+ cEvtFracFragCl[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of Clusters that are Fragments per event");
+ cEvtFracHitsUncl[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of Hits Unclustered per event");
+ cEvtFracEPrim[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of measusred E in Primaries per event");
+ cEvtFracEFrag[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of measusred E in Fragments per event");
+ cEvtFracEUncl[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of measusred E Unclustered per event");
+ cEvtFracGenEvis[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of Generated E Visible per event");
+ cEvtFracGenEfound[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of visible Gen E Found per event");
+ cEvtNClusters[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":# Clusters per event");
+ tree.mkdir(TreeName+"/GeneratedParticleProperties/"+typeNames[i]);
+ tree.cd(TreeName+"/GeneratedParticleProperties/"+typeNames[i]);
+ cPartGenEall[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Generated Energy of all particles");
+ cPartGenEvis[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Generated Energy of visible particles");
+ cPartGenEfound[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Generated Energy of found particles");
+ hPartGenEall[i] = histogramFactory.createHistogram1D(
+ typeNames[i]+": Generated Energy of all particles",100,0.,100.);
+ hPartGenEvis[i] = histogramFactory.createHistogram1D(
+ typeNames[i]+": Generated Energy of visible particles",100,0.,100.);
+ hPartGenEfound[i] = histogramFactory.createHistogram1D(
+ typeNames[i]+": Generated Energy of found particles",100,0.,100.);
+ hPartCosall[i] = histogramFactory.createHistogram1D(
+ typeNames[i]+":Cos theta of all particles",100,-1.,1.);
+ hPartCosvis[i] = histogramFactory.createHistogram1D(
+ typeNames[i]+":Cos theta of visible particles",100,-1.,1.);
+ hPartCosfound[i] = histogramFactory.createHistogram1D(
+ typeNames[i]+":Cos theta of found particles",100,-1.,1.);
+ tree.mkdir(TreeName+"/ClusterProperties/"+typeNames[i]);
+ tree.cd(TreeName+"/ClusterProperties/"+typeNames[i]);
+ cClPrimFracPartE[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of particle E in Primary cluster");
+ cClPrimFracPartN[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of particle Hits in Primary cluster");
+ cClPrimFracClE[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of Primary cluster E from particle");
+ cClPrimFracClN[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of Primary cluster Hits from particle");
+ cClFragFracPartE[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of particle E in Fragment cluster");
+ cClFragFracPartN[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of particle Hits in Fragment cluster");
+ cClFragFracClE[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of Fragment cluster E from particle");
+ cClFragFracClN[i] = histogramFactory.createCloud1D(
+ typeNames[i]+":Fraction of Fragment cluster Hits from particle");
+ pClPrimFracPartE[i] = histogramFactory.createProfile1D(
+ typeNames[i]+": Fraction of particle E in Primary cluster",100,0.,100.);
+ pClPrimFracPartN[i] = histogramFactory.createProfile1D(
+ typeNames[i]+": Fraction of particle Hits in Primary cluster",100,0.,100.);
+ pClPrimFracClE[i] = histogramFactory.createProfile1D(
+ typeNames[i]+": Fraction of Primary cluster E from particle",100,0.,100.);
+ pClPrimFracClN[i] = histogramFactory.createProfile1D(
+ typeNames[i]+": Fraction of Primary cluster Hits from particle",100,0.,100.);
+ pClFragFracPartE[i] = histogramFactory.createProfile1D(
+ typeNames[i]+": Fraction of particle E in Fragment cluster",40,0.,10.);
+ pClFragFracPartN[i] = histogramFactory.createProfile1D(
+ typeNames[i]+": Fraction of particle Hits in Fragment cluster",40,0.,10.);
+ pClFragFracClE[i] = histogramFactory.createProfile1D(
+ typeNames[i]+": Fraction of Fragment cluster E from particle",40,0.,10.);
+ pClFragFracClN[i] = histogramFactory.createProfile1D(
+ typeNames[i]+": Fraction of Fragment cluster Hits from particle",40,0.,10.);
+ }
+
+ }
+ public void analyzeEvent(EventHeader event)
+ {
+
+ //Zero Event sums
+ for(int i=0;i<ntypes+1;i++)
+ {
+ EvtGenEall[i] = 0.;
+ EvtGenEvis[i] = 0.;
+ EvtGenEfound[i] = 0.;
+ EvtVisEvis[i] = 0.;
+ EvtVisEfound[i] = 0.;
+ EvtMeasEprim[i] = 0.;
+ EvtMeasEfrag[i] = 0.;
+ EvtMeasEuncl[i] = 0.;
+ EvtNhitsAll[i] = 0.;
+ EvtNhitsClustered[i] = 0.;
+ EvtNhitsPrim[i] = 0.;
+ EvtNhitsFrag[i] = 0.;
+ EvtNClusters[i] = 0.;
+ EvtNFragClusters[i] = 0.;
+ EvtNhitsUncl[i] = 0.;
+ }
+ //
+
+ List<SimCalorimeterHit> unclustered = new ArrayList<SimCalorimeterHit>();
+ List<MCPClusterInfo> pc = (List<MCPClusterInfo>) event.get(MCPCname);
+ List<ClusterMCPInfo> cp = (List<ClusterMCPInfo>) event.get(CMCPname);
+ for(MCPClusterInfo pci:pc)
+ {
+ List<SimCalorimeterHit> uncl = pci.getUnclusteredHits();
+ for(SimCalorimeterHit h:uncl)
+ {
+ if(!unclustered.contains(h))unclustered.add(h);
+ }
+ MCParticle p = pci.getMCParticle();
+ int itype = getType(p);
+ double pE = p.getEnergy();
+ double vE = pci.getTotalEcorrected();
+ Hep3Vector pp = p.getMomentum();
+ double ptot = pp.magnitude();
+ double ct = pp.z()/ptot;
+ boolean vis = pci.getTotalHits() > 0;
+ boolean found = false;
+ if(pci.getNClusters() > 0)
+ {
+ if(!(pci.getMaxCMCP() == null))
+ {
+ found = pci.getMaxCMCP().getMaxContributer() == p;
+//
+// Allow an exception for overlapping gammas from the same pi0
+//
+ if(!found)
+ {
+ if(p.getPDGID() == 22)
+ {
+ if(pci.getMaxCMCP().getMaxContributer().getPDGID() == 22)
+ {
+ if(p.getParents().get(0) ==
+ pci.getMaxCMCP().getMaxContributer().getParents().get(0))
+ {
+ found = true;
+ }
+ }
+ }
+ }
+ }
+ }
+ if(!pci.isFinalState())
+ {
+ continue;
+ }
+ cPartGenEall[itype].fill(pE);
+ cPartGenEall[ntypes].fill(pE);
+ hPartGenEall[itype].fill(pE);
+ hPartGenEall[ntypes].fill(pE);
+ hPartCosall[itype].fill(ct);
+ hPartCosall[ntypes].fill(ct);
+ EvtGenEall[itype] += pE;
+ EvtGenEall[ntypes] += pE;
+ if(vis)
+ {
+ cPartGenEvis[itype].fill(pE);
+ cPartGenEvis[ntypes].fill(pE);
+ hPartGenEvis[itype].fill(pE);
+ hPartGenEvis[ntypes].fill(pE);
+ hPartCosvis[itype].fill(ct);
+ hPartCosvis[ntypes].fill(ct);
+ EvtGenEvis[itype] += pE;
+ EvtGenEvis[ntypes] += pE;
+ EvtVisEvis[itype] += vE;
+ EvtVisEvis[ntypes] += vE;
+ }
+ if(found)
+ {
+ cPartGenEfound[itype].fill(pE);
+ cPartGenEfound[ntypes].fill(pE);
+ hPartGenEfound[itype].fill(pE);
+ hPartGenEfound[ntypes].fill(pE);
+ hPartCosfound[itype].fill(ct);
+ hPartCosfound[ntypes].fill(ct);
+ EvtGenEfound[itype] += pE;
+ EvtGenEfound[ntypes] += pE;
+ EvtVisEfound[itype] += vE;
+ EvtVisEfound[ntypes] += vE;
+ }
+ }
+ for(ClusterMCPInfo cpi:cp)
+ {
+ Cluster c = cpi.getCluster();
+ boolean isprim = cpi.getMaxMCPC().getMaxCluster() == c;
+ int itype = getType(cpi.getMaxMCPC().getMCParticle());
+ double cE = c.getEnergy();
+ int cN = c.getSize();
+ double eeff = cpi.getMaxCorrectedEnergy()/cpi.getMaxMCPC().getTotalEcorrected();
+ double neff = (double)(cpi.getMaxNHits())/cpi.getMaxMCPC().getTotalHits();
+ double epur = cpi.getMaxCorrectedEnergy()/c.getEnergy();
+ double npur = (double)(cpi.getMaxNHits())/c.getCalorimeterHits().size();
+ EvtNhitsClustered[itype] += cN;
+ EvtNhitsClustered[ntypes] += cN;
+ EvtNClusters[itype] += 1.;
+ EvtNClusters[ntypes] += 1.;
+ if(isprim)
+ {
+ EvtMeasEprim[itype] += cE;
+ EvtMeasEprim[ntypes] += cE;
+ EvtNhitsPrim[itype] += cN;
+ EvtNhitsPrim[ntypes] += cN;
+ cClPrimFracPartE[itype].fill(eeff);
+ cClPrimFracPartN[itype].fill(neff);
+ cClPrimFracClE[itype].fill(epur);
+ cClPrimFracClN[itype].fill(npur);
+ cClPrimFracPartE[ntypes].fill(eeff);
+ cClPrimFracPartN[ntypes].fill(neff);
+ cClPrimFracClE[ntypes].fill(epur);
+ cClPrimFracClN[ntypes].fill(npur);
+ pClPrimFracPartE[itype].fill(cE,eeff);
+ pClPrimFracPartN[itype].fill(cE,neff);
+ pClPrimFracClE[itype].fill(cE,epur);
+ pClPrimFracClN[itype].fill(cE,npur);
+ pClPrimFracPartE[ntypes].fill(cE,eeff);
+ pClPrimFracPartN[ntypes].fill(cE,neff);
+ pClPrimFracClE[ntypes].fill(cE,epur);
+ pClPrimFracClN[ntypes].fill(cE,npur);
+ }
+ else
+ {
+ EvtNFragClusters[itype] += 1.;
+ EvtNFragClusters[ntypes] += 1.;
+ EvtMeasEfrag[itype] += cE;
+ EvtMeasEfrag[ntypes] += cE;
+ EvtNhitsFrag[itype] += cN;
+ EvtNhitsFrag[ntypes] += cN;
+ cClFragFracPartE[itype].fill(eeff);
+ cClFragFracPartN[itype].fill(neff);
+ cClFragFracClE[itype].fill(epur);
+ cClFragFracClN[itype].fill(epur);
+ cClFragFracPartE[ntypes].fill(eeff);
+ cClFragFracPartN[ntypes].fill(neff);
+ cClFragFracClE[ntypes].fill(epur);
+ cClFragFracClN[ntypes].fill(epur);
+ pClFragFracPartE[itype].fill(cE,eeff);
+ pClFragFracPartN[itype].fill(cE,neff);
+ pClFragFracClE[itype].fill(cE,epur);
+ pClFragFracClN[itype].fill(cE,epur);
+ pClFragFracPartE[ntypes].fill(cE,eeff);
+ pClFragFracPartN[ntypes].fill(cE,neff);
+ pClFragFracClE[ntypes].fill(cE,epur);
+ pClFragFracClN[ntypes].fill(cE,epur);
+ }
+ }
+ for(SimCalorimeterHit h:unclustered)
+ {
+ int itype = getType(h.getMCParticle(0));
+ EvtMeasEuncl[itype] += h.getCorrectedEnergy();
+ EvtNhitsUncl[itype] += 1.;
+ EvtMeasEuncl[ntypes] += h.getCorrectedEnergy();
+ EvtNhitsUncl[ntypes] += 1.;
+ }
+ for(int i=0;i<ntypes+1;i++)
+ {
+ EvtNhitsClustered[i] = EvtNhitsPrim[i] + EvtNhitsFrag[i];
+ EvtNhitsAll[i] = EvtNhitsClustered[i] + EvtNhitsUncl[i];
+ cEvtGenFSE[i].fill(EvtGenEall[i]);
+ cEvtGenFSEvis[i].fill(EvtGenEvis[i]);
+ cEvtGenFSEfound[i].fill(EvtGenEfound[i]);
+ cEvtVisEvis[i].fill(EvtVisEvis[i]);
+ cEvtVisEfound[i].fill(EvtVisEfound[i]);
+ cEvtMeasEprim[i].fill(EvtMeasEprim[i]);
+ cEvtMeasEfrag[i].fill(EvtMeasEfrag[i]);
+ cEvtMeasEuncl[i].fill(EvtMeasEuncl[i]);
+ cEvtNhitsAll[i].fill(EvtNhitsAll[i]);
+ cEvtNhitsClustered[i].fill(EvtNhitsClustered[i]);
+ cEvtNhitsPrim[i].fill(EvtNhitsPrim[i]);
+ cEvtNhitsFrag[i].fill(EvtNhitsFrag[i]);
+ cEvtNhitsUncl[i].fill(EvtNhitsUncl[i]);
+ cEvtFracHitsUncl[i].fill(EvtNhitsUncl[i]/EvtNhitsAll[i]);
+ double totE = EvtMeasEprim[i]+EvtMeasEfrag[i]+EvtMeasEuncl[i];
+ cEvtFracEPrim[i].fill(EvtMeasEprim[i]/totE);
+ cEvtFracEFrag[i].fill(EvtMeasEfrag[i]/totE);
+ cEvtFracEUncl[i].fill(EvtMeasEuncl[i]/totE);
+ cEvtFracGenEvis[i].fill(EvtGenEvis[i]/EvtGenEall[i]);
+ cEvtFracGenEfound[i].fill(EvtGenEfound[i]/EvtGenEvis[i]);
+ cEvtNClusters[i].fill(EvtNClusters[i]);
+ cEvtFracFragCl[i].fill(EvtNFragClusters[i]/EvtNClusters[i]);
+ tree.cd(TreeName+"/GeneratedParticleProperties/"+typeNames[i]);
+ hPartGenEvisoverall[i] = histogramFactory.divide(
+ typeNames[i]+":Generated Energy - visible over all",
+ hPartGenEvis[i], hPartGenEall[i]);
+ hPartGenEfoundovervis[i] = histogramFactory.divide(
+ typeNames[i]+":Generated Energy - found over visible",
+ hPartGenEfound[i], hPartGenEvis[i]);
+ }
+ ievt ++;
+ }
+ private int getType(MCParticle p)
+ {
+ int itype = -1;
+ int pdgid = Math.abs(p.getPDGID());
+ double charge = p.getCharge();
+ if( (pdgid == 12)||(pdgid == 14)||(pdgid == 16) )
+ {
+// neutrinos
+ itype = 0;
+ }
+ else if( (pdgid == 13) )
+ {
+// Muon
+ itype = 1;
+ }
+ else if( (pdgid == 11) )
+ {
+// Electron
+ itype = 2;
+ }
+ else if( (pdgid == 22) )
+ {
+// Gamma
+ itype = 3;
+ }
+ else if( charge == 0 )
+ {
+// neutral Had
+ itype = 4;
+ }
+ else
+ {
+// charged Had
+ itype = 5;
+ }
+ return itype;
+ }
+}
lcsim/src/org/lcsim/recon/cluster/analysis
diff -u -r1.3 -r1.4
--- ClusterMCPInfo.java 18 Oct 2005 13:48:11 -0000 1.3
+++ ClusterMCPInfo.java 17 Feb 2006 17:11:36 -0000 1.4
@@ -7,6 +7,7 @@
import org.lcsim.event.SimCalorimeterHit;
import org.lcsim.event.CalorimeterHit;
import org.lcsim.event.MCParticle;
+import org.lcsim.event.Cluster;
import org.lcsim.recon.cluster.util.*;
/**
@@ -21,7 +22,7 @@
private int maxNhits;
private double maxErawContribution;
private double maxEcorrectedContribution;
- private BasicCluster thisCluster;
+ private Cluster thisCluster;
private MCParticle[] mcparts;
private int nPrimary;
private int nFragment;
@@ -33,7 +34,7 @@
*@param c - cluster for which this object is created
*@param fsmap - map of MCParticle to final state MCParticle
*/
- public ClusterMCPInfo(BasicCluster c, Map<MCParticle,MCParticle> fsmap)
+ public ClusterMCPInfo(Cluster c, Map<MCParticle,MCParticle> fsmap)
{
thisCluster = c;
contributers = new HashMap();
@@ -131,7 +132,7 @@
// Store the maximum contributer to this cluster
//
int nh = co.getNhits();
- if(nh > maxNhits)
+ if(co.getEraw() > maxErawContribution)
{
maxNhits = nh;
maxErawContribution = co.getEraw();
@@ -139,17 +140,6 @@
maxMCP = p;
maxMCPC = ci;
}
- else if(nh == maxNhits)
- {
- if(co.getEraw() > maxErawContribution)
- {
- maxNhits = nh;
- maxErawContribution = co.getEraw();
- maxEcorrectedContribution = co.getEcorrected();
- maxMCP = p;
- maxMCPC = ci;
- }
- }
}
primaries = new MCParticle[nPrimary];
fragments = new MCParticle[nFragment];
@@ -256,7 +246,7 @@
{
return maxEcorrectedContribution;
}
- public BasicCluster getCluster()
+ public Cluster getCluster()
{return thisCluster;}
}
lcsim/src/org/lcsim/recon/cluster/analysis
diff -u -r1.4 -r1.5
--- CreateClusterAnalysisLists.java 19 Oct 2005 10:33:05 -0000 1.4
+++ CreateClusterAnalysisLists.java 17 Feb 2006 17:11:36 -0000 1.5
@@ -8,17 +8,18 @@
import org.lcsim.event.MCParticle;
import org.lcsim.event.CalorimeterHit;
import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.event.Cluster;
import org.lcsim.util.Driver;
import org.lcsim.event.EventHeader;
import hep.physics.vec.*;
/**
- * Driver to create and write to event list of objects referencing MCParticles
+ * Class to create and write to event list of objects referencing MCParticles
* to Clusters, and Clusters to MCParticles.
* @author Ron Cassell
- * @version $Id: CreateClusterAnalysisLists.java,v 1.4 2005/10/19 10:33:05 jstrube Exp $
+ * @version $Id: CreateClusterAnalysisLists.java,v 1.5 2006/02/17 17:11:36 cassell Exp $
*/
-public class CreateClusterAnalysisLists extends Driver
+public class CreateClusterAnalysisLists
{
String[] hitlistnames;
String[] clusterlistnames;
@@ -27,48 +28,58 @@
String cplistname;
/**
*
+ * @param event
* @param calhitlists
* @param clusterlists
* @param fslist
* @param pclist
* @param cplist
*/
- public CreateClusterAnalysisLists(String[] calhitlists,String[] clusterlists,
- String fslist,String pclist,String cplist)
+ public CreateClusterAnalysisLists(String[] calhitlists,
+ String[] clusterlists, String fslist,String pclist,String cplist)
{
hitlistnames = calhitlists;
clusterlistnames = clusterlists;
fslistname = fslist;
pclistname = pclist;
cplistname = cplist;
+ String hclist = hitlistnames[0];
+ for(int i=1;i<hitlistnames.length;i++)
+ {
+ hclist += ", "+hitlistnames[i];
+ }
+ String cllist = clusterlistnames[0];
+ for(int i=1;i<clusterlistnames.length;i++)
+ {
+ cllist += ", "+clusterlistnames[i];
+ }
+ System.out.println("CreateClusterAnalysisLists initialized");
+ System.out.println("Hit lists = "+hclist);
+ System.out.println("Cluster lists = "+cllist);
+ System.out.println("FSParticle list = "+fslistname);
}
- /**
- * Process an event.
- *
- * @param event
- */
- protected void process(EventHeader event)
+ public void CreateLists(EventHeader event)
{
//
// Get the Lists from event
//
List<MCParticle> mcpall = null;
List<MCParticle> mcpfs = null;
- List<BasicCluster> clusters = new ArrayList<BasicCluster>();
+ List<Cluster> clusters = new ArrayList<Cluster>();
List<List<MCParticle>> mclists = event.get(MCParticle.class);
for(List<MCParticle> mclist:mclists)
{
if(event.getMetaData(mclist).getName().compareTo("MCParticle") == 0)mcpall = mclist;
if(event.getMetaData(mclist).getName().compareTo(fslistname) == 0)mcpfs = mclist;
}
- List<List<BasicCluster>> cllists = event.get(BasicCluster.class);
- for(List<BasicCluster> cllist:cllists)
+ List<List<Cluster>> cllists = event.get(Cluster.class);
+ for(List<Cluster> cllist:cllists)
{
for(int i=0;i<clusterlistnames.length;i++)
{
if(event.getMetaData(cllist).getName().compareTo(clusterlistnames[i]) == 0)
{
- for(BasicCluster c:cllist)
+ for(Cluster c:cllist)
{
clusters.add(c);
}
@@ -129,7 +140,7 @@
//
// Loop over all the clusters
//
- for(BasicCluster c:clusters)
+ for(Cluster c:clusters)
{
//
// For each cluster, create a ClusterMCPInfo
lcsim/src/org/lcsim/recon/cluster/analysis
diff -u -r1.3 -r1.4
--- MCPClusterInfo.java 18 Oct 2005 13:48:11 -0000 1.3
+++ MCPClusterInfo.java 17 Feb 2006 17:11:36 -0000 1.4
@@ -6,6 +6,7 @@
import org.lcsim.event.SimCalorimeterHit;
import org.lcsim.event.CalorimeterHit;
import org.lcsim.event.MCParticle;
+import org.lcsim.event.Cluster;
import org.lcsim.recon.cluster.util.*;
/**
@@ -15,11 +16,11 @@
public class MCPClusterInfo
{
private int NClusters;
- private BasicCluster MaxCluster;
+ private Cluster MaxCluster;
private double MaxEraw;
private double MaxEcorrected;
private int MaxNHits;
- private List<BasicCluster> Clusters;
+ private List<Cluster> Clusters;
private List<Double> ErawContribution;
private List<Double> EcorrectedContribution;
private List<Integer> NhitContribution;
@@ -50,7 +51,7 @@
TotalHits = 0;
TotalEraw = 0.;
TotalEcorrected = 0.;
- Clusters = new ArrayList<BasicCluster>();
+ Clusters = new ArrayList<Cluster>();
ErawContribution = new ArrayList<Double>();
EcorrectedContribution = new ArrayList<Double>();
NhitContribution = new ArrayList<Integer>();
@@ -59,7 +60,7 @@
/**
* Add a cluster this particle contributes to
*/
- public void AddCluster(BasicCluster c)
+ public void AddCluster(Cluster c)
{
Clusters.add(c);
NClusters++;
@@ -97,23 +98,13 @@
//
// See if this is the maximum cluster contributing to this particle
//
- if(nhits > MaxNHits)
+ if(ecorr > MaxEcorrected)
{
MaxNHits = nhits;
MaxCluster = c;
MaxEraw = eraw;
MaxEcorrected = ecorr;
}
- else if(nhits == MaxNHits)
- {
- if(ecorr > MaxEcorrected)
- {
- MaxNHits = nhits;
- MaxCluster = c;
- MaxEraw = eraw;
- MaxEcorrected = ecorr;
- }
- }
}
/**
* Add an unclustered hit this particle contributes to
@@ -168,7 +159,7 @@
/**
* Return the cluster with the maximum number of hits contributed by this particle
*/
- public BasicCluster getMaxCluster()
+ public Cluster getMaxCluster()
{return MaxCluster;}
/**
* Return the raw energy contributed by this particle to the cluster with the
@@ -209,7 +200,7 @@
/**
* Return a list of clusters to which this particle contributes
*/
- public List<BasicCluster> getClusters()
+ public List<Cluster> getClusters()
{return Clusters;}
/**
* Return an array of raw energy contributions of this particle to the list of clusters
lcsim/src/org/lcsim/recon/cluster/analysis
diff -N ExampleClusterAnalysis.java
--- ExampleClusterAnalysis.java 20 Oct 2005 12:35:24 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,472 +0,0 @@
-package org.lcsim.recon.cluster.analysis;
-
-import java.util.*;
-import hep.aida.*;
-import org.lcsim.event.MCParticle;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.CalorimeterHit;
-import org.lcsim.event.SimCalorimeterHit;
-import org.lcsim.event.util.*;
-import org.lcsim.util.Driver;
-import org.lcsim.recon.cluster.analysis.*;
-import org.lcsim.recon.cluster.util.*;
-import org.lcsim.recon.cluster.cheat.*;
-import hep.physics.vec.*;
-
-/**
- * Driver to generate a multitude of plots analyzing
- * clusters.
- * @author Ron Cassell
- */
-public class ExampleClusterAnalysis extends Driver
-{
- int ievt = 0;
- String[] hitcollnames;
- String[] clustercollnames;
- String fsParticleListName;
- String MCPCname;
- String CMCPname;
- String Treename;
- int ntypes = 6;
- String[] typeNames = {"neutrino","muon","electron","gamma","nHad","cHad","All"};
- ICloud1D[] cEventFSE;
- ICloud1D[] cEventVE;
- ICloud1D[] cEventFE;
- ICloud1D[] cEventAVE;
- ICloud1D[] cEventAFE;
- IHistogram1D[] hPerPull;
- IHistogram1D[] hPerRes;
- IHistogram1D[] hActPull;
- IHistogram1D[] hActRes;
- ICloud1D[] cAFEnergy;
- ICloud1D[] cVFEnergy;
- ICloud1D[] cFFEnergy;
- ICloud1D[] cAFcosTheta;
- ICloud1D[] cVFcosTheta;
- ICloud1D[] cFFcosTheta;
- ICloud1D[] cAFEnergyw;
- ICloud1D[] cVFEnergyw;
- ICloud1D[] cFFEnergyw;
- ICloud1D[] cAXEnergy;
- ICloud1D[] cVXEnergy;
- ICloud1D[] cFXEnergy;
- ICloud1D[] cAXEnergyw;
- ICloud1D[] cVXEnergyw;
- ICloud1D[] cFXEnergyw;
- IProfile1D[][] pFEff;
- IProfile1D[][] pFPur;
- IProfile1D[][] pFEffv;
- ICloud1D[] cPclEnergy;
- IProfile1D[] pPclEffh;
- IProfile1D[] pPclPurh;
- IProfile1D[] pPclEffe;
- IProfile1D[] pPclPure;
- ICloud1D[] cPclPurh;
- ICloud1D[] cPclEffh;
- ICloud1D[] cPclPure;
- ICloud1D[] cPclEffe;
- ICloud1D[] cFclEnergy;
- IProfile1D[] pFclEffh;
- IProfile1D[] pFclPurh;
- IProfile1D[] pFclEffe;
- IProfile1D[] pFclPure;
- ICloud1D[] cFclPurh;
- ICloud1D[] cFclEffh;
- ICloud1D[] cFclPure;
- ICloud1D[] cFclEffe;
- /**
- *
- */
- public ExampleClusterAnalysis(String[] hcnames,String[] ccnames,String fsPLName,
- String mcpcname,String cmcpname,String treename)
- {
- hitcollnames = hcnames;
- clustercollnames = ccnames;
- fsParticleListName = fsPLName;
- MCPCname = mcpcname;
- CMCPname = cmcpname;
- Treename = treename;
- IAnalysisFactory analysisFactory = IAnalysisFactory.create();
- ITree tree = analysisFactory.createTreeFactory().create();
- IHistogramFactory histogramFactory = analysisFactory.createHistogramFactory(tree);
- String TreeName = "/"+Treename;
- tree.mkdir(TreeName);
- tree.mkdir(TreeName+"/Particles");
- tree.mkdir(TreeName+"/Event");
- tree.mkdir(TreeName+"/Clusters");
- tree.mkdir(TreeName+"/Clusters/Fragment");
- tree.mkdir(TreeName+"/Clusters/Primary");
- cEventFSE = new ICloud1D[ntypes+1];
- cEventVE = new ICloud1D[ntypes+1];
- cEventFE = new ICloud1D[ntypes+1];
- cEventAVE = new ICloud1D[ntypes+1];
- cEventAFE = new ICloud1D[ntypes+1];
- hPerPull = new IHistogram1D[ntypes];
- hPerRes = new IHistogram1D[ntypes];
- hActPull = new IHistogram1D[ntypes];
- hActRes = new IHistogram1D[ntypes];
- cAFEnergy = new ICloud1D[ntypes+1];
- cFFEnergy = new ICloud1D[ntypes+1];
- cVFEnergy = new ICloud1D[ntypes+1];
- cAFcosTheta = new ICloud1D[ntypes+1];
- cFFcosTheta = new ICloud1D[ntypes+1];
- cVFcosTheta = new ICloud1D[ntypes+1];
- cAFEnergyw = new ICloud1D[ntypes+1];
- cFFEnergyw = new ICloud1D[ntypes+1];
- cVFEnergyw = new ICloud1D[ntypes+1];
- cAXEnergy = new ICloud1D[ntypes+1];
- cFXEnergy = new ICloud1D[ntypes+1];
- cVXEnergy = new ICloud1D[ntypes+1];
- cAXEnergyw = new ICloud1D[ntypes+1];
- cFXEnergyw = new ICloud1D[ntypes+1];
- cVXEnergyw = new ICloud1D[ntypes+1];
- pFEff = new IProfile1D[ntypes+1][4];
- pFPur = new IProfile1D[ntypes+1][4];
- pFEffv = new IProfile1D[ntypes+1][4];
- cPclEnergy = new ICloud1D[ntypes+1];
- cFclEnergy = new ICloud1D[ntypes+1];
- pPclEffh = new IProfile1D[ntypes+1];
- pPclEffe = new IProfile1D[ntypes+1];
- pPclPurh = new IProfile1D[ntypes+1];
- pPclPure = new IProfile1D[ntypes+1];
- cPclEffh = new ICloud1D[ntypes+1];
- cPclEffe = new ICloud1D[ntypes+1];
- cPclPurh = new ICloud1D[ntypes+1];
- cPclPure = new ICloud1D[ntypes+1];
- pFclEffh = new IProfile1D[ntypes+1];
- pFclEffe = new IProfile1D[ntypes+1];
- pFclPurh = new IProfile1D[ntypes+1];
- pFclPure = new IProfile1D[ntypes+1];
- cFclEffh = new ICloud1D[ntypes+1];
- cFclEffe = new ICloud1D[ntypes+1];
- cFclPurh = new ICloud1D[ntypes+1];
- cFclPure = new ICloud1D[ntypes+1];
- int[] nbins = {2000,400,200,40};
- for(int i=0;i<ntypes+1;i++)
- {
- tree.mkdir(TreeName+"/Event/"+typeNames[i]);
- tree.cd(TreeName+"/Event/"+typeNames[i]);
- cEventFSE[i] = histogramFactory.createCloud1D(
- typeNames[i]+":Final State Energy per event");
- cEventVE[i] = histogramFactory.createCloud1D(
- typeNames[i]+":Visible Generated Energy per event");
- cEventFE[i] = histogramFactory.createCloud1D(
- typeNames[i]+":Generated Energy of found particles per event");
- cEventAVE[i] = histogramFactory.createCloud1D(
- typeNames[i]+":Energy deposited per event");
- cEventAFE[i] = histogramFactory.createCloud1D(
- typeNames[i]+":Energy Measured per event");
- tree.mkdir(TreeName+"/Particles/"+typeNames[i]);
- tree.cd(TreeName+"/Particles/"+typeNames[i]);
- if(i<ntypes)
- {
- hPerPull[i] = histogramFactory.createHistogram1D(typeNames[i]+
- "Total Edep - E over E ",100, -.5, .5);
- hPerRes[i] = histogramFactory.createHistogram1D(typeNames[i]+
- "Total Edep - E over sqrtE ",100, -4, 4);
- hActPull[i] = histogramFactory.createHistogram1D(typeNames[i]+
- "Cluster Edep - E over E ",100, -.5, .5);
- hActRes[i] = histogramFactory.createHistogram1D(typeNames[i]+
- "Cluster Edep - E over sqrtE ",100, -4, 4);
- }
- cAFEnergy[i] = histogramFactory.createCloud1D("All FS "+typeNames[i]+": Energy");
- cVFEnergy[i] = histogramFactory.createCloud1D("Vis FS "+typeNames[i]+": Energy");
- cFFEnergy[i] = histogramFactory.createCloud1D("Found FS "+typeNames[i]+": Energy");
- cAFcosTheta[i] = histogramFactory.createCloud1D("All FS "+typeNames[i]+": cosTheta");
- cVFcosTheta[i] = histogramFactory.createCloud1D("Vis FS "+typeNames[i]+": cosTheta");
- cFFcosTheta[i] = histogramFactory.createCloud1D("Found FS "+typeNames[i]+": cosTheta");
- cAFEnergyw[i] = histogramFactory.createCloud1D("All FS "+typeNames[i]+": wted Energy");
- cVFEnergyw[i] = histogramFactory.createCloud1D("Vis FS "+typeNames[i]+": wted Energy");
- cFFEnergyw[i] = histogramFactory.createCloud1D("Found FS "+typeNames[i]+": wted Energy");
- cAXEnergy[i] = histogramFactory.createCloud1D("All nonFS "+typeNames[i]+": Energy");
- cVXEnergy[i] = histogramFactory.createCloud1D("Vis nonFS "+typeNames[i]+": Energy");
- cFXEnergy[i] = histogramFactory.createCloud1D("Found nonFS "+typeNames[i]+": Energy");
- cAXEnergyw[i] = histogramFactory.createCloud1D("All nonFS "+typeNames[i]+": wted Energy");
- cVXEnergyw[i] = histogramFactory.createCloud1D("Vis nonFS "+typeNames[i]+": wted Energy");
- cFXEnergyw[i] = histogramFactory.createCloud1D("Found nonFS "+typeNames[i]+": wted Energy");
- for(int j=0;j<4;j++)
- {
- pFEff[i][j] = histogramFactory.createProfile1D(
- "FS "+typeNames[i]+" Eff"+j,"Efficiency",nbins[j],0.,200.);
- pFPur[i][j] = histogramFactory.createProfile1D(
- "FS "+typeNames[i]+" Pur","Purity",nbins[j],0.,200.);
- pFEffv[i][j] = histogramFactory.createProfile1D(
- "Visible "+typeNames[i]+" Eff"+j,"Efficiency",nbins[j],0.,200.);
- }
- tree.mkdir(TreeName+"/Clusters/Primary/"+typeNames[i]);
- tree.cd(TreeName+"/Clusters/Primary/"+typeNames[i]);
- cPclEnergy[i] = histogramFactory.createCloud1D("Primary cluster Energy: "+typeNames[i]);
- pPclEffh[i] = histogramFactory.createProfile1D("Primary cluster hit efficiency",typeNames[i],200,0.,200.);
- pPclPurh[i] = histogramFactory.createProfile1D("Primary cluster hit purity",typeNames[i],200,0.,200.);
- pPclEffe[i] = histogramFactory.createProfile1D("Primary cluster energy efficiency",typeNames[i],200,0.,200.);
- pPclPure[i] = histogramFactory.createProfile1D("Primary cluster energy purity",typeNames[i],200,0.,200.);
- cPclEffh[i] = histogramFactory.createCloud1D(typeNames[i] + "Primary cluster hit efficiency");
- cPclPurh[i] = histogramFactory.createCloud1D(typeNames[i]+"Primary cluster hit purity");
- cPclEffe[i] = histogramFactory.createCloud1D(typeNames[i]+"Primary cluster energy efficiency");
- cPclPure[i] = histogramFactory.createCloud1D(typeNames[i]+"Primary cluster energy purity");
- tree.mkdir(TreeName+"/Clusters/Fragment/"+typeNames[i]);
- tree.cd(TreeName+"/Clusters/Fragment/"+typeNames[i]);
- cFclEnergy[i] = histogramFactory.createCloud1D("Fragment cluster Energy: "+typeNames[i]);
- pFclEffh[i] = histogramFactory.createProfile1D("Fragment cluster hit efficiency",typeNames[i],200,0.,200.);
- pFclPurh[i] = histogramFactory.createProfile1D("Fragment cluster hit purity",typeNames[i],200,0.,200.);
- pFclEffe[i] = histogramFactory.createProfile1D("Fragment cluster energy efficiency",typeNames[i],200,0.,200.);
- pFclPure[i] = histogramFactory.createProfile1D("Fragment cluster energy purity",typeNames[i],200,0.,200.);
- cFclEffh[i] = histogramFactory.createCloud1D(typeNames[i] + "Fragment cluster hit efficiency");
- cFclPurh[i] = histogramFactory.createCloud1D(typeNames[i]+"Fragment cluster hit purity");
- cFclEffe[i] = histogramFactory.createCloud1D(typeNames[i]+"Fragment cluster energy efficiency");
- cFclPure[i] = histogramFactory.createCloud1D(typeNames[i]+"Fragment cluster energy purity");
- }
-
- }
- public ExampleClusterAnalysis(String[] hcnames,String[] ccnames,String fsPLName,
- String mcpcname,String cmcpname)
- {
- this(hcnames,ccnames,fsPLName,mcpcname,cmcpname,"default");
- }
- protected void process(EventHeader event)
- {
- double[] cheatFSE = new double[ntypes+1];
- double[] cheatVE = new double[ntypes+1];
- double[] cheatFE = new double[ntypes+1];
- double[] realVE = new double[ntypes+1];
- double[] realFE = new double[ntypes+1];
- List<CalorimeterHit> unclustered = new ArrayList<CalorimeterHit>();
- for(int i=0;i<ntypes+1;i++)
- {
- cheatFSE[i] = 0.;
- cheatVE[i] = 0.;
- cheatFE[i] = 0.;
- realVE[i] = 0.;
- realFE[i] = 0.;
- }
- List<MCPClusterInfo> pc = (List<MCPClusterInfo>) event.get(MCPCname);
- List<ClusterMCPInfo> cp = (List<ClusterMCPInfo>) event.get(CMCPname);
- for(MCPClusterInfo pci:pc)
- {
- MCParticle p = pci.getMCParticle();
- int itype = getType(p);
- double pE = p.getEnergy();
- double vE = pci.getTotalEcorrected();
- Hep3Vector pp = p.getMomentum();
- double ptot = pp.magnitude();
- double ct = pp.z()/ptot;
- boolean vis = pci.getTotalHits() > 0;
- boolean found = false;
- if(pci.getNClusters() > 0)
- {
- if(!(pci.getMaxCMCP() == null))
- {
- found = pci.getMaxCMCP().getMaxContributer() == p;
-//
-// Allow an exception for overlapping gammas from the same pi0
-//
- if(!found)
- {
- if(p.getPDGID() == 22)
- {
- if(pci.getMaxCMCP().getMaxContributer().getPDGID() == 22)
- {
- if(p.getParents().get(0) ==
- pci.getMaxCMCP().getMaxContributer().getParents().get(0))
- {
- found = true;
- }
- }
- }
- }
- }
- }
- if(!pci.isFinalState())
- {
- cAXEnergy[itype].fill(pE);
- cAXEnergyw[itype].fill(pE, pE);
- cAXEnergy[ntypes].fill(pE);
- cAXEnergyw[ntypes].fill(pE, pE);
- if(vis)
- {
- cVXEnergy[itype].fill(pE);
- cVXEnergyw[itype].fill(pE, pE);
- cVXEnergy[ntypes].fill(pE);
- cVXEnergyw[ntypes].fill(pE,pE);
- }
- if(found)
- {
- cFXEnergy[itype].fill(pE);
- cFXEnergyw[itype].fill(pE, pE);
- cFXEnergy[ntypes].fill(pE);
- cFXEnergyw[ntypes].fill(pE, pE);
- }
- continue;
- }
- cAFEnergy[itype].fill(pE);
- cAFcosTheta[itype].fill(ct);
- cAFEnergyw[itype].fill(pE, pE);
- cAFEnergy[ntypes].fill(pE);
- cAFcosTheta[ntypes].fill(ct);
- cAFEnergyw[ntypes].fill(pE, pE);
- cheatFSE[itype] += pE;
- cheatFSE[ntypes] += pE;
- if(vis)
- {
- cVFEnergy[itype].fill(pE);
- cVFcosTheta[itype].fill(ct);
- cVFEnergyw[itype].fill(pE, pE);
- cVFEnergy[ntypes].fill(pE);
- cVFcosTheta[ntypes].fill(ct);
- cVFEnergyw[ntypes].fill(pE,pE);
- cheatVE[itype] += pE;
- cheatVE[ntypes] += pE;
- realVE[itype] += vE;
- realVE[ntypes] += vE;
- hPerPull[itype].fill((vE-pE)/pE);
- hPerRes[itype].fill((vE-pE)/Math.sqrt(pE));
- }
- if(found)
- {
- cFFEnergy[itype].fill(pE);
- cFFcosTheta[itype].fill(ct);
- cFFEnergyw[itype].fill(pE, pE);
- cFFEnergy[ntypes].fill(pE);
- cFFcosTheta[ntypes].fill(ct);
- cFFEnergyw[ntypes].fill(pE, pE);
- cheatFE[itype] += pE;
- cheatFE[ntypes] += pE;
- realFE[itype] += vE;
- realFE[ntypes] += vE;
- double mE = pci.getMaxCMCP().getCluster().getEnergy();
- hActPull[itype].fill((mE-pE)/pE);
- hActRes[itype].fill((mE-pE)/Math.sqrt(pE));
- for(int j=0;j<4;j++)
- {
- pFEff[itype][j].fill(pE,1.);
- pFEffv[itype][j].fill(pE,1.);
- pFPur[itype][j].fill(pE,1.);
- pFEff[ntypes][j].fill(pE,1.);
- pFEffv[ntypes][j].fill(pE,1.);
- pFPur[ntypes][j].fill(pE,1.);
- }
- }
- else
- {
- for(int j=0;j<4;j++)
- {
- pFEff[itype][j].fill(pE,0.);
- pFEff[ntypes][j].fill(pE,0.);
- if(vis)
- {
- pFEffv[itype][j].fill(pE,0.);
- pFEffv[ntypes][j].fill(pE,0.);
- }
- }
- }
- }
- for(ClusterMCPInfo cpi:cp)
- {
- BasicCluster c = cpi.getCluster();
- boolean isprim = cpi.getMaxMCPC().getMaxCluster() == c;
- int itype = getType(cpi.getMaxMCPC().getMCParticle());
- double cE = c.getEnergy();
- double eeff = cpi.getMaxCorrectedEnergy()/cpi.getMaxMCPC().getTotalEcorrected();
- double neff = (double)(cpi.getMaxNHits())/cpi.getMaxMCPC().getTotalHits();
- double epur = cpi.getMaxCorrectedEnergy()/c.getEnergy();
- double npur = (double)(cpi.getMaxNHits())/c.getCalorimeterHits().size();
- if(isprim)
- {
- double Escale = cpi.getMaxContributer().getEnergy();
- cPclEnergy[itype].fill(cE);
- pPclEffh[itype].fill(Escale,neff);
- pPclEffe[itype].fill(Escale,eeff);
- pPclPurh[itype].fill(Escale,npur);
- pPclPure[itype].fill(Escale,epur);
- cPclEffh[itype].fill(neff);
- cPclEffe[itype].fill(eeff);
- cPclPurh[itype].fill(npur);
- cPclPure[itype].fill(epur);
- cPclEnergy[ntypes].fill(cE);
- pPclEffh[ntypes].fill(Escale,neff);
- pPclEffe[ntypes].fill(Escale,eeff);
- pPclPurh[ntypes].fill(Escale,npur);
- pPclPure[ntypes].fill(Escale,epur);
- cPclEffh[ntypes].fill(neff);
- cPclEffe[ntypes].fill(eeff);
- cPclPurh[ntypes].fill(npur);
- cPclPure[ntypes].fill(epur);
- }
- else
- {
- cFclEnergy[itype].fill(cE);
- pFclEffh[itype].fill(cE,neff);
- pFclEffe[itype].fill(cE,eeff);
- pFclPurh[itype].fill(cE,npur);
- pFclPure[itype].fill(cE,epur);
- cFclEffh[itype].fill(neff);
- cFclEffe[itype].fill(eeff);
- cFclPurh[itype].fill(npur);
- cFclPure[itype].fill(epur);
- for(int j=0;j<4;j++)
- {
- pFPur[itype][j].fill(cE,0.);
- }
- cFclEnergy[ntypes].fill(cE);
- pFclEffh[ntypes].fill(cE,neff);
- pFclEffe[ntypes].fill(cE,eeff);
- pFclPurh[ntypes].fill(cE,npur);
- pFclPure[ntypes].fill(cE,epur);
- cFclEffh[ntypes].fill(neff);
- cFclEffe[ntypes].fill(eeff);
- cFclPurh[ntypes].fill(npur);
- cFclPure[ntypes].fill(epur);
- for(int j=0;j<4;j++)
- {
- pFPur[ntypes][j].fill(cE,0.);
- }
- }
- }
- for(int i=0;i<ntypes+1;i++)
- {
- cEventFSE[i].fill(cheatFSE[i]);
- cEventVE[i].fill(cheatVE[i]);
- cEventFE[i].fill(cheatFE[i]);
- cEventAVE[i].fill(realVE[i]);
- cEventAFE[i].fill(realFE[i]);
- }
- ievt ++;
- }
- private int getType(MCParticle p)
- {
- int itype = -1;
- int pdgid = Math.abs(p.getPDGID());
- double charge = p.getCharge();
- if( (pdgid == 12)||(pdgid == 14)||(pdgid == 16) )
- {
-// neutrinos
- itype = 0;
- }
- else if( (pdgid == 13) )
- {
-// Muon
- itype = 1;
- }
- else if( (pdgid == 11) )
- {
-// Electron
- itype = 2;
- }
- else if( (pdgid == 22) )
- {
-// Gamma
- itype = 3;
- }
- else if( charge == 0 )
- {
-// neutral Had
- itype = 4;
- }
- else
- {
-// charged Had
- itype = 5;
- }
- return itype;
- }
-}
CVSspam 0.2.8