lcsim-cal-calib/src/org/lcsim/cal/calib
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