lcsim-cal-calib/src/org/lcsim/cal/calib
diff -u -r1.2 -r1.3
--- PandoraPfoSamplingFractionAnalysislDriver.java 4 Nov 2010 00:20:12 -0000 1.2
+++ PandoraPfoSamplingFractionAnalysislDriver.java 10 Feb 2012 15:26:26 -0000 1.3
@@ -1,7 +1,7 @@
/*
* PandoraPfoSamplingFractionAnalysislDriver
*
- * $Id: PandoraPfoSamplingFractionAnalysislDriver.java,v 1.2 2010/11/04 00:20:12 ngraf Exp $
+ * $Id: PandoraPfoSamplingFractionAnalysislDriver.java,v 1.3 2012/02/10 15:26:26 grefe Exp $
*/
package org.lcsim.cal.calib;
@@ -99,7 +99,19 @@
// 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); // this only works if particle was generated by stdhep, does not work with GEANT gps
+ MCParticle mcpart = null;
+ // Look for the most energetic final state particle
+ for (MCParticle myMCP : mcparts) {
+ if (myMCP.getGeneratorStatus() == MCParticle.FINAL_STATE) {
+ if (mcpart == null) {
+ mcpart = myMCP;
+ } else if (mcpart.getEnergy() < myMCP.getEnergy()) {
+ mcpart = myMCP;
+ }
+ }
+ }
+
String particleType = mcpart.getType().getName();
double mcEnergy = mcpart.getEnergy();
long mcIntegerEnergy = Math.round(mcEnergy);
@@ -163,7 +175,11 @@
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 = mcEnergy;
+ // In case of protons and neutrons we expect only the kinetic energy in the shower
+ if (mcpart.getPDGID() == 2212 || mcpart.getPDGID() == 2112) {
+ expectedEnergy = mcEnergy - mcMass;
+ }
// System.out.println(mcpart.getType().getName()+" "+expectedEnergy);
aida.cloud1D("measured - predicted energy").fill(clusterEnergy - expectedEnergy);
@@ -189,7 +205,7 @@
// System.out.println(detectorName);
int caltype = 0;
// FIXME Hard-coded name.
- if (detectorName.startsWith("Ecal"))
+ if (detectorName.toUpperCase().startsWith("ECAL"))
{
for (int i = 1; i < _ecalLayering.length + 1; ++i)
{
@@ -198,7 +214,7 @@
}
}
// FIXME Hard-coded name.
- if (detectorName.startsWith("Hcal"))
+ if (detectorName.toUpperCase().startsWith("HCAL"))
{
caltype = 3;
nHCalHits +=1;
lcsim-cal-calib/src/org/lcsim/cal/calib
diff -u -r1.3 -r1.4
--- SamplingFractionAnalysisPolyCalDriver.java 3 Nov 2010 15:19:51 -0000 1.3
+++ SamplingFractionAnalysisPolyCalDriver.java 10 Feb 2012 15:26:26 -0000 1.4
@@ -3,7 +3,7 @@
*
* Created on May 19, 2008, 11:54 AM
*
- * $Id: SamplingFractionAnalysisPolyCalDriver.java,v 1.3 2010/11/03 15:19:51 ngraf Exp $
+ * $Id: SamplingFractionAnalysisPolyCalDriver.java,v 1.4 2012/02/10 15:26:26 grefe Exp $
*/
package org.lcsim.cal.calib;
@@ -150,7 +150,19 @@
// 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); // this only works if particle was generated by stdhep, does not work with GEANT gps
+ MCParticle mcpart = null;
+ // Look for the most energetic final state particle
+ for (MCParticle myMCP : mcparts) {
+ if (myMCP.getGeneratorStatus() == MCParticle.FINAL_STATE) {
+ if (mcpart == null) {
+ mcpart = myMCP;
+ } else if (mcpart.getEnergy() < myMCP.getEnergy()) {
+ mcpart = myMCP;
+ }
+ }
+ }
+
String particleType = mcpart.getType().getName();
double mcEnergy = mcpart.getEnergy();
long mcIntegerEnergy = Math.round(mcEnergy);
@@ -244,7 +256,11 @@
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 = mcEnergy;
+ // In case of protons and neutrons we expect only the kinetic energy in the shower
+ if (mcpart.getPDGID() == 2212 || mcpart.getPDGID() == 2112) {
+ expectedEnergy = mcEnergy - mcMass;
+ }
// System.out.println(mcpart.getType().getName()+" "+expectedEnergy);
aida.cloud1D("measured - predicted energy").fill(clusterEnergy - expectedEnergy);
@@ -265,7 +281,7 @@
String detectorName = decoder.getSubdetector().getName();
int type = 0;
// FIXME Hard-coded name.
- if (detectorName.startsWith("Ecal"))
+ if (detectorName.toUpperCase().startsWith("ECAL"))
{
if (layer >= firstEmStartLayer && layer < secondEmStartLayer)
{
@@ -276,7 +292,7 @@
}
}
// FIXME Hard-coded name.
- if (detectorName.startsWith("Hcal"))
+ if (detectorName.toUpperCase().startsWith("HCAL"))
{
type = 2;
}
@@ -287,10 +303,10 @@
// expectedEnergy = a*E1 + b*E2 +c*E3
for (int j = 0; j < 3; ++j)
{
- _vec[j] += expectedEnergy * vals[j];
+ _vec[j] += vals[j];
for (int k = 0; k < 3; ++k)
{
- _acc[j][k] += vals[j] * vals[k];
+ _acc[j][k] += (vals[j] * vals[k]) / expectedEnergy ;
}
}
} // end of single cluster cut