lcsim/src/org/lcsim/contrib/uiowa
diff -N FuzzyNeutralHadronClusterEnergyCalculator.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FuzzyNeutralHadronClusterEnergyCalculator.java 19 Nov 2007 21:50:06 -0000 1.1
@@ -0,0 +1,421 @@
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.conditions.*;
+import org.lcsim.recon.cluster.util.*;
+
+import org.lcsim.contrib.uiowa.FuzzyCalorimeterHit;
+
+/**
+ * Custom calibration designed for neutral hadrons which allows
+ * for hits to be weighted -- i.e. to be shared among multiple
+ * clusters such that the total weight adds up to 1. This is
+ * an extension of Ron's DetailedNeutralHadronClusterEnergyCalculator
+ * class.
+ *
+ * @version $Id: FuzzyNeutralHadronClusterEnergyCalculator.java,v 1.1 2007/11/19 21:50:06 mcharles Exp $
+ */
+
+public class FuzzyNeutralHadronClusterEnergyCalculator extends MatDetailedNeutralHadronClusterEnergyCalculator
+{
+ /** Simple constructor. */
+ public FuzzyNeutralHadronClusterEnergyCalculator() {
+ super();
+ }
+
+ /** Find the energy of this cluster. */
+ public double getEnergy(Cluster c)
+ {
+ if (_mgr==null) {
+ init();
+ }
+
+ double EmeasEst = 0.;
+//
+// First get the sums for each subdetector
+//
+ double[] unweightedEnergySums = {0.,0.,0.,0.,0.,0.};
+ double[] weightedEnergySums = {0.,0.,0.,0.,0.,0.};
+ int[] unweightedHitSums = {0,0,0,0,0,0};
+ double[] weightedHitSums = {0,0,0,0,0,0};
+ for(CalorimeterHit hit:c.getCalorimeterHits())
+ {
+ double hitWeight = 1.0;
+ if (hit instanceof FuzzyCalorimeterHit) {
+ hitWeight = ((FuzzyCalorimeterHit)(hit)).getWeight();
+ }
+ int index = -1;
+ double timecut = 0.;
+ double Ecut = 0.;
+ IDDecoder idc = hit.getIDDecoder();
+ idc.setID(hit.getCellID());
+ int detector_index = idc.getValue("system");
+ if(detector_index == 2)
+ {
+ int layer = idc.getValue("layer");
+ if(layer < 20)
+ {
+ index = 0;
+ timecut = EM1BarTcut;
+ Ecut = EM1BarEcut;
+ }
+ else
+ {
+ index = 1;
+ timecut = EM2BarTcut;
+ Ecut = EM2BarEcut;
+ }
+ }
+ else if(detector_index == 6)
+ {
+ int layer = idc.getValue("layer");
+ if(layer < 20)
+ {
+ index = 2;
+ timecut = EM1ECTcut;
+ Ecut = EM1ECEcut;
+ }
+ else
+ {
+ index = 3;
+ timecut = EM2ECTcut;
+ Ecut = EM2ECEcut;
+ }
+ }
+ else if(detector_index == 3)
+ {
+ index = 4;
+ timecut = HADBarTcut;
+ Ecut = HADBarEcut;
+ }
+ else if(detector_index == 7)
+ {
+ index = 5;
+ timecut = HADECTcut;
+ Ecut = HADECEcut;
+ }
+ if(index > -1)
+ {
+ if(hit.getTime() < timecut)
+ {
+ if(hit.getRawEnergy() > Ecut)
+ {
+ unweightedHitSums[index]++;
+ weightedHitSums[index] += hitWeight;
+ unweightedEnergySums[index] += hit.getRawEnergy();
+ weightedEnergySums[index] += hitWeight * hit.getRawEnergy();
+ }
+ }
+ }
+ }
+//
+// Use the position to calculate an angle, assuming origin at 0
+//
+ double[] pos = c.getPosition();
+ double R = Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
+ double ctheta = Math.abs(pos[2])/R;
+ double stheta = Math.sqrt(1.-ctheta*ctheta);
+ int ndetcont = 0;
+//
+// Estimate the energy by adding up contributions from each subdetector,
+// using isolated detector response calibration
+//
+ if(weightedHitSums[0] > 0)
+ {
+ ndetcont++;
+ if(weightedHitSums[0] < EM1BarDtoAcut)
+ {
+ EmeasEst += EcontEst( weightedHitSums[0],EM1BarEnergy,EM1BarDigResponse,EM1BarDigAngCorr,stheta);
+ }
+ else
+ {
+ EmeasEst += EcontEst( weightedEnergySums[0],EM1BarEnergy,EM1BarAnaResponse,EM1BarAnaAngCorr,stheta);
+ }
+ }
+ if(weightedHitSums[1] > 0)
+ {
+ ndetcont++;
+ if(weightedHitSums[1] < EM2BarDtoAcut)
+ {
+ EmeasEst += EcontEst( weightedHitSums[1],EM2BarEnergy,EM2BarDigResponse,EM2BarDigAngCorr,stheta);
+ }
+ else
+ {
+ EmeasEst += EcontEst( weightedEnergySums[1],EM2BarEnergy,EM2BarAnaResponse,EM2BarAnaAngCorr,stheta);
+ }
+ }
+ if(weightedHitSums[2] > 0)
+ {
+ ndetcont++;
+ if(weightedHitSums[2] < EM1ECDtoAcut)
+ {
+ EmeasEst += EcontEst( weightedHitSums[2],EM1ECEnergy,EM1ECDigResponse,EM1ECDigAngCorr,ctheta);
+ }
+ else
+ {
+ EmeasEst += EcontEst( weightedEnergySums[2],EM1ECEnergy,EM1ECAnaResponse,EM1ECAnaAngCorr,ctheta);
+ }
+ }
+ if(weightedHitSums[3] > 0)
+ {
+ ndetcont++;
+ if(weightedHitSums[3] < EM2ECDtoAcut)
+ {
+ EmeasEst += EcontEst( weightedHitSums[3],EM2ECEnergy,EM2ECDigResponse,EM2ECDigAngCorr,ctheta);
+ }
+ else
+ {
+ EmeasEst += EcontEst( weightedEnergySums[3],EM2ECEnergy,EM2ECAnaResponse,EM2ECAnaAngCorr,ctheta);
+ }
+ }
+ if(weightedHitSums[4] > 0)
+ {
+ ndetcont++;
+ EmeasEst += EcontEst( weightedHitSums[4],HADBarEnergy,HADBarResponse,HADBarAngCorr,stheta);
+ }
+ if(weightedHitSums[5] > 0)
+ {
+ ndetcont++;
+ EmeasEst += EcontEst( weightedHitSums[5],HADECEnergy,HADECResponse,HADECAngCorr,ctheta);
+ }
+//
+// If only one subdetector contributed, we're done with the estimate
+//
+ if(ndetcont < 2)
+ {
+ Emeas = EmeasEst;
+ }
+//
+// Otherwise, refine the estimate by treating the subdetector contribution as a
+// fraction of the response of a neutral with the energy estimate
+//
+ else
+ {
+ Emeas = 0.;
+ if(weightedHitSums[0] > 0)
+ {
+ if(weightedHitSums[0] < EM1BarDtoAcut)
+ {
+ Emeas += Econt( weightedHitSums[0],EM1BarEnergy,EM1BarDigResponse,EM1BarDigAngCorr,stheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( weightedEnergySums[0],EM1BarEnergy,EM1BarAnaResponse,EM1BarAnaAngCorr,stheta,EmeasEst);
+ }
+ }
+ if(weightedHitSums[1] > 0)
+ {
+ if(weightedHitSums[1] < EM2BarDtoAcut)
+ {
+ Emeas += Econt( weightedHitSums[1],EM2BarEnergy,EM2BarDigResponse,EM2BarDigAngCorr,stheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( weightedEnergySums[1],EM2BarEnergy,EM2BarAnaResponse,EM2BarAnaAngCorr,stheta,EmeasEst);
+ }
+ }
+ if(weightedHitSums[2] > 0)
+ {
+ if(weightedHitSums[2] < EM1ECDtoAcut)
+ {
+ Emeas += Econt( weightedHitSums[2],EM1ECEnergy,EM1ECDigResponse,EM1ECDigAngCorr,ctheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( weightedEnergySums[2],EM1ECEnergy,EM1ECAnaResponse,EM1ECAnaAngCorr,ctheta,EmeasEst);
+ }
+ }
+ if(weightedHitSums[3] > 0)
+ {
+ if(weightedHitSums[3] < EM2ECDtoAcut)
+ {
+ Emeas += Econt( weightedHitSums[3],EM2ECEnergy,EM2ECDigResponse,EM2ECDigAngCorr,ctheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( weightedEnergySums[3],EM2ECEnergy,EM2ECAnaResponse,EM2ECAnaAngCorr,ctheta,EmeasEst);
+ }
+ }
+ if(weightedHitSums[4] > 0)
+ {
+ Emeas += Econt( weightedHitSums[4],HADBarEnergy,HADBarResponse,HADBarAngCorr,stheta,EmeasEst);
+ }
+ if(weightedHitSums[5] > 0)
+ {
+ Emeas += Econt( weightedHitSums[5],HADECEnergy,HADECResponse,HADECAngCorr,ctheta,EmeasEst);
+ }
+//
+// Iterate on the energy estimate
+//
+ EmeasEst = Emeas;
+ Emeas = 0.;
+ if(weightedHitSums[0] > 0)
+ {
+ if(weightedHitSums[0] < EM1BarDtoAcut)
+ {
+ Emeas += Econt( weightedHitSums[0],EM1BarEnergy,EM1BarDigResponse,EM1BarDigAngCorr,stheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( weightedEnergySums[0],EM1BarEnergy,EM1BarAnaResponse,EM1BarAnaAngCorr,stheta,EmeasEst);
+ }
+ }
+ if(weightedHitSums[1] > 0)
+ {
+ if(weightedHitSums[1] < EM2BarDtoAcut)
+ {
+ Emeas += Econt( weightedHitSums[1],EM2BarEnergy,EM2BarDigResponse,EM2BarDigAngCorr,stheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( weightedEnergySums[1],EM2BarEnergy,EM2BarAnaResponse,EM2BarAnaAngCorr,stheta,EmeasEst);
+ }
+ }
+ if(weightedHitSums[2] > 0)
+ {
+ if(weightedHitSums[2] < EM1ECDtoAcut)
+ {
+ Emeas += Econt( weightedHitSums[2],EM1ECEnergy,EM1ECDigResponse,EM1ECDigAngCorr,ctheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( weightedEnergySums[2],EM1ECEnergy,EM1ECAnaResponse,EM1ECAnaAngCorr,ctheta,EmeasEst);
+ }
+ }
+ if(weightedHitSums[3] > 0)
+ {
+ if(weightedHitSums[3] < EM2ECDtoAcut)
+ {
+ Emeas += Econt( weightedHitSums[3],EM2ECEnergy,EM2ECDigResponse,EM2ECDigAngCorr,ctheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( weightedEnergySums[3],EM2ECEnergy,EM2ECAnaResponse,EM2ECAnaAngCorr,ctheta,EmeasEst);
+ }
+ }
+ if(weightedHitSums[4] > 0)
+ {
+ Emeas += Econt( weightedHitSums[4],HADBarEnergy,HADBarResponse,HADBarAngCorr,stheta,EmeasEst);
+ }
+ if(weightedHitSums[5] > 0)
+ {
+ Emeas += Econt( weightedHitSums[5],HADECEnergy,HADECResponse,HADECAngCorr,ctheta,EmeasEst);
+ }
+ }
+//
+// Energy estimate done, now use inversion calibration data to get best value of
+// neutral energy from energy estimate
+//
+ int nC = fE.length;
+ double Eret = 0.;
+ int j = 0;
+ if(Emeas > fE[nC-1])
+ {
+ Eret = Emeas;
+ Eret *= fC[nC-1]*fN;
+ }
+ else
+ {
+ while(Emeas > fE[j])
+ {j++;}
+ if(j == 0)
+ {
+ Eret = Emeas;
+ Eret *= fC[0]*fN;
+ }
+ else
+ {
+ double fc = fC[j-1] + (fC[j] - fC[j-1])*(Emeas - fE[j-1])/(fE[j] - fE[j-1]);
+ Eret = Emeas;
+ Eret *= fc*fN;
+ }
+ }
+//
+// Put a lower bound on Energy and return
+//
+ if (Eret > minimumAllowedEnergy) {
+ return Eret;
+ } else {
+ return minimumAllowedEnergy;
+ }
+ }
+ public double EcontEst(double val,double[] Energy,double[] res,double[] angc,double ang)
+ {
+ double alpha;
+ int nE = Energy.length;
+ int i=0;
+ while((i<nE)&&(val > res[i]))
+ {i++;}
+ if(i==0)alpha = angc[0];
+ else if(i == nE)alpha = angc[nE-1];
+ else alpha = angc[i-1] + (angc[i] - angc[i-1])*(val - res[i-1])/(res[i] - res[i-1]);
+ double st = ang;
+ if(st < .707)st = .707;
+ double corr = 1. + alpha*(st - 1.);
+ double veff = val/corr;
+ i=0;
+ while((i<nE)&&(veff > res[i]))
+ {i++;}
+ double E0;
+ double E1;
+ double v0;
+ double v1;
+ if(i == 0)
+ {
+ E0 = 0.;
+ E1 = Energy[0];
+ v0 = 0;
+ v1 = res[0];
+ }
+ else if(i == nE)
+ {
+ E0 = Energy[nE-2];
+ E1 = Energy[nE-1];
+ v0 = res[nE-2];
+ v1 = res[nE-1];
+ }
+ else
+ {
+ E0 = Energy[i-1];
+ E1 = Energy[i];
+ v0 = res[i-1];
+ v1 = res[i];
+ }
+ return E0 + (E1 - E0)*(veff - v0)/(v1 - v0);
+ }
+ public double Econt(double val,double[] Energy,double[] res,double[] angc,double ang,double Eest)
+ {
+ double alpha;
+ double st = ang;
+ if(st < .707)st = .707;
+ int nE = Energy.length;
+ int i=0;
+ while((i<nE)&&(Eest > Energy[i]))
+ {i++;}
+ if(i==0)
+ {
+ alpha = angc[0];
+ double response = Eest/Energy[0]*res[0]*(1.+alpha*(st - 1.));
+ return val/response*Eest;
+ }
+ else if(i == nE)
+ {
+ alpha = angc[nE-1];
+ double response = (res[nE-2] + (res[nE-1] - res[nE-2])*(Eest - Energy[nE-2])/(Energy[nE-1] - Energy[nE-2]))*(1.+alpha*(st - 1.));
+ return val*Eest/response;
+ }
+ else
+ {
+ alpha = angc[i-1] + (angc[i] - angc[i-1])*(val - res[i-1])/(res[i] - res[i-1]);
+ double response = (res[i-1] + (res[i] - res[i-1])*(Eest - Energy[i-1])/(Energy[i] - Energy[i-1]))*(1.+alpha*(st - 1.));
+ return val*Eest/response;
+ }
+ }
+ /** Specify the minimum energy that will be returned. Default is 1.0 GeV. */
+ public void setMinimumAllowedEnergy(double eMin) {
+ minimumAllowedEnergy = eMin;
+ }
+}