lcsim/src/org/lcsim/contrib/uiowa
diff -N ModifiedDetailedNeutralHadronClusterEnergyCalculator.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ModifiedDetailedNeutralHadronClusterEnergyCalculator.java 16 May 2007 20:21:54 -0000 1.1
@@ -0,0 +1,525 @@
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.conditions.*;
+
+/** This duplicates most of the code in
+ org.lcsim.recon.cluster.util.DetailedNeutralHadronClusterEnergyCalculator
+ and is not a good idea in the long term -- just a temporary fix. */
+
+public class ModifiedDetailedNeutralHadronClusterEnergyCalculator implements ClusterEnergyCalculator
+{
+ private ConditionsManager _mgr;
+ double Emeas;
+// mean value of the measured Energy per correction bin
+ double[] fE;
+// correction value for each bin
+ double[] fC;
+// Overall normalization
+ double fN;
+// Energy of the calibration points for each detector
+ double[] EM1BarEnergy;
+ double[] EM1ECEnergy;
+ double[] EM2BarEnergy;
+ double[] EM2ECEnergy;
+ double[] HADBarEnergy;
+ double[] HADECEnergy;
+// Mean responses for each detector - need both analog and digital for EM
+ double[] EM1BarAnaResponse;
+ double[] EM1BarDigResponse;
+ double[] EM1ECAnaResponse;
+ double[] EM1ECDigResponse;
+ double[] EM2BarAnaResponse;
+ double[] EM2BarDigResponse;
+ double[] EM2ECAnaResponse;
+ double[] EM2ECDigResponse;
+ double[] HADBarResponse;
+ double[] HADECResponse;
+// Angle corrections for each detector
+ double[] EM1BarAnaAngCorr;
+ double[] EM1BarDigAngCorr;
+ double[] EM1ECAnaAngCorr;
+ double[] EM1ECDigAngCorr;
+ double[] EM2BarAnaAngCorr;
+ double[] EM2BarDigAngCorr;
+ double[] EM2ECAnaAngCorr;
+ double[] EM2ECDigAngCorr;
+ double[] HADBarAngCorr;
+ double[] HADECAngCorr;
+// hit energy cuts
+ double EM1BarEcut;
+ double EM1ECEcut;
+ double EM2BarEcut;
+ double EM2ECEcut;
+ double HADBarEcut;
+ double HADECEcut;
+// hit time cuts
+ double EM1BarTcut;
+ double EM1ECTcut;
+ double EM2BarTcut;
+ double EM2ECTcut;
+ double HADBarTcut;
+ double HADECTcut;
+// EM #hit cut to switch from digital to analog
+ double EM1BarDtoAcut;
+ double EM1ECDtoAcut;
+ double EM2BarDtoAcut;
+ double EM2ECDtoAcut;
+ public ModifiedDetailedNeutralHadronClusterEnergyCalculator()
+ {
+ }
+ private void init() {
+//
+// Get the calibration constants
+//
+ _mgr = ConditionsManager.defaultInstance();
+ ConditionsSet cs = _mgr.getConditions("hadronCalibration/HADBarrel");
+ HADBarEcut = cs.getDouble("hitEnergycut");
+ HADBarTcut = cs.getDouble("hitTimecut");
+ HADBarEnergy = cs.getDoubleArray("calEnergies");
+ HADBarResponse = cs.getDoubleArray("calResponses");
+ HADBarAngCorr = cs.getDoubleArray("calAngleCorrections");
+ cs = _mgr.getConditions("hadronCalibration/HADEndcap");
+ HADECEcut = cs.getDouble("hitEnergycut");
+ HADECTcut = cs.getDouble("hitTimecut");
+ HADECEnergy = cs.getDoubleArray("calEnergies");
+ HADECResponse = cs.getDoubleArray("calResponses");
+ HADECAngCorr = cs.getDoubleArray("calAngleCorrections");
+ cs = _mgr.getConditions("hadronCalibration/EMBarrel");
+ EM1BarEcut = cs.getDouble("hitEnergycut1");
+ EM1BarTcut = cs.getDouble("hitTimecut1");
+ EM1BarDtoAcut = cs.getDouble("hitDtoAcut1");
+ EM1BarEnergy = cs.getDoubleArray("calEnergies1");
+ EM1BarAnaResponse = cs.getDoubleArray("calAnaResponses1");
+ EM1BarAnaAngCorr = cs.getDoubleArray("calAnaAngleCorrections1");
+ EM1BarDigResponse = cs.getDoubleArray("calDigResponses1");
+ EM1BarDigAngCorr = cs.getDoubleArray("calDigAngleCorrections1");
+ EM2BarEcut = cs.getDouble("hitEnergycut2");
+ EM2BarTcut = cs.getDouble("hitTimecut2");
+ EM2BarDtoAcut = cs.getDouble("hitDtoAcut2");
+ EM2BarEnergy = cs.getDoubleArray("calEnergies2");
+ EM2BarAnaResponse = cs.getDoubleArray("calAnaResponses2");
+ EM2BarAnaAngCorr = cs.getDoubleArray("calAnaAngleCorrections2");
+ EM2BarDigResponse = cs.getDoubleArray("calDigResponses2");
+ EM2BarDigAngCorr = cs.getDoubleArray("calDigAngleCorrections2");
+ cs = _mgr.getConditions("hadronCalibration/EMEndcap");
+ EM1ECEcut = cs.getDouble("hitEnergycut1");
+ EM1ECTcut = cs.getDouble("hitTimecut1");
+ EM1ECDtoAcut = cs.getDouble("hitDtoAcut1");
+ EM1ECEnergy = cs.getDoubleArray("calEnergies1");
+ EM1ECAnaResponse = cs.getDoubleArray("calAnaResponses1");
+ EM1ECAnaAngCorr = cs.getDoubleArray("calAnaAngleCorrections1");
+ EM1ECDigResponse = cs.getDoubleArray("calDigResponses1");
+ EM1ECDigAngCorr = cs.getDoubleArray("calDigAngleCorrections1");
+ EM2ECEcut = cs.getDouble("hitEnergycut2");
+ EM2ECTcut = cs.getDouble("hitTimecut2");
+ EM2ECDtoAcut = cs.getDouble("hitDtoAcut2");
+ EM2ECEnergy = cs.getDoubleArray("calEnergies2");
+ EM2ECAnaResponse = cs.getDoubleArray("calAnaResponses2");
+ EM2ECAnaAngCorr = cs.getDoubleArray("calAnaAngleCorrections2");
+ EM2ECDigResponse = cs.getDoubleArray("calDigResponses2");
+ EM2ECDigAngCorr = cs.getDoubleArray("calDigAngleCorrections2");
+ cs = _mgr.getConditions("hadronCalibration/ZpoleInversion");
+ fN = cs.getDouble("normalization");
+ fE = cs.getDoubleArray("meanEnergy");
+ fC = cs.getDoubleArray("corrections");
+
+ }
+ public double getEnergy(Cluster c)
+ {
+ if (_mgr==null) {
+ init();
+ }
+
+ double EmeasEst = 0.;
+//
+// First get the sums for each subdetector
+//
+ double[] esums = {0.,0.,0.,0.,0.,0.};
+ int[] hitsums = {0,0,0,0,0,0};
+ for(CalorimeterHit hit:c.getCalorimeterHits())
+ {
+ 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)
+ {
+ hitsums[index]++;
+ esums[index] += 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(hitsums[0] > 0)
+ {
+ ndetcont++;
+ if(hitsums[0] < EM1BarDtoAcut)
+ {
+ EmeasEst += EcontEst( (double) hitsums[0],EM1BarEnergy,EM1BarDigResponse,EM1BarDigAngCorr,stheta);
+ }
+ else
+ {
+ EmeasEst += EcontEst( esums[0],EM1BarEnergy,EM1BarAnaResponse,EM1BarAnaAngCorr,stheta);
+ }
+ }
+ if(hitsums[1] > 0)
+ {
+ ndetcont++;
+ if(hitsums[1] < EM2BarDtoAcut)
+ {
+ EmeasEst += EcontEst( (double) hitsums[1],EM2BarEnergy,EM2BarDigResponse,EM2BarDigAngCorr,stheta);
+ }
+ else
+ {
+ EmeasEst += EcontEst( esums[1],EM2BarEnergy,EM2BarAnaResponse,EM2BarAnaAngCorr,stheta);
+ }
+ }
+ if(hitsums[2] > 0)
+ {
+ ndetcont++;
+ if(hitsums[2] < EM1ECDtoAcut)
+ {
+ EmeasEst += EcontEst( (double) hitsums[2],EM1ECEnergy,EM1ECDigResponse,EM1ECDigAngCorr,ctheta);
+ }
+ else
+ {
+ EmeasEst += EcontEst( esums[2],EM1ECEnergy,EM1ECAnaResponse,EM1ECAnaAngCorr,ctheta);
+ }
+ }
+ if(hitsums[3] > 0)
+ {
+ ndetcont++;
+ if(hitsums[3] < EM2ECDtoAcut)
+ {
+ EmeasEst += EcontEst( (double) hitsums[3],EM2ECEnergy,EM2ECDigResponse,EM2ECDigAngCorr,ctheta);
+ }
+ else
+ {
+ EmeasEst += EcontEst( esums[3],EM2ECEnergy,EM2ECAnaResponse,EM2ECAnaAngCorr,ctheta);
+ }
+ }
+ if(hitsums[4] > 0)
+ {
+ ndetcont++;
+ EmeasEst += EcontEst( (double) hitsums[4],HADBarEnergy,HADBarResponse,HADBarAngCorr,stheta);
+ }
+ if(hitsums[5] > 0)
+ {
+ ndetcont++;
+ EmeasEst += EcontEst( (double) hitsums[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(hitsums[0] > 0)
+ {
+ if(hitsums[0] < EM1BarDtoAcut)
+ {
+ Emeas += Econt( (double) hitsums[0],EM1BarEnergy,EM1BarDigResponse,EM1BarDigAngCorr,stheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( esums[0],EM1BarEnergy,EM1BarAnaResponse,EM1BarAnaAngCorr,stheta,EmeasEst);
+ }
+ }
+ if(hitsums[1] > 0)
+ {
+ if(hitsums[1] < EM2BarDtoAcut)
+ {
+ Emeas += Econt( (double) hitsums[1],EM2BarEnergy,EM2BarDigResponse,EM2BarDigAngCorr,stheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( esums[1],EM2BarEnergy,EM2BarAnaResponse,EM2BarAnaAngCorr,stheta,EmeasEst);
+ }
+ }
+ if(hitsums[2] > 0)
+ {
+ if(hitsums[2] < EM1ECDtoAcut)
+ {
+ Emeas += Econt( (double) hitsums[2],EM1ECEnergy,EM1ECDigResponse,EM1ECDigAngCorr,ctheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( esums[2],EM1ECEnergy,EM1ECAnaResponse,EM1ECAnaAngCorr,ctheta,EmeasEst);
+ }
+ }
+ if(hitsums[3] > 0)
+ {
+ if(hitsums[3] < EM2ECDtoAcut)
+ {
+ Emeas += Econt( (double) hitsums[3],EM2ECEnergy,EM2ECDigResponse,EM2ECDigAngCorr,ctheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( esums[3],EM2ECEnergy,EM2ECAnaResponse,EM2ECAnaAngCorr,ctheta,EmeasEst);
+ }
+ }
+ if(hitsums[4] > 0)
+ {
+ Emeas += Econt( (double) hitsums[4],HADBarEnergy,HADBarResponse,HADBarAngCorr,stheta,EmeasEst);
+ }
+ if(hitsums[5] > 0)
+ {
+ Emeas += Econt( (double) hitsums[5],HADECEnergy,HADECResponse,HADECAngCorr,ctheta,EmeasEst);
+ }
+//
+// Iterate on the energy estimate
+//
+ EmeasEst = Emeas;
+ Emeas = 0.;
+ if(hitsums[0] > 0)
+ {
+ if(hitsums[0] < EM1BarDtoAcut)
+ {
+ Emeas += Econt( (double) hitsums[0],EM1BarEnergy,EM1BarDigResponse,EM1BarDigAngCorr,stheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( esums[0],EM1BarEnergy,EM1BarAnaResponse,EM1BarAnaAngCorr,stheta,EmeasEst);
+ }
+ }
+ if(hitsums[1] > 0)
+ {
+ if(hitsums[1] < EM2BarDtoAcut)
+ {
+ Emeas += Econt( (double) hitsums[1],EM2BarEnergy,EM2BarDigResponse,EM2BarDigAngCorr,stheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( esums[1],EM2BarEnergy,EM2BarAnaResponse,EM2BarAnaAngCorr,stheta,EmeasEst);
+ }
+ }
+ if(hitsums[2] > 0)
+ {
+ if(hitsums[2] < EM1ECDtoAcut)
+ {
+ Emeas += Econt( (double) hitsums[2],EM1ECEnergy,EM1ECDigResponse,EM1ECDigAngCorr,ctheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( esums[2],EM1ECEnergy,EM1ECAnaResponse,EM1ECAnaAngCorr,ctheta,EmeasEst);
+ }
+ }
+ if(hitsums[3] > 0)
+ {
+ if(hitsums[3] < EM2ECDtoAcut)
+ {
+ Emeas += Econt( (double) hitsums[3],EM2ECEnergy,EM2ECDigResponse,EM2ECDigAngCorr,ctheta,EmeasEst);
+ }
+ else
+ {
+ Emeas += Econt( esums[3],EM2ECEnergy,EM2ECAnaResponse,EM2ECAnaAngCorr,ctheta,EmeasEst);
+ }
+ }
+ if(hitsums[4] > 0)
+ {
+ Emeas += Econt( (double) hitsums[4],HADBarEnergy,HADBarResponse,HADBarAngCorr,stheta,EmeasEst);
+ }
+ if(hitsums[5] > 0)
+ {
+ Emeas += Econt( (double) hitsums[5],HADECEnergy,HADECResponse,HADECAngCorr,ctheta,EmeasEst);
+ }
+ }
+//
+// Energy estimate done, now use inversion calibration data to get best value of
+// neutral energy from energy estimate
+//
+ double Eret = 0.;
+ if (m_doInversion) {
+ int nC = fE.length;
+ 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;
+ }
+ }
+ } else {
+ Eret = Emeas;
+ }
+//
+// Put a lower bound of (default: 1 GeV) on Energy and return
+//
+ if (Eret > m_eMin) {
+ return Eret;
+ } else {
+ return m_eMin;
+ }
+ }
+ 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;
+ }
+ }
+
+ // Mat makes it more flexible
+ protected double m_eMin = 1.0;
+ public void setMinimumEnergy(double eMin) {
+ m_eMin = eMin;
+ }
+ protected boolean m_doInversion = true;
+ public void setDoInversion(boolean doInversion) {
+ m_doInversion = doInversion;
+ }
+}
