lcsim/src/org/lcsim/recon/pfa/output
diff -u -r1.3 -r1.4
--- MassPlots.java 5 Jun 2008 00:18:02 -0000 1.3
+++ MassPlots.java 10 Feb 2009 19:03:00 -0000 1.4
@@ -52,57 +52,95 @@
ICloud1D m_h2trueResPassingCut2;
String m_jetListName;
String m_inputListName;
- List<Double> m_list_PassingCut2;
- List<Double> m_list_trueResPassingCut2;
+ String m_outputFileName;
+ List<Double> m_listMass_PassingCut2;
+ List<Double> m_listMass_trueResPassingCut2;
+ List<Double> m_listEnergy_PassingCut2;
+ List<Double> m_listEnergy_trueResPassingCut2;
+
+ protected double m_cosThetaMin;
+ protected double m_cosThetaMax;
+
+ public void setCosThetaMin(double cut) { m_cosThetaMin = cut; }
+ public void setCosThetaMax(double cut) { m_cosThetaMax = cut; }
+ public void setInputListName(String s) {
+ m_inputListName = s;
+ m_jetListName = "jetOutput__" + s;
+ jNonTrivial.setInputCollectionName(m_inputListName);
+ jNonTrivial.setOutputCollectionName(m_jetListName);
+ jNonTrivial.setFinder(twojet);
+ }
+ public void setOutputFileName(String s) { m_outputFileName = s; }
+
- public MassPlots(String inputList, String outputFilename) {
- m_inputListName = inputList;
- m_jetListName = new String();
- m_jetListName += "jetOutput__";
- m_jetListName += inputList;
- JetFinder twojet = new FixNumberOfJetsFinder(2);
- JetDriver jNonTrivial = new JetDriver();
- jNonTrivial.setInputCollectionName(inputList);
- jNonTrivial.setOutputCollectionName(m_jetListName);
- jNonTrivial.setFinder(twojet);
+ JetFinder twojet = new FixNumberOfJetsFinder(2);
+ JetDriver jNonTrivial = new JetDriver();
+
+ public MassPlots() {
+ this("ReconstructedParticles", "outplot.aida", 0.0, 0.8);
+ }
+
+ public MassPlots(String inputList, String outputFilename, double cosThetaMin, double cosThetaMax) {
+
+ m_cosThetaMin = cosThetaMin;
+ m_cosThetaMax = cosThetaMax;
+
+ m_inputListName = inputList;
+ m_outputFileName = outputFilename;
+ m_jetListName = "jetOutput__" + m_inputListName;
+ jNonTrivial.setInputCollectionName(m_inputListName);
+ jNonTrivial.setOutputCollectionName(m_jetListName);
+ jNonTrivial.setFinder(twojet);
+
add(jNonTrivial);
m_eventCount = 0;
- // Init plots
- IAnalysisFactory af = IAnalysisFactory.create();
- try {
- m_tree = af.createTreeFactory().create(outputFilename, "xml", false, true);
- m_histoFactory = af.createHistogramFactory(m_tree);
- m_h1 = m_histoFactory.createCloud1D("eventMass");
- m_h2 = m_histoFactory.createCloud1D("eventEnergy");
- m_h3 = m_histoFactory.createCloud1D("eventMomentum");
- m_h1trueRes = m_histoFactory.createCloud1D("eventMassResidualsToTruth");
- 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");
- 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 startOfData()
+ {
+ // Init plots
+ IAnalysisFactory af = IAnalysisFactory.create();
+ try {
+ m_tree = af.createTreeFactory().create(m_outputFileName, "xml", false, true);
+ m_histoFactory = af.createHistogramFactory(m_tree);
+ m_h1 = m_histoFactory.createCloud1D("eventMass");
+ m_h2 = m_histoFactory.createCloud1D("eventEnergy");
+ m_h3 = m_histoFactory.createCloud1D("eventMomentum");
+ m_h1trueRes = m_histoFactory.createCloud1D("eventMassResidualsToTruth");
+ m_h1PassingCut = m_histoFactory.createCloud1D("eventMass");
+ m_h1trueResPassingCut = m_histoFactory.createCloud1D("eventMassResidualsToTruth");
+ m_h1PassingCut2 = m_histoFactory.createCloud1D("eventMass2");
+ m_h2PassingCut2 = m_histoFactory.createCloud1D("eventEnergy2");
+ m_h1trueResPassingCut2 = m_histoFactory.createCloud1D("eventMassResidualsToTruth2");
+ m_h2trueResPassingCut2 = m_histoFactory.createCloud1D("eventEnergyResidualsToTruth2");
+ } catch (IOException ioe1) {
+ ioe1.printStackTrace();
+ }
+
+ }
+
+
protected void process(EventHeader event)
{
try {
super.process(event);
} catch(FixNumberOfJetsFinder.NumJetsNotFoundException e) {
- System.out.println("NumJetsNoFoundException: "+e.getMessage());
+ System.out.println("NumJetsNoFoundException: "+e.getMessage()+" when trying to make mass/energy plots for '"+m_inputListName+"'");
return;
- }
+ } catch(java.lang.IllegalArgumentException e) {
+ System.out.println("IllegalArgumentException: "+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
@@ -138,13 +176,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; }
}
}
}
@@ -197,7 +236,7 @@
double sign = 1.;
if(masssq < 0.)sign = -1.;
double signedMass = sign*Math.sqrt(sign*masssq);
- if( (Math.abs(ct) < .8)&&(Math.abs(ct1) < .8) ) {
+ if( (Math.abs(ct) > m_cosThetaMin)&&(Math.abs(ct) < m_cosThetaMax)&&(Math.abs(ct1) > m_cosThetaMin)&&(Math.abs(ct1)<m_cosThetaMax) ) {
//System.out.println("EVENT PASSES CUT, since cos(theta) = "+ct+", "+ct1+". Jet-jet mass = "+signedMass);
m_h1PassingCut.fill(evtmass);
if (partZToqq != null) {
@@ -214,8 +253,10 @@
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);
}
@@ -225,7 +266,6 @@
//System.out.println("DEBUG: Checkpoint at "+m_eventCount);
try { m_tree.commit(); } catch(IOException ioe1) { ioe1.printStackTrace(); }
}
-
}
public void suspend() {
@@ -235,16 +275,25 @@
ioe1.printStackTrace();
}
+ System.out.println("--> "+ m_cosThetaMin + " < cos theta < " + m_cosThetaMax );
+ 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;
@@ -255,14 +304,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;
}
@@ -272,10 +321,19 @@
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);
+ }
}
