Print

Print


Commit in lcsim-cal-calib/src/org/lcsim/cal/calib on MAIN
ClusterEnergyAnalysis.java+61-191.2 -> 1.3
organize ths histogram tree

lcsim-cal-calib/src/org/lcsim/cal/calib
ClusterEnergyAnalysis.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- ClusterEnergyAnalysis.java	15 Jul 2008 01:45:14 -0000	1.2
+++ ClusterEnergyAnalysis.java	15 Jul 2008 01:58:25 -0000	1.3
@@ -3,12 +3,13 @@
  *
  * Created on July 14, 2008, 6:18 PM
  *
- * $Id: ClusterEnergyAnalysis.java,v 1.2 2008/07/15 01:45:14 ngraf Exp $
+ * $Id: ClusterEnergyAnalysis.java,v 1.3 2008/07/15 01:58:25 ngraf Exp $
  */
 
 package org.lcsim.cal.calib;
 
 import hep.aida.ITree;
+import hep.physics.vec.Hep3Vector;
 import java.util.List;
 import org.lcsim.conditions.ConditionsManager;
 import org.lcsim.conditions.ConditionsManager.ConditionsSetNotFoundException;
@@ -28,8 +29,10 @@
 import static java.lang.Math.sin;
 import static java.lang.Math.PI;
 import static java.lang.Math.abs;
+import static java.lang.Math.sqrt;
 import java.util.HashMap;
 import java.util.Map;
+import org.lcsim.event.MCParticle;
 
 /**
  *
@@ -44,6 +47,9 @@
     private double[] _ecalLayering;
     boolean _useFirstLayer;
     
+    private double emCalInnerRadius = 0.;
+    private double emCalInnerZ = 0.;
+    
     private Map<String, Double> _fitParameters = new HashMap<String, Double>();
     
     
@@ -78,7 +84,7 @@
             
             _ecalLayering = _cond.getDoubleArray("ECalLayering");
             _useFirstLayer = _cond.getDouble("IsFirstEmLayerSampling")==1.;
-
+            
             // photons
             String photonFitParametersList = _cond.getString("PhotonFitParameters");
             String[]  photonFitParameters = photonFitParametersList.split(",\\s");
@@ -96,23 +102,59 @@
             
             _initialized = true;
         }
-        String processedHitsName = _cond.getString("ProcessedHitsCollectionName");
-        List<CalorimeterHit> hitsToCluster = _collectionmanager.getList(processedHitsName);
-        // cluster the hits
-        List<Cluster> clusters = _fcc.createClusters(hitsToCluster);
-        String type = "gamma";
-        for(Cluster c : clusters)
+        
+        //start processing the event.
+        // organize the histogram tree by species and energy
+        List<MCParticle> mcparts = event.getMCParticles();
+        MCParticle mcpart = mcparts.get(mcparts.size()-1);
+        String particleType = mcpart.getType().getName();
+        double mcEnergy = mcpart.getEnergy();
+        long mcIntegerEnergy = Math.round(mcEnergy);
+        boolean meV = false;
+        if(mcEnergy<.99)
         {
-            aida.cloud1D("uncorrected cluster energy for all clusters").fill(c.getEnergy());
-            double e = correctClusterEnergy(c, type);
-            aida.cloud1D("corrected cluster energy for all clusters").fill(e);
-            SpacePoint p = new CartesianPoint(c.getPosition());
-            double clusterTheta = p.theta();
-            aida.cloud2D("uncorrected cluster energy for all clusters vs theta").fill(clusterTheta, c.getEnergy());
-            aida.cloud2D("corrected cluster energy for all clusters vs theta").fill(clusterTheta, e);
-            aida.cloud2D("raw-corrected cluster energy for all clusters vs theta").fill(clusterTheta, c.getEnergy()-e);
+            mcIntegerEnergy = Math.round(mcEnergy*1000);
+            meV = true;
         }
         
+        _tree.mkdirs(particleType);
+        _tree.cd(particleType);
+        _tree.mkdirs(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
+        _tree.cd(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
+        
+        // this analysis is intended for single particle calorimeter response.
+        // let's make sure that the primary particle did not interact in the
+        // tracker...
+        Hep3Vector endpoint = mcpart.getEndPoint();
+        // this is just crap. Why not use SpacePoint?
+        double radius = sqrt(endpoint.x()*endpoint.x()+endpoint.y()*endpoint.y());
+        double z = endpoint.z();
+//        System.out.println("Input MCParticle endpoint: r="+radius+" z= "+z);
+        
+        boolean doit = true;
+        if(radius<emCalInnerRadius && abs(z) < emCalInnerZ) doit = false;
+        if(doit)
+        {
+            String processedHitsName = _cond.getString("ProcessedHitsCollectionName");
+            List<CalorimeterHit> hitsToCluster = _collectionmanager.getList(processedHitsName);
+            // cluster the hits
+            List<Cluster> clusters = _fcc.createClusters(hitsToCluster);
+            String type = "gamma";
+            for(Cluster c : clusters)
+            {
+                aida.cloud1D("uncorrected cluster energy for all clusters").fill(c.getEnergy());
+                double e = correctClusterEnergy(c, type);
+                aida.cloud1D("corrected cluster energy for all clusters").fill(e);
+                SpacePoint p = new CartesianPoint(c.getPosition());
+                double clusterTheta = p.theta();
+                aida.cloud2D("uncorrected cluster energy for all clusters vs theta").fill(clusterTheta, c.getEnergy());
+                aida.cloud2D("corrected cluster energy for all clusters vs theta").fill(clusterTheta, e);
+                aida.cloud2D("raw-corrected cluster energy for all clusters vs theta").fill(clusterTheta, c.getEnergy()-e);
+            }
+        }// end of check on decays outside tracker volume
+        //reset our histogram tree
+         _tree.cd("/");
+        
     }
     
     private double correctClusterEnergy(Cluster c, String type)
@@ -175,9 +217,9 @@
             }
             
             double correctionFactor = a + b*sin(normalTheta);
-            aida.cloud2D(name+" "+layer+"hit Theta vs correction factor ").fill(hitTheta, correctionFactor);
-            aida.cloud2D(name+" hit Theta vs correction factor ").fill(hitTheta, correctionFactor);
-            if(isEM) aida.cloud2D("EM layer vs caltype").fill(layer, caltype);
+//            aida.cloud2D(name+" "+layer+"hit Theta vs correction factor ").fill(hitTheta, correctionFactor);
+//            aida.cloud2D(name+" hit Theta vs correction factor ").fill(hitTheta, correctionFactor);
+//            if(isEM) aida.cloud2D("EM layer vs caltype").fill(layer, caltype);
             
             // now apply the correction
             
CVSspam 0.2.8