Commit in java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal on MAIN
EcalClusterIC.java+94-461133 -> 1134
HPSEcalClusterIC.java+81-161133 -> 1134
+175-62
2 modified files
update correction functions for energy and position

java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal
EcalClusterIC.java 1133 -> 1134
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/EcalClusterIC.java	2014-10-02 15:48:16 UTC (rev 1133)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/EcalClusterIC.java	2014-10-03 20:25:15 UTC (rev 1134)
@@ -1,5 +1,4 @@
 package org.hps.recon.ecal;
-
 import hep.physics.vec.BasicHep3Vector;
 import hep.physics.vec.Hep3Vector;
 import hep.physics.vec.VecOp;
@@ -16,6 +15,7 @@
 import java.util.Map;
 import java.util.Set;
 
+import org.hps.recon.ecal.HPSEcalClusterIC;
 import org.lcsim.detector.IGeometryInfo;
 import org.lcsim.detector.solids.Trd;
 import org.lcsim.event.CalorimeterHit;
@@ -212,10 +212,10 @@
         	hitList.add(r);
         }
         
-        // Create a list to store the newly created clusters in.
+        //Create a list to store the newly created clusters in.
         ArrayList<HPSEcalClusterIC> clusterList = new ArrayList<HPSEcalClusterIC>();
         
-        // Create a list to store the rejected hits in.
+        //Create a list to store the rejected hits in.
         ArrayList<CalorimeterHit> rejectedHitList = new ArrayList<CalorimeterHit>();
         
         // Sort the list of hits by energy.
@@ -239,18 +239,17 @@
         	else { continue; }
         }
         
-    	// Create a map to connect the cell ID of a calorimeter crystal
-        // to the hit which occurred in that crystal.
+    	//Create a map to connect the cell ID of a calorimeter crystal to the hit which occurred in that crystal.
     	HashMap<Long, CalorimeterHit> hitMap = new HashMap<Long, CalorimeterHit>();
         for (CalorimeterHit hit : hitList) { hitMap.put(hit.getCellID(), hit); }
         
-        // Map a crystal to a list of all clusters in which it is a member.
+        //Map a crystal to a list of all clusters in which it is a member.
         Map<CalorimeterHit, List<CalorimeterHit>> commonHits = new HashMap<CalorimeterHit, List<CalorimeterHit>>();
-        
-        // Map a crystal to the seed of the cluster of which it is a member.
+
+        //Map a crystal to the seed of the cluster of which it is a member.
         HashMap<CalorimeterHit, CalorimeterHit> hitSeedMap = new HashMap<CalorimeterHit, CalorimeterHit>();
         
-      	// Set containing hits immediately around a seed hit.
+      	//Set containing hits immediately around a seed hit.
       	HashSet<CalorimeterHit> surrSeedSet = new HashSet<CalorimeterHit>();
         
         // Loop through all calorimeter hits to locate seeds and perform
@@ -478,8 +477,8 @@
         	CalorimeterHit commonCell = entry1.getKey();
         	CalorimeterHit seedA = entry1.getValue().get(0);
         	CalorimeterHit seedB = entry1.getValue().get(1);    	
-        	double eFractionA = seedEnergy.get(seedA)/(seedEnergy.get(seedA)+seedEnergy.get(seedB));
-        	double eFractionB = seedEnergy.get(seedB)/(seedEnergy.get(seedA)+seedEnergy.get(seedB));
+        	double eFractionA = Math.log(seedEnergy.get(seedA))/Math.log((seedEnergy.get(seedA)+seedEnergy.get(seedB)));
+        	double eFractionB = Math.log(seedEnergy.get(seedB))/Math.log((seedEnergy.get(seedA)+seedEnergy.get(seedB)));
         	double currEnergyA = seedEnergyTot.get(seedA);
         	double currEnergyB = seedEnergyTot.get(seedB);
         	currEnergyA += eFractionA * commonCell.getCorrectedEnergy();
@@ -497,22 +496,19 @@
             // Iterate through known clusters with energies and apply correction.
             for (Map.Entry<CalorimeterHit, Double> entryC : seedEnergyTot.entrySet()) {
                 double rawEnergy = entryC.getValue();
-                // Energy correction for electron:
-                double corrEnergy = rawEnergy / (-0.0027 * rawEnergy - 0.06 / (Math.sqrt(rawEnergy)) + 0.95);
-
-                // Energy correction for positron:
-//                double corrEnergy = rawEnergy / (-0.0096 * rawEnergy - 0.042 / (Math.sqrt(rawEnergy)) + 0.94);
-                // Energy correction for photon:
-//                double corrEnergy = rawEnergy / (0.0015 * rawEnergy - 0.047 / (Math.sqrt(rawEnergy)) + 0.94);
-
                 
+                // Energy correction for initial guess of electron:
+                int pdg = 11;
+                double corrEnergy = enCorrection(pdg, rawEnergy);
+              
                 seedEnergyCorr.put(entryC.getKey(), corrEnergy);    
             }// end of energy corrections
         
                 
         // Cluster Position as per HPS Note 2014-001
         // Create map with seed as key to position/centroid value
-        Map<CalorimeterHit, double[]> seedPosition = new HashMap<CalorimeterHit, double[]>();
+        Map<CalorimeterHit, double[]> rawSeedPosition = new HashMap<CalorimeterHit, double[]>();
+        Map<CalorimeterHit, double[]> corrSeedPosition = new HashMap<CalorimeterHit, double[]>();
         
         // top level iterates through seeds
         for (Map.Entry<CalorimeterHit, Double> entryS : seedEnergyTot.entrySet()) {
@@ -570,31 +566,28 @@
         	xCl = eNumX/eDen;
             yCl = eNumY/eDen;
             
+            double[] rawPosition = new double[3];
+            rawPosition[0] = xCl;
+            rawPosition[1] = yCl;
+            rawPosition[2] = correctedPositionMap.get(seedP)[2];
+            
             // Apply position correction factors:
             // Position correction for electron:
-            double xCorr = xCl-(0.0066/Math.sqrt(seedEnergyTot.get(seedP))-0.03)*xCl-
-            		(0.028*seedEnergyTot.get(seedP)-0.45/Math.sqrt(seedEnergyTot.get(seedP))+0.465);
-            
-            // Position correction for positron:
-//            double xCorr = xCl-(0.0072/Math.sqrt(seedEnergyTot.get(seedP))-0.031)*xCl-
-//            		(0.007*seedEnergyTot.get(seedP)+0.342/Math.sqrt(seedEnergyTot.get(seedP))+0.108);
-            
-            // Position correction for photon:
-//            double xCorr = xCl-(0.005/Math.sqrt(seedEnergyTot.get(seedP))-0.032)*xCl-
-//            		(0.011*seedEnergyTot.get(seedP)-0.037/Math.sqrt(seedEnergyTot.get(seedP))+0.294);
-            
-            
-            
+            int pdg = 11;
+            double xCorr = posCorrection(pdg, xCl, seedEnergyTot.get(seedP));
+           
             double[] corrPosition = new double[2];
             corrPosition[0] = xCorr;
             corrPosition[1] = yCl;
-            seedPosition.put(seedP, corrPosition);
+            corrPosition[2] = correctedPositionMap.get(seedP)[2];
+                        
+            corrSeedPosition.put(seedP, corrPosition);
+            rawSeedPosition.put(seedP, rawPosition);
+
         	
         }// end of cluster position calculation
 
-        
-              
-        
+                
         /*
          * Outputs results to cluster collection. 
          */
@@ -629,16 +622,20 @@
                 		
                     for (Map.Entry<CalorimeterHit, List<CalorimeterHit>> entry4 : commonHits.entrySet()) {
                         if (entry4.getValue().contains(entry2.getKey())) {                       	
-                        	// Added in shared hits for energy distribution between clusters, changed by HS 02JUN14
+                        	// Add shared hits for energy distribution between clusters
                             cluster.addSharedHit(entry4.getKey()); 
                         }
                     }
                                         
-                    //Overwrites cluster energy with corrected cluster energy
-            		if (seedEnergyCorr.values().size() > 0){cluster.setEnergy(seedEnergyCorr.get(entry2.getKey()));}
+                    //Input both raw and corrected cluster energies
+            		if (seedEnergyCorr.values().size() > 0){
+            			cluster.setEnergy(seedEnergyCorr.get(entry2.getKey()));
+            			cluster.setRawEnergy(seedEnergyTot.get(entry2.getKey()));
+            			}
 
-            		//Overwrites cluster position with corrected cluster position
-            		cluster.setPosition(seedPosition.get(entry2.getKey()));
+            		//Input both uncorrected and corrected cluster positions. 
+            		cluster.setCorrPosition(corrSeedPosition.get(entry2.getKey()));
+            		cluster.setRawPosition(rawSeedPosition.get(entry2.getKey()));
                   
                     
                 	}// End checking thresholds and write out.
@@ -647,7 +644,7 @@
 //            System.out.println("Number of clusters: "+clusterList.size());    
 
             
-        } //End event display out loop.
+        } //End event output loop.
         int flag = 1 << LCIOConstants.CLBIT_HITS;
         event.put(clusterCollectionName, clusterList, HPSEcalClusterIC.class, flag);
         event.put(rejectedHitName, rejectedHitList, CalorimeterHit.class, flag);
@@ -730,7 +727,12 @@
      
     
 
-    // Handles pathological case where multiple neighboring crystals have EXACTLY the same energy.
+    /**
+     * Handles pathological case where multiple neighboring crystals have EXACTLY the same energy.
+     * @param hit
+     * @param neighbor Neighbor to hit
+     * @return boolean value of if the hit is a seed
+     */
     private boolean equalEnergies(CalorimeterHit hit, CalorimeterHit neighbor){
     	boolean isSeed = true;
     	
@@ -751,8 +753,54 @@
     	}
     	return isSeed;	
     }
+    /**
+     * Calculates energy correction based on cluster raw energy and particle type as per HPS Note 2014-001
+     * @param pdg Particle id as per PDG
+     * @param rawEnergy Raw Energy of the cluster (sum of hits with shared hit distribution)
+     * @return Corrected Energy
+     */    
+    public double enCorrection(int pdg, double rawEnergy){
+  	   if (pdg == 11) { // Particle is electron
+  		   double corrEnergy = rawEnergy / (-0.0027 * rawEnergy - 0.06 / (Math.sqrt(rawEnergy)) + 0.95);
+  		  return corrEnergy;}
+  	   else if (pdg == -11) { //Particle is positron
+  		   double corrEnergy = rawEnergy / (-0.0096 * rawEnergy - 0.042 / (Math.sqrt(rawEnergy)) + 0.94);
+  		  return corrEnergy;}
+  	   else if (pdg == 22) { //Particle is photon
+  		   double corrEnergy = rawEnergy / (0.0015 * rawEnergy - 0.047 / (Math.sqrt(rawEnergy)) + 0.94);
+  		  return corrEnergy;}
+  	   else { //Unknown 
+  		   double corrEnergy = rawEnergy;
+  		  return corrEnergy;}
+  	   
+     }   
     
+    /**
+     * Calculates position correction based on cluster raw energy, x calculated position, 
+     * and particle type as per HPS Note 2014-001
+     * @param pdg Particle id as per PDG
+     * @param xCl Calculated x centroid position of the cluster, uncorrected, at face
+     * @param rawEnergy Raw energy of the cluster (sum of hits with shared hit distribution)
+     * @return Corrected x position
+     */
+    public double posCorrection(int pdg, double xCl, double rawEnergy){
+    	if (pdg == 11) { //Particle is electron
+    		double xCorr = xCl-(0.0066/Math.sqrt(rawEnergy)-0.03)*xCl-
+    				(0.028*rawEnergy-0.45/Math.sqrt(rawEnergy)+0.465);
+    		return xCorr;}
+    	else if (pdg == -11) {// Particle is positron
+    		double xCorr = xCl-(0.0072/Math.sqrt(rawEnergy)-0.031)*xCl-
+    				(0.007*rawEnergy+0.342/Math.sqrt(rawEnergy)+0.108);
+    		return xCorr;}
+    	else if (pdg == 22) {// Particle is photon
+    		double xCorr = xCl-(0.005/Math.sqrt(rawEnergy)-0.032)*xCl-
+    				(0.011*rawEnergy-0.037/Math.sqrt(rawEnergy)+0.294);
+    		return xCorr;}
+    	else { //Unknown 
+    		double xCorr = xCl;
+    		return xCorr;}
+    	}
     
     
-
-}    
\ No newline at end of file
+    	
+ }    
\ No newline at end of file

java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal
HPSEcalClusterIC.java 1133 -> 1134
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/HPSEcalClusterIC.java	2014-10-02 15:48:16 UTC (rev 1133)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/HPSEcalClusterIC.java	2014-10-03 20:25:15 UTC (rev 1134)
@@ -1,14 +1,8 @@
 package org.hps.recon.ecal;
 
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.Hep3Vector;
-import hep.physics.vec.VecOp;
-
 import java.util.ArrayList;
 import java.util.List;
 
-import org.lcsim.detector.IGeometryInfo;
-import org.lcsim.detector.solids.Trd;
 import org.lcsim.event.CalorimeterHit;
 import org.lcsim.event.base.BaseCluster;
 
@@ -25,7 +19,7 @@
     private CalorimeterHit seedHit = null;
     private long cellID;
     private ArrayList<CalorimeterHit> sharedHitList = new ArrayList<CalorimeterHit>(); 
-    private double[] position = new double[2];
+    private double[] rawPosition = new double[2];
 
     
     
@@ -56,26 +50,97 @@
         }
         return seedHit;
     }
-    
+    /**
+     * Input shared hits between two clusters. 
+     */
     public void addSharedHit(CalorimeterHit sharedHit) {
     	sharedHitList.add(sharedHit);
     }
-    
+    /**
+     * Return shared hit list between two clusters. 
+     */
     public List<CalorimeterHit> getSharedHits() {
     	return sharedHitList;
+    }  
+    /**
+     * Inputs the uncorrected x,y,z position of the cluster.
+     */
+    public void setRawPosition(double[] Position) {
+    	rawPosition = Position;
+    }  
+    /**
+     * Returns the uncorrected x,y,z position of the cluster.
+     */
+    @Override
+    public double[] getPosition(){
+    	return this.rawPosition;
+    }   
+    /**
+     * Do an external calculation of the raw energy and set it. Includes shared hit distribution.
+     */
+    public void setRawEnergy(double rawEnergy){
+    	raw_energy = rawEnergy;
     }
-    
-    
-    public void setPosition(double[] Position) {
+    /**
+     * Inputs the corrected position of the cluster, see HPS Note 2014-001.
+     */
+    public void setCorrPosition(double[] Position) {
     	position = Position;
-    }
-    
-    @Override
-    public double[] getPosition(){
+    }    
+    /**
+     * Returns the corrected position of the cluster. 
+     */
+    public double[] getCorrPosition(){
     	return this.position;
     }
     
+    /**
+     * Calculates energy correction based on cluster raw energy and particle type as per HPS Note 2014-001
+     * @param pdg Particle id as per PDG
+     * @param rawEnergy Raw Energy of the cluster (sum of hits with shared hit distribution)
+     * @return Corrected Energy
+     */    
+    public double energyCorrection(int pdg, double rawEnergy){
+  	   if (pdg == 11) { // Particle is electron
+  		   double corrEnergy = rawEnergy / (-0.0027 * rawEnergy - 0.06 / (Math.sqrt(rawEnergy)) + 0.95);
+  		  return corrEnergy;}
+  	   else if (pdg == -11) { //Particle is positron
+  		   double corrEnergy = rawEnergy / (-0.0096 * rawEnergy - 0.042 / (Math.sqrt(rawEnergy)) + 0.94);
+  		  return corrEnergy;}
+  	   else if (pdg == 22) { //Particle is photon
+  		   double corrEnergy = rawEnergy / (0.0015 * rawEnergy - 0.047 / (Math.sqrt(rawEnergy)) + 0.94);
+  		  return corrEnergy;}
+  	   else { //Unknown 
+  		   double corrEnergy = rawEnergy;
+  		  return corrEnergy;}
+  	   
+     }   
     
+    /**
+     * Calculates position correction based on cluster raw energy, x calculated position, 
+     * and particle type as per HPS Note 2014-001
+     * @param pdg Particle id as per PDG
+     * @param xCl Calculated x centroid position of the cluster, uncorrected, at face
+     * @param rawEnergy Raw energy of the cluster (sum of hits with shared hit distribution)
+     * @return Corrected x position
+     */
+    public double positionCorrection(int pdg, double xCl, double rawEnergy){
+    	if (pdg == 11) { //Particle is electron
+    		double xCorr = xCl-(0.0066/Math.sqrt(rawEnergy)-0.03)*xCl-
+    				(0.028*rawEnergy-0.45/Math.sqrt(rawEnergy)+0.465);
+    		return xCorr;}
+    	else if (pdg == -11) {// Particle is positron
+    		double xCorr = xCl-(0.0072/Math.sqrt(rawEnergy)-0.031)*xCl-
+    				(0.007*rawEnergy+0.342/Math.sqrt(rawEnergy)+0.108);
+    		return xCorr;}
+    	else if (pdg == 22) {// Particle is photon
+    		double xCorr = xCl-(0.005/Math.sqrt(rawEnergy)-0.032)*xCl-
+    				(0.011*rawEnergy-0.037/Math.sqrt(rawEnergy)+0.294);
+    		return xCorr;}
+    	else { //Unknown 
+    		double xCorr = xCl;
+    		return xCorr;}
+    	}
     
     
  /*   @Override
SVNspam 0.1