lcsim/src/org/lcsim/contrib/uiowa

diff -u -r1.1 -r1.2
--- MassPlots.java	27 Aug 2007 21:27:39 -0000	1.1
+++ MassPlots.java	26 Nov 2007 00:39:15 -0000	1.2
@@ -49,10 +49,15 @@

     ICloud1D m_h1PassingCut;
     ICloud1D m_h1trueResPassingCut;
     ICloud1D m_h1PassingCut2;

+    ICloud1D m_h2PassingCut2;

     ICloud1D m_h1trueResPassingCut2;
     String m_jetListName;

+    String m_inputListName;
+    List<Double> m_list_PassingCut2;
+    List<Double> m_list_trueResPassingCut2;

 
     public MassPlots(String inputList, String outputFilename) {

+	m_inputListName = inputList;

 	m_jetListName = new String();
 	m_jetListName += "jetOutput__";
 	m_jetListName += inputList;

@@ -77,12 +82,18 @@

 	    m_h1PassingCut = m_histoFactory.createCloud1D("eventMassInBarrel");
 	    m_h1trueResPassingCut = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel");
 	    m_h1PassingCut2 = m_histoFactory.createCloud1D("eventMassInBarrel2");

+	    m_h2PassingCut2 = m_histoFactory.createCloud1D("eventEnergyInBarrel2");

 	    m_h1trueResPassingCut2 = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel2");
 	} catch (IOException ioe1) {
             ioe1.printStackTrace(); 
 	}

+
+	// For rms90 calculation:
+	m_list_PassingCut2 = new Vector<Double>();
+	m_list_trueResPassingCut2 = new Vector<Double>();

     }
 

+    boolean m_doCheckpoints = false;

     int m_eventCount;
     protected void process(EventHeader event)
     {

@@ -124,7 +135,7 @@

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

-		System.out.println("DEBUG: Mass of Z -> qq = "+partZToqq.getMass());

+		//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());

@@ -160,7 +171,7 @@

 	double psum = Math.sqrt(pxsum*pxsum+pysum*pysum+pzsum*pzsum);
 	double evtmass = Math.sqrt(esum*esum - psum*psum);
 	if (partZToqq != null) {

-	    System.out.println("  => Event residuals: mass-truth="+(evtmass-partZToqq.getMass()));

+	    //System.out.println("  => Event residuals: mass-truth="+(evtmass-partZToqq.getMass()));

 	}
 
 	m_h1.fill(evtmass);

@@ -182,7 +193,7 @@

 		if(masssq < 0.)sign = -1.;
 		double signedMass = sign*Math.sqrt(sign*masssq);
 		if( (Math.abs(ct) < .8)&&(Math.abs(ct1) < .8) ) {

-		    System.out.println("EVENT PASSES CUT, since cos(theta) = "+ct+", "+ct1+". Jet-jet mass = "+signedMass);

+		    //System.out.println("EVENT PASSES CUT, since cos(theta) = "+ct+", "+ct1+". Jet-jet mass = "+signedMass);

 		    m_h1PassingCut.fill(evtmass);
 		    if (partZToqq != null) {
 			m_h1trueResPassingCut.fill(evtmass-partZToqq.getMass());

@@ -192,26 +203,71 @@

 	}
 
 	if (passesTruthAcceptanceCut) {

-	    System.out.println("EVENT PASSES CUT2. Jet-jet mass = "+evtmass);

+	    //System.out.println("EVENT PASSES CUT2. Jet-jet mass = "+evtmass);

 	    m_h1PassingCut2.fill(evtmass);

+	    m_h2PassingCut2.fill(esum);

 	    m_h1trueResPassingCut2.fill(evtmass-partZToqq.getMass());

+	    //System.out.println(evtmass-partZToqq.getMass());
+	    m_list_trueResPassingCut2.add(evtmass-partZToqq.getMass());
+	    m_list_PassingCut2.add(evtmass);

 	}
 	
 
 	m_eventCount++;

-	if (m_eventCount % 50 == 0) {

+	if (m_doCheckpoints && m_eventCount % 50 == 0) {

 	    // Checkpoint

-	    System.out.println("DEBUG: Checkpoint at "+m_eventCount);

+	    //System.out.println("DEBUG: Checkpoint at "+m_eventCount);

 	    try { m_tree.commit(); } catch(IOException ioe1) { ioe1.printStackTrace(); }
 	}
     }
 

-        public void suspend() {

+    public void suspend() {

         try {
             m_tree.commit();
         } catch(IOException ioe1) {
             ioe1.printStackTrace(); 
         }

+
+	// Compute rms90
+	double targetCentralFraction = 0.9;
+	Collections.sort(m_list_trueResPassingCut2);
+
+	int nPoints = m_list_trueResPassingCut2.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+")");
+
+	double rms_min = 100000.0;
+	double mean_min = 0.0;
+	double output_min = -1;
+	double output_max = -1;
+	for (int i=0; i<nTail || (nTail==0 && i==0); i++) {
+	    double xm = 0.0;
+	    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));
+		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);
+		rms_min = rms;
+		mean_min = xm/xn;
+	    }
+	}
+
+	System.out.println("rms90:    "+rms_min);
+	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();
     }
 }