lcsim/src/org/lcsim/recon/cluster/analysis
diff -u -r1.1 -r1.2
--- ExampleClusterAnalysis.java 12 Oct 2005 17:02:42 -0000 1.1
+++ ExampleClusterAnalysis.java 20 Oct 2005 12:35:24 -0000 1.2
@@ -27,11 +27,23 @@
String MCPCname;
String CMCPname;
String Treename;
- int ntypes;
- String[] typeNames = {"neutrino","muon","EM","nHad","cHad","All"};
+ 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;
@@ -74,19 +86,31 @@
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+"/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];
@@ -120,11 +144,37 @@
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");
@@ -175,12 +225,31 @@
}
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)
@@ -188,82 +257,93 @@
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;
+ }
+ }
+ }
+ }
}
}
- 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);
+ cAXEnergy[ntypes].fill(pE);
+ cAXEnergyw[ntypes].fill(pE, pE);
if(vis)
{
cVXEnergy[itype].fill(pE);
cVXEnergyw[itype].fill(pE, pE);
- cVXEnergy[5].fill(pE);
- cVXEnergyw[5].fill(pE,pE);
+ cVXEnergy[ntypes].fill(pE);
+ cVXEnergyw[ntypes].fill(pE,pE);
}
if(found)
{
cFXEnergy[itype].fill(pE);
cFXEnergyw[itype].fill(pE, pE);
- cFXEnergy[5].fill(pE);
- cFXEnergyw[5].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[5].fill(pE);
- cAFEnergyw[5].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[5].fill(pE);
- cVFEnergyw[5].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[5].fill(pE);
- cFFEnergyw[5].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[5][j].fill(pE,1.);
- pFEffv[5][j].fill(pE,1.);
- pFPur[5][j].fill(pE,1.);
+ pFEff[ntypes][j].fill(pE,1.);
+ pFEffv[ntypes][j].fill(pE,1.);
+ pFPur[ntypes][j].fill(pE,1.);
}
}
else
@@ -271,11 +351,11 @@
for(int j=0;j<4;j++)
{
pFEff[itype][j].fill(pE,0.);
- pFEff[5][j].fill(pE,0.);
+ pFEff[ntypes][j].fill(pE,0.);
if(vis)
{
pFEffv[itype][j].fill(pE,0.);
- pFEffv[5][j].fill(pE,0.);
+ pFEffv[ntypes][j].fill(pE,0.);
}
}
}
@@ -284,34 +364,7 @@
{
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;
- }
+ 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();
@@ -329,15 +382,15 @@
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);
+ 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
{
@@ -354,21 +407,66 @@
{
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);
+ 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[5][j].fill(cE,0.);
+ 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;
+ }
}
