Print

Print


Commit in lcsim-cal-calib/src/org/lcsim/cal/calib on MAIN
PandoraPfoSamplingFractionAnalysislDriver.java+21-51.2 -> 1.3
SamplingFractionAnalysisPolyCalDriver.java+23-71.3 -> 1.4
+44-12
2 modified files
Fixed logic to identify the initial particle to also work on samples created by GEANT particle gun
Corrected expected energy calculation

lcsim-cal-calib/src/org/lcsim/cal/calib
PandoraPfoSamplingFractionAnalysislDriver.java 1.2 -> 1.3
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
SamplingFractionAnalysisPolyCalDriver.java 1.3 -> 1.4
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
CVSspam 0.2.12


Use REPLY-ALL to reply to list

To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1