lcsim/src/org/lcsim/contrib/SteveKuhlmann
diff -u -r1.1 -r1.2
--- BTrMipClusterID.java 17 Aug 2005 17:20:57 -0000 1.1
+++ BTrMipClusterID.java 19 Aug 2005 14:13:22 -0000 1.2
@@ -62,7 +62,7 @@
private HMatrixBuilder _hmb;
private HMatrix _hmx;
private boolean debug = false;
- private boolean trhist = true;
+ private boolean trhist = false;
private boolean phohist = true;
private boolean readytotryClusterID = true;
@@ -83,7 +83,7 @@
public BTrMipClusterID()
{
- this(new NetworkBasic());
+ this(new NetworkBasic(),HMatrixTask.ANALYZE);
//FIXME: This should depend on detector
myNet.setResourceStream("/hep/lcd/detector/sdmnocalgaps/nn/");
myNet.configureNetwork("network.net");
@@ -99,18 +99,18 @@
// public BTrMipClusterID(HMatrixTask task)
// {
-// _tree = aida.tree();
-// _task = task;
-//
-// if(task==HMatrixTask.ANALYZE)
-// {
-// // The HMatrix could possibly change, so be sensitive to this.
-// getConditionsManager().registerConditionsConverter(new HMatrixConditionsConverter());
-// }
// }
- public BTrMipClusterID(NetworkBasic net)
+ public BTrMipClusterID(NetworkBasic net, HMatrixTask task)
{
myNet = net;
+ _tree = aida.tree();
+ _task = task;
+
+ if(task==HMatrixTask.ANALYZE)
+ {
+ // The HMatrix could possibly change, so be sensitive to this.
+ getConditionsManager().registerConditionsConverter(new HMatrixConditionsConverter());
+ }
}
protected void process(EventHeader event)
@@ -768,12 +768,15 @@
// PHOTON-Finder here (use this QAD one for now)
// check these clusters and try to ID photons - look for track-cluster matches at shower max
// require that EM cluster has hits at this layer.
- if(trhist) aida.cloud1D("Num EM B NNClusters").fill(ebclusters.size());
+ if(phohist) aida.cloud1D("Num EM B NNClusters").fill(ebclusters.size());
List<BasicCluster> phoclusters = new ArrayList();
List<BasicCluster> ebshowclusters = new ArrayList();
- int nGoodClusters = 0;
+ int nGoodClusters = 0; double pho1=0; double pho2=0;
for (BasicCluster ebcluster : ebclusters)
{
+ double energycor = ebcluster.getEnergy();
+ if(energycor >= pho1) pho1=energycor;
+ if(energycor < pho1 && energycor>pho2) pho2 = energycor;
boolean badphoton = true;
// Get theta, phi for cluster
double ecp[] = ebcluster.getPosition();
@@ -818,17 +821,18 @@
double firstDdiff = (firstL + 1)*1.0/length;// *1.0 to force int to double
int maxL=firstL+5;
if(maxL>30) maxL=30;
- for(int i=firstL; i<maxL; i++) aveLE5 += layerE[i];
+ for(int i=firstL; i<maxL; i++) aveLE5 += layerE[i]*ebcluster.getEnergy();
aveLE5 = aveLE5/5.;
int ncells = ebcluster.getCalorimeterHits().size();
aveHits = ncells/length;
- if(ncells>10)
+ if(ncells>4)
{
nGoodClusters++;
double energy = ebcluster.getRawEnergy();
if(phohist) aida.cloud1D("Number of cells in cluster").fill(ncells);
- if(phohist) aida.cloud1D("energy").fill(energy);
+ if(phohist) aida.cloud1D("EM Cluster raw energy").fill(energy);
+ if(phohist) aida.cloud1D("EM Cluster corrected energy").fill(energycor);
// Add the correlation to the log of the cluster energy
//We want the common logarithm (log 10) of the energy
_vals[_logEIndex]=Math.log(energy)*_log10inv;
@@ -848,9 +852,6 @@
// if(phohist) aida.cloud1D("ChisqD Probability").fill(chiprobD);
// double logchiprobD = Math.log(chiprobD)*_log10inv;
// if(phohist) aida.cloud1D("Log ChisqD Probability").fill(logchiprobD);
-//
-// //System.out.println("chisqD "+chisqD+" prob "+ChisqProb.gammq(_nmeas,chisqD));
-//
// double chisqND = chisq - chisqD;
// if(phohist) aida.cloud1D("ChisqND").fill(chisqND);
// if (logchiprobD > -9.5)
@@ -875,24 +876,37 @@
for ( int i = 0; i < 6 ; i++ )
{
sum = sum + prob[i];
- }
-
+ }
+ //for(int i=0; i<15; i++) {
+ // System.out.println("input index "+i+" input value "+input[i]);
+ //}
+ probs[0] = prob[0]/sum;
+ probs[1] = prob[1]/sum;
+ probs[2] = prob[2]/sum;
+ probs[3] = Math.max(prob[3]/sum,prob[4]/sum);
double max = 0;
int nmax = 0;
- for ( int i = 0; i<6 ; i++ )
+ for ( int i = 0; i<4 ; i++ )
{
- if ( prob[i] > max )
+ if ( probs[i] > max )
{
- max = prob[i];
+ max = probs[i];
nmax = i;
}
}
max = max/sum;
- probs[0] = prob[0]/sum;
- probs[1] = prob[1]/sum;
- probs[2] = prob[2]/sum;
- probs[3] = Math.max(prob[3]/sum,prob[4]/sum);
- System.out.println(probs[0]+" "+probs[1]+" "+probs[2]+" "+probs[3]);
+ if(phohist) aida.cloud1D("Maximum probability type").fill(nmax);
+ if(phohist) aida.cloud1D("ClusterID Photon Prob").fill(probs[0]);
+ if(phohist) aida.cloud1D("ClusterID ChPion Prob").fill(probs[1]);
+ if(phohist) aida.cloud1D("ClusterID Neutron Prob").fill(probs[2]);
+ if(phohist) aida.cloud1D("ClusterID Fragment Prob").fill(probs[3]);
+ //System.out.println(probs[0]+" "+probs[1]+" "+probs[2]+" "+probs[3]);
+ if (nmax==0)
+ {
+ phoclusters.add(ebcluster);
+ badphoton = false;
+ //System.out.println("Real Photon Candidate "+energycor);
+ }
}
}
@@ -901,7 +915,11 @@
} // end loop over photon candidates
if(phohist) aida.cloud1D("number of found photon candidates (above hits threshold)").fill(nGoodClusters);
- if (phoclusters.size() > 0) event.put("PhoBClusters",phoclusters);
+ if(phohist) aida.cloud1D("Largest Energy Photon Cluster").fill(pho1);
+ if(phohist) aida.cloud1D("2nd Largest Energy Photon Cluster").fill(pho2);
+ if(phohist) aida.cloud1D("Cluster2 Divided by Cluster1").fill(pho2/pho1);
+ if(phohist) aida.cloud2D("Cluster2 vs Cluster1").fill(pho1,pho2);
+ if(phoclusters.size() > 0) event.put("PhoBClusters",phoclusters);
double PhoClusE = 0;
if(phohist) aida.cloud1D("Num Photon B Clusters").fill(phoclusters.size());
@@ -1106,6 +1124,8 @@
{
npho++;
phoE += mcpE;
+ //System.out.println("MC Photon "+mcpE);
+
}
// if neutron or Klong, count up
if (pid == 2112 | pid==130 || pid==310 || pid==311)
@@ -1132,6 +1152,8 @@
if(phohist) aida.cloud1D("Num MC Photons").fill(npho);
if(phohist) aida.cloud1D("True Photon E").fill(phoE);
if(phohist) aida.cloud1D("Total Photon Cluster E - True Photon E").fill(PhoClusE-phoE);
+ //System.out.println("Real Photon E "+PhoClusE+" MC Photon E "+phoE);
+
if(phohist) aida.cloud1D("Num True Photons - Photon Clusters").fill(npho-phoclusters.size());
if(trhist) aida.cloud1D("Num Neutrals").fill(nneu);
if(trhist) aida.cloud1D("Neutral E Sum").fill(neuE);