lcsim/src/org/lcsim/contrib/uiowa

diff -u -r1.3 -r1.4
--- MassPlots.java	20 Dec 2007 04:08:20 -0000	1.3
+++ MassPlots.java	25 Aug 2008 22:02:05 -0000	1.4
@@ -51,10 +51,16 @@

     ICloud1D m_h1PassingCut2;
     ICloud1D m_h2PassingCut2;
     ICloud1D m_h1trueResPassingCut2;

+    ICloud1D m_h2trueResPassingCut2;

     String m_jetListName;
     String m_inputListName;

-    List<Double> m_list_PassingCut2;
-    List<Double> m_list_trueResPassingCut2;

+    List<Double> m_listMass_PassingCut2;
+    List<Double> m_listMass_trueResPassingCut2;
+    List<Double> m_listEnergy_PassingCut2;
+    List<Double> m_listEnergy_trueResPassingCut2;
+
+    protected double m_cosThetaMin = 0.0;
+    protected double m_cosThetaMax = 0.8;

 
     public MassPlots(String inputList, String outputFilename) {
 	m_inputListName = inputList;

@@ -84,20 +90,28 @@

 	    m_h1PassingCut2 = m_histoFactory.createCloud1D("eventMassInBarrel2");
 	    m_h2PassingCut2 = m_histoFactory.createCloud1D("eventEnergyInBarrel2");
 	    m_h1trueResPassingCut2 = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel2");

+	    m_h2trueResPassingCut2 = m_histoFactory.createCloud1D("eventEnergyResidualsToTruthInBarrel2");

 	} catch (IOException ioe1) {
             ioe1.printStackTrace(); 
 	}
 
 	// For rms90 calculation:

-	m_list_PassingCut2 = new Vector<Double>();
-	m_list_trueResPassingCut2 = new Vector<Double>();

+	m_listMass_PassingCut2 = new Vector<Double>();
+	m_listMass_trueResPassingCut2 = new Vector<Double>();
+	m_listEnergy_PassingCut2 = new Vector<Double>();
+	m_listEnergy_trueResPassingCut2 = new Vector<Double>();

     }
 
     boolean m_doCheckpoints = false;
     int m_eventCount;
     protected void process(EventHeader event)
     {

-	super.process(event);

+	try {
+	    super.process(event);
+        } catch(FixNumberOfJetsFinder.NumJetsNotFoundException e) {
+	    System.out.println("NumJetsNoFoundException: "+e.getMessage()+" when trying to make mass/energy plots for '"+m_inputListName+"'");
+	    return;
+        }

 
 	// Event cut: Both reconstructed jets have cos(theta) < 0.8
 	// Event cut 2: Both quarks have cos(theta) < 0.8

@@ -133,13 +147,14 @@

 		//System.out.println("ERROR: no Z -> vv cand");
 	    }
 	    if (partZToqq == null) {

-		System.out.println("ERROR: no Z -> qq cand");

+		//System.out.println("ERROR: no Z -> qq cand");

 	    } else {
 		//System.out.println("DEBUG: Mass of Z -> qq = "+partZToqq.getMass());
 		passesTruthAcceptanceCut = true;
 		for (MCParticle dau : partZToqq_daughters) {
 		    double cosTheta = Math.abs(dau.getMomentum().z() / dau.getMomentum().magnitude());

-		    if (cosTheta>0.8) { passesTruthAcceptanceCut = false; }

+		    if (cosTheta>m_cosThetaMax) { passesTruthAcceptanceCut = false; }
+		    if (cosTheta<m_cosThetaMin) { passesTruthAcceptanceCut = false; }

 		}
 	    }
 	}

@@ -207,9 +222,12 @@

 	    m_h1PassingCut2.fill(evtmass);
 	    m_h2PassingCut2.fill(esum);
 	    m_h1trueResPassingCut2.fill(evtmass-partZToqq.getMass());

+	    m_h2trueResPassingCut2.fill(esum-partZToqq.getEnergy());

 	    //System.out.println(evtmass-partZToqq.getMass());

-	    m_list_trueResPassingCut2.add(evtmass-partZToqq.getMass());
-	    m_list_PassingCut2.add(evtmass);

+	    m_listMass_trueResPassingCut2.add(evtmass-partZToqq.getMass());
+	    m_listMass_PassingCut2.add(evtmass);
+	    m_listEnergy_trueResPassingCut2.add(esum-partZToqq.getEnergy());
+	    m_listEnergy_PassingCut2.add(esum);

 	}
 	
 

@@ -228,16 +246,24 @@

             ioe1.printStackTrace(); 
         }
 

+	printRMS90(m_listMass_trueResPassingCut2, "Mass residuals: ");
+	printRMS90(m_listEnergy_trueResPassingCut2, "Energy residuals: ");
+
+	// Done
+        super.suspend();
+    }
+
+    protected void printRMS90(List<Double> data, String info) {

 	// Compute rms90
 	double targetCentralFraction = 0.9;

-	Collections.sort(m_list_trueResPassingCut2);

+	Collections.sort(data);

-	int nPoints = m_list_trueResPassingCut2.size();

+	int nPoints = data.size();

 	int nTail = (int)((1.0 - targetCentralFraction) * nPoints);
 	int nNonTail = nPoints - nTail;
 	double centralFraction = ((double)(nNonTail))/((double)(nPoints));
 

-	System.out.println("Computing rms90 for "+m_inputListName+": nPoints = "+nPoints+" (of which "+nTail+" tail and "+nNonTail+" non-tail for a central fraction of "+centralFraction+")");

+	System.out.println("Computing rms90 for "+m_inputListName+": "+info+"nPoints = "+nPoints+" (of which "+nTail+" tail and "+nNonTail+" non-tail for a central fraction of "+centralFraction+")");

 
 	double rms_min = 100000.0;
 	double mean_min = 0.0;

@@ -248,14 +274,14 @@

 	    double xm2 = 0.0;
 	    double xn = 0.0;
 	    for (int j=i; j<i+nNonTail; j++) {

-		xm += m_list_trueResPassingCut2.get(j);
-		xm2 += (m_list_trueResPassingCut2.get(j) * m_list_trueResPassingCut2.get(j));

+		xm += data.get(j);
+		xm2 += (data.get(j) * data.get(j));

 		xn += 1.0;
 	    }
 	    double rms = Math.sqrt(xm2/xn - (xm/xn)*(xm/xn));
 	    if (rms<rms_min) {

-		output_min = m_list_trueResPassingCut2.get(i);
-		output_max = m_list_trueResPassingCut2.get(i+nNonTail-1);

+		output_min = data.get(i);
+		output_max = data.get(i+nNonTail-1);

 		rms_min = rms;
 		mean_min = xm/xn;
 	    }

@@ -265,9 +291,22 @@

 	System.out.println("mean90:   "+mean_min);
 	System.out.println("min90:    "+output_min);
 	System.out.println("max90:    "+output_max);

-	System.out.println("Full RMS: "+m_h1trueResPassingCut2.rms());

-	// Done
-        super.suspend();

+	// Check: Compute full RMS
+	double full_xm = 0.0;
+	double full_xm2 = 0.0;
+	double full_xn = 0.0;
+	for (Double point : data) {
+	    double value = point.doubleValue();
+	    full_xm += value;
+	    full_xm2 += (value * value);
+	    full_xn += 1.0;
+	}
+	double fullRMS = Math.sqrt(full_xm2/full_xn - (full_xm/full_xn)*(full_xm/full_xn));
+	System.out.println("Full RMS: "+fullRMS);
+

     }

+
+    public void setCosThetaMin(double cut) { m_cosThetaMin = cut; }
+    public void setCosThetaMax(double cut) { m_cosThetaMax = cut; }

 }