lcsim-cal-calib/src/org/lcsim/cal/calib
diff -u -r1.2 -r1.3
--- SamplingFractionAnalysisPolyCalDriver.java 28 May 2010 18:44:37 -0000 1.2
+++ SamplingFractionAnalysisPolyCalDriver.java 3 Nov 2010 15:19:51 -0000 1.3
@@ -3,9 +3,8 @@
*
* Created on May 19, 2008, 11:54 AM
*
- * $Id: SamplingFractionAnalysisPolyCalDriver.java,v 1.2 2010/05/28 18:44:37 ngraf Exp $
+ * $Id: SamplingFractionAnalysisPolyCalDriver.java,v 1.3 2010/11/03 15:19:51 ngraf Exp $
*/
-
package org.lcsim.cal.calib;
import static java.lang.Math.abs;
@@ -39,6 +38,7 @@
*/
public class SamplingFractionAnalysisPolyCalDriver extends Driver
{
+
private ConditionsSet _cond;
private CollectionManager _collectionmanager = CollectionManager.defaultInstance();
// we will accumulate the raw energy values in three depths:
@@ -48,73 +48,77 @@
//
private double[][] _acc = new double[3][3];
private double[] _vec = new double[3];
-
// let's use a clusterer to remove effects of calorimeter cells hit far, far away.
// use the only cross-detector clusterer we have:
private FixedConeClusterer _fcc;
-
private AIDA aida = AIDA.defaultInstance();
private ITree _tree;
-
private boolean _initialized = false;
private boolean _debug = false;
-
-
// TODO fix this dependence on EM calorimeter geometry
boolean skipFirstLayer = false;
int firstEmStartLayer = 0;
int secondEmStartLayer = 20;
-
- double emCalInnerRadius = 0.;
+
+ private double[] _ecalLayering;
+ boolean _useFirstLayer;
+
+// double emCalInnerRadius = 0.;
double emCalInnerZ = 0.;
-
+
/** Creates a new instance of SamplingFractionAnalysisDriver */
public SamplingFractionAnalysisPolyCalDriver()
{
_tree = aida.tree();
}
-
-
+
protected void process(EventHeader event)
{
super.process(event);
//System.out.println("processing SamplingFractionAnalysisDriver");
// TODO make these values runtime definable
- String[] det = {"EcalBarrel","EcalEndcap"};
+ String[] det =
+ {
+ "EcalBarrel", "EcalEndcap"
+ };
// String[] collNames = {"EcalBarrHits", "EcalEndcapHits", "HcalBarrHits", "HcalEndcapHits"};
// double[] mipHistMaxBinVal = {.0005, .0005, .005, .005};
// double timeCut = 100.; // cut on energy depositions coming more than 100 ns late
// double ECalMipCut = .0001/3.; // determined from silicon using muons at normal incidence
// double HCalMipCut = .0008/3.; // determined from scintillator using muons at normal incidence
// double[] mipCut = {ECalMipCut, ECalMipCut, HCalMipCut, HCalMipCut};
-
-
- if(!_initialized)
+
+
+ if (!_initialized)
{
ConditionsManager mgr = ConditionsManager.defaultInstance();
try
{
_cond = mgr.getConditions("CalorimeterCalibration");
- }
- catch(ConditionsSetNotFoundException e)
+ } catch (ConditionsSetNotFoundException e)
{
- System.out.println("ConditionSet CalorimeterCalibration not found for detector "+mgr.getDetector());
+ System.out.println("ConditionSet CalorimeterCalibration not found for detector " + mgr.getDetector());
System.out.println("Please check that this properties file exists for this detector ");
}
double radius = .5;
double seed = 0.;//.1;
double minE = .05; //.25;
_fcc = new FixedConeClusterer(radius, seed, minE, FixedConeDistanceMetric.DPHIDTHETA);
-
+
// detector geometries here...
// barrel
//CylindricalCalorimeter calsubBarrel = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get(det[0]);
System.out.println("looking up subdet: " + det[0]);
- Subdetector subdet = event.getDetector().getSubdetectors().get(det[0]);
+ Subdetector subdet = event.getDetector().getSubdetectors().get(det[0]);
System.out.println("proc subdet: " + subdet.getName());
- AbstractPolyhedraCalorimeter calsubBarrel = (AbstractPolyhedraCalorimeter)event.getDetector().getSubdetectors().get(det[0]);
+
+ _ecalLayering = _cond.getDoubleArray("ECalLayering");
+ _useFirstLayer = _cond.getDouble("IsFirstEmLayerSampling")==1.;
+
+
+ AbstractPolyhedraCalorimeter calsubBarrel = (AbstractPolyhedraCalorimeter) event.getDetector().getSubdetectors().get(det[0]);
// TODO remove this hardcoded dependence on the first layer
- if(calsubBarrel.getLayering().getLayer(0).getSlices().get(0).isSensitive())
+ if (calsubBarrel.getLayering().getLayer(0).getSlices().get(0).isSensitive())
{
skipFirstLayer = true;
firstEmStartLayer += 1;
@@ -132,47 +136,49 @@
// System.out.println("Layer "+i+" slice "+j+" is "+ slice.getMaterial().getName() +" and "+(slice.isSensitive() ? " is sensitive" : ""));
// }
// }
- emCalInnerRadius = calsubBarrel.getInnerRadius();
+// emCalInnerRadius = calsubBarrel.getInnerRadius();
//endcap
- AbstractPolyhedraCalorimeter calsubEndcap = (AbstractPolyhedraCalorimeter)event.getDetector().getSubdetectors().get(det[1]);
- emCalInnerZ = abs(calsubEndcap.getZMin());
- if(skipFirstLayer) System.out.println("processing "+event.getDetectorName()+" with an em calorimeter with a massless first gap");
- System.out.println("Calorimeter bounds: r= "+emCalInnerRadius+ " z= "+emCalInnerZ);
+// AbstractPolyhedraCalorimeter calsubEndcap = (AbstractPolyhedraCalorimeter) event.getDetector().getSubdetectors().get(det[1]);
+// emCalInnerZ = abs(calsubEndcap.getZMin());
+ if (skipFirstLayer)
+ System.out.println("processing " + event.getDetectorName() + " with an em calorimeter with a massless first gap");
+// System.out.println("Calorimeter bounds: r= " + emCalInnerRadius + " z= " + emCalInnerZ);
System.out.println("initialized...");
-
+
_initialized = true;
}
-
+
// organize the histogram tree by species and energy
List<MCParticle> mcparts = event.getMCParticles();
- MCParticle mcpart = mcparts.get(mcparts.size()-1);
+ 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)
+ if (mcEnergy < .99)
{
- mcIntegerEnergy = Math.round(mcEnergy*1000);
+ 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"));
-
+ _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)
+// 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)
+ if (mcpart.getSimulatorStatus().isDecayedInCalorimeter())
{
// // now let's check the em calorimeters...
// // get all of the calorimeter hits...
@@ -207,7 +213,8 @@
String processedHitsName = _cond.getString("ProcessedHitsCollectionName");
List<CalorimeterHit> hitsToCluster = _collectionmanager.getList(processedHitsName);//event.get(CalorimeterHit.class, processedHitsName);
- if(_debug) System.out.println("clustering "+hitsToCluster.size()+" hits");
+ if (_debug)
+ System.out.println("clustering " + hitsToCluster.size() + " hits");
// quick check
// for(CalorimeterHit hit : hitsToCluster)
// {
@@ -216,30 +223,31 @@
// }
// cluster the hits
List<Cluster> clusters = _fcc.createClusters(hitsToCluster);
- if(_debug) System.out.println("found "+clusters.size()+" clusters");
+ if (_debug)
+ System.out.println("found " + clusters.size() + " clusters");
aida.histogram1D("number of found clusters", 10, -0.5, 9.5).fill(clusters.size());
- for(Cluster c : clusters)
+ for (Cluster c : clusters)
{
// System.out.println(c);
aida.cloud1D("cluster energy for all clusters").fill(c.getEnergy());
}
-
+
// proceed only if we found a single cluster above threshold
// too restrictive! simply take the highest energy cluster
- if(clusters.size()>0)
+ if (clusters.size() > 0)
{
Cluster c = clusters.get(0);
-
+
aida.cloud1D("Highest energy cluster energy").fill(c.getEnergy());
aida.cloud1D("Highest energy cluster number of cells").fill(c.getCalorimeterHits().size());
-
+
double clusterEnergy = c.getEnergy();
double mcMass = mcpart.getType().getMass();
// subtract the mass to get kinetic energy...
- double expectedEnergy = sqrt(mcEnergy*mcEnergy - mcMass*mcMass);
+ double expectedEnergy = sqrt(mcEnergy * mcEnergy - mcMass * mcMass);
// System.out.println(mcpart.getType().getName()+" "+expectedEnergy);
aida.cloud1D("measured - predicted energy").fill(clusterEnergy - expectedEnergy);
-
+
// let's now break down the cluster by component.
// this analysis uses:
// 1.) first 20 EM layers
@@ -248,7 +256,7 @@
List<CalorimeterHit> hits = c.getCalorimeterHits();
double[] vals = new double[3];
double clusterRawEnergy = 0.;
- for(CalorimeterHit hit : hits)
+ for (CalorimeterHit hit : hits)
{
long id = hit.getCellID();
IDDecoder decoder = hit.getIDDecoder();
@@ -257,19 +265,18 @@
String detectorName = decoder.getSubdetector().getName();
int type = 0;
// FIXME Hard-coded name.
- if(detectorName.startsWith("Ecal"))
+ if (detectorName.startsWith("Ecal"))
{
- if(layer>=firstEmStartLayer && layer<secondEmStartLayer)
+ if (layer >= firstEmStartLayer && layer < secondEmStartLayer)
{
type = 0;
- }
- else
+ } else
{
type = 1;
}
}
// FIXME Hard-coded name.
- if(detectorName.startsWith("Hcal"))
+ if (detectorName.startsWith("Hcal"))
{
type = 2;
}
@@ -278,62 +285,60 @@
} // end of loop over hits in cluster
// set up linear least squares:
// expectedEnergy = a*E1 + b*E2 +c*E3
- for(int j=0; j<3; ++j)
+ for (int j = 0; j < 3; ++j)
{
- _vec[j] += expectedEnergy*vals[j];
- for(int k=0; k<3; ++k)
+ _vec[j] += expectedEnergy * vals[j];
+ for (int k = 0; k < 3; ++k)
{
- _acc[j][k] += vals[j]*vals[k];
+ _acc[j][k] += vals[j] * vals[k];
}
}
} // end of single cluster cut
-
+
// event.put("All Calorimeter Hits",allHits);
// event.put("Hits To Cluster",hitsToCluster);
- event.put("Found Clusters",clusters);
+ event.put("Found Clusters", clusters);
}// end of check on decays outside tracker volume
_tree.cd("/");
}
-
+
protected void endOfData()
{
System.out.println("done! endOfData.");
// calculate the sampling fractions...
Matrix A = new Matrix(_acc, 3, 3);
- A.print(6,4);
- Matrix b = new Matrix(3,1);
- for(int i=0; i<3; ++i)
+ A.print(6, 4);
+ Matrix b = new Matrix(3, 1);
+ for (int i = 0; i < 3; ++i)
{
- b.set(i,0,_vec[i]);
+ b.set(i, 0, _vec[i]);
}
- b.print(6,4);
+ b.print(6, 4);
try
{
Matrix x = A.solve(b);
x.print(6, 4);
- System.out.println("SamplingFractions: "+ (skipFirstLayer ? "1., ":"")+1./x.get(0,0)+ ", " + 1./x.get(1,0)+", "+1./x.get(2,0));
- }
- catch(Exception e)
+ System.out.println("SamplingFractions: " + (skipFirstLayer ? "1., " : "") + 1. / x.get(0, 0) + ", " + 1. / x.get(1, 0) + ", " + 1. / x.get(2, 0));
+ } catch (Exception e)
{
e.printStackTrace();
System.out.println("try reducing dimensionality...");
Matrix Ap = new Matrix(_acc, 2, 2);
- Ap.print(6,4);
- Matrix bp = new Matrix(2,1);
- for(int i=0; i<2; ++i)
+ Ap.print(6, 4);
+ Matrix bp = new Matrix(2, 1);
+ for (int i = 0; i < 2; ++i)
{
- bp.set(i,0,_vec[i]);
+ bp.set(i, 0, _vec[i]);
}
- bp.print(6,4);
+ bp.print(6, 4);
try
{
Matrix x = Ap.solve(bp);
x.print(6, 4);
- }
- catch(Exception e2)
+ } catch (Exception e2)
{
e2.printStackTrace();
}
}
}
-}
\ No newline at end of file
+}
