lcsim/src/org/lcsim/recon/cluster/analysis
diff -N ExampleClusterAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ExampleClusterAnalysis.java 12 Oct 2005 17:02:42 -0000 1.1
@@ -0,0 +1,374 @@
+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;
+ String[] typeNames = {"neutrino","muon","EM","nHad","cHad","All"};
+ ICloud1D[] cAFEnergy;
+ ICloud1D[] cVFEnergy;
+ ICloud1D[] cFFEnergy;
+ 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;
+ ntypes = 5;
+ 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+"/Clusters");
+ tree.mkdir(TreeName+"/Clusters/Fragment");
+ tree.mkdir(TreeName+"/Clusters/Primary");
+ cAFEnergy = new ICloud1D[ntypes+1];
+ cFFEnergy = new ICloud1D[ntypes+1];
+ cVFEnergy = 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+"/Particles/"+typeNames[i]);
+ tree.cd(TreeName+"/Particles/"+typeNames[i]);
+ 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");
+ 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)
+ {
+ List<MCPClusterInfo> pc = (List<MCPClusterInfo>) event.get(MCPCname);
+ List<ClusterMCPInfo> cp = (List<ClusterMCPInfo>) event.get(CMCPname);
+ for(MCPClusterInfo pci:pc)
+ {
+ MCParticle p = pci.getMCParticle();
+ double pE = p.getEnergy();
+ boolean vis = pci.getTotalHits() > 0;
+ boolean found = false;
+ if(pci.getNClusters() > 0)
+ {
+ if(!(pci.getMaxCMCP() == null))
+ {
+ found = pci.getMaxCMCP().getMaxContributer() == p;
+ }
+ }
+ int pdgid = Math.abs(p.getPDGID());
+ int itype = -1;
+ if( (pdgid == 12)||(pdgid == 14)||(pdgid == 16) )
+ {
+// neutrinos
+ itype = 0;
+ }
+ else if( (pdgid == 13) )
+ {
+// Muon
+ itype = 1;
+ }
+ else if( (pdgid == 11)||(pdgid == 22) )
+ {
+// EM
+ itype = 2;
+ }
+ else if( p.getCharge() == 0 )
+ {
+// neutral Had
+ itype = 3;
+ }
+ else
+ {
+// charged Had
+ itype = 4;
+ }
+ if(!pci.isFinalState())
+ {
+ cAXEnergy[itype].fill(pE);
+ cAXEnergyw[itype].fill(pE, pE);
+ cAXEnergy[5].fill(pE);
+ cAXEnergyw[5].fill(pE, pE);
+ if(vis)
+ {
+ cVXEnergy[itype].fill(pE);
+ cVXEnergyw[itype].fill(pE, pE);
+ cVXEnergy[5].fill(pE);
+ cVXEnergyw[5].fill(pE,pE);
+ }
+ if(found)
+ {
+ cFXEnergy[itype].fill(pE);
+ cFXEnergyw[itype].fill(pE, pE);
+ cFXEnergy[5].fill(pE);
+ cFXEnergyw[5].fill(pE, pE);
+ }
+ continue;
+ }
+ cAFEnergy[itype].fill(pE);
+ cAFEnergyw[itype].fill(pE, pE);
+ cAFEnergy[5].fill(pE);
+ cAFEnergyw[5].fill(pE, pE);
+ if(vis)
+ {
+ cVFEnergy[itype].fill(pE);
+ cVFEnergyw[itype].fill(pE, pE);
+ cVFEnergy[5].fill(pE);
+ cVFEnergyw[5].fill(pE,pE);
+ }
+ if(found)
+ {
+ cFFEnergy[itype].fill(pE);
+ cFFEnergyw[itype].fill(pE, pE);
+ cFFEnergy[5].fill(pE);
+ cFFEnergyw[5].fill(pE, 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[5][j].fill(pE,1.);
+ pFEffv[5][j].fill(pE,1.);
+ pFPur[5][j].fill(pE,1.);
+ }
+ }
+ else
+ {
+ for(int j=0;j<4;j++)
+ {
+ pFEff[itype][j].fill(pE,0.);
+ pFEff[5][j].fill(pE,0.);
+ if(vis)
+ {
+ pFEffv[itype][j].fill(pE,0.);
+ pFEffv[5][j].fill(pE,0.);
+ }
+ }
+ }
+ }
+ for(ClusterMCPInfo cpi:cp)
+ {
+ BasicCluster c = cpi.getCluster();
+ boolean isprim = cpi.getMaxMCPC().getMaxCluster() == c;
+ int pdgid = Math.abs(cpi.getMaxContributer().getPDGID());
+ double charge = cpi.getMaxContributer().getCharge();
+ int itype = -1;
+ if( (pdgid == 12)||(pdgid == 14)||(pdgid == 16) )
+ {
+// neutrinos
+ itype = 0;
+ }
+ else if( (pdgid == 13) )
+ {
+// Muon
+ itype = 1;
+ }
+ else if( (pdgid == 11)||(pdgid == 22) )
+ {
+// EM
+ itype = 2;
+ }
+ else if( charge == 0 )
+ {
+// neutral Had
+ itype = 3;
+ }
+ else
+ {
+// charged Had
+ itype = 4;
+ }
+ 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[5].fill(cE);
+ pPclEffh[5].fill(Escale,neff);
+ pPclEffe[5].fill(Escale,eeff);
+ pPclPurh[5].fill(Escale,npur);
+ pPclPure[5].fill(Escale,epur);
+ cPclEffh[5].fill(neff);
+ cPclEffe[5].fill(eeff);
+ cPclPurh[5].fill(npur);
+ cPclPure[5].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[5].fill(cE);
+ pFclEffh[5].fill(cE,neff);
+ pFclEffe[5].fill(cE,eeff);
+ pFclPurh[5].fill(cE,npur);
+ pFclPure[5].fill(cE,epur);
+ cFclEffh[5].fill(neff);
+ cFclEffe[5].fill(eeff);
+ cFclPurh[5].fill(npur);
+ cFclPure[5].fill(epur);
+ for(int j=0;j<4;j++)
+ {
+ pFPur[5][j].fill(cE,0.);
+ }
+ }
+ }
+ ievt ++;
+ }
+}