lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.2 -r1.3
--- MassAnalysis.java 10 Apr 2007 20:43:08 -0000 1.2
+++ MassAnalysis.java 23 Apr 2007 21:33:38 -0000 1.3
@@ -5,6 +5,8 @@
import hep.physics.vec.*;
import hep.physics.jet.*;
import hep.aida.*;
+
+import org.lcsim.event.EventHeader.LCMetaData;
import org.lcsim.util.*;
import org.lcsim.util.decision.*;
@@ -56,14 +58,20 @@
ICloud1D m_h1resAdHoc;
ICloud1D m_h2resAdHoc;
ICloud1D m_h3resAdHoc;
+ ICloud1D m_h1resNoMerge;
+ ICloud1D m_h2resNoMerge;
+ ICloud1D m_h3resNoMerge;
ICloud1D m_h1trueRes;
ICloud1D m_h1trueResAdHoc;
+ ICloud1D m_h1trueResNoMerge;
ICloud1D m_h1PassingCut;
ICloud1D m_h1trueResPassingCut;
ICloud1D m_h1trueResPassingCutAdHoc;
+ ICloud1D m_h1trueResPassingCutNoMerge;
ICloud1D m_h1PassingCut2;
ICloud1D m_h1trueResPassingCut2;
ICloud1D m_h1trueResPassingCut2AdHoc;
+ ICloud1D m_h1trueResPassingCut2NoMerge;
public MassAnalysis() {
// PFA
@@ -83,6 +91,11 @@
jNonTrivialAdHoc.setOutputCollectionName("jetOutputAdHoc");
jNonTrivialAdHoc.setFinder(twojet);
add(jNonTrivialAdHoc);
+ JetDriver jNonTrivialNoMerge = new JetDriver();
+ jNonTrivialNoMerge.setInputCollectionName("all particles [no frag merge] (ron calib)");
+ jNonTrivialNoMerge.setOutputCollectionName("jetOutputNoMerge");
+ jNonTrivialNoMerge.setFinder(twojet);
+ add(jNonTrivialNoMerge);
JetDriver jTrivial = new JetDriver();
jTrivial.setInputCollectionName("perfect particles");
jTrivial.setOutputCollectionName("jetOutputTrivial");
@@ -105,17 +118,23 @@
m_h1resAdHoc = m_histoFactory.createCloud1D("eventMassResidualsAdHoc");
m_h2resAdHoc = m_histoFactory.createCloud1D("eventEnergyResidualsAdHoc");
m_h3resAdHoc = m_histoFactory.createCloud1D("eventMomentumResidualsAdHoc");
+ m_h1resNoMerge = m_histoFactory.createCloud1D("eventMassResidualsNoMerge");
+ m_h2resNoMerge = m_histoFactory.createCloud1D("eventEnergyResidualsNoMerge");
+ m_h3resNoMerge = m_histoFactory.createCloud1D("eventMomentumResidualsNoMerge");
m_h1perf = m_histoFactory.createCloud1D("eventMassPerfect");
m_h2perf = m_histoFactory.createCloud1D("eventEnergyPerfect");
m_h3perf = m_histoFactory.createCloud1D("eventMomentumPerfect");
m_h1trueRes = m_histoFactory.createCloud1D("eventMassResidualsToTruth");
m_h1trueResAdHoc = m_histoFactory.createCloud1D("eventMassResidualsToTruthAdHoc");
+ m_h1trueResNoMerge = m_histoFactory.createCloud1D("eventMassResidualsToTruthNoMerge");
m_h1PassingCut = m_histoFactory.createCloud1D("eventMassInBarrel");
m_h1trueResPassingCut = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel");
m_h1trueResPassingCutAdHoc = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrelAdHoc");
+ m_h1trueResPassingCutNoMerge = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrelNoMerge");
m_h1PassingCut2 = m_histoFactory.createCloud1D("eventMassInBarrel2");
m_h1trueResPassingCut2 = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel2");
m_h1trueResPassingCut2AdHoc = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel2AdHoc");
+ m_h1trueResPassingCut2NoMerge = m_histoFactory.createCloud1D("eventMassResidualsToTruthInBarrel2NoMerge");
} catch (IOException ioe1) {
ioe1.printStackTrace();
}
@@ -126,6 +145,15 @@
{
super.process(event);
+ {
+ org.lcsim.util.event.BaseLCSimEvent tmpEvent = (org.lcsim.util.event.BaseLCSimEvent)(event);
+ java.util.Collection<LCMetaData> metaData = tmpEvent.MatGetMetaData();
+ for (LCMetaData meta : metaData) {
+ System.out.println("DEBUG: Event contains list '"+meta.getName()+"' of type "+meta.getType());
+ }
+ }
+
+
// Event cut: Both reconstructed jets have cos(theta) < 0.8
// Event cut 2: Both quarks have cos(theta) < 0.8
boolean passesTruthAcceptanceCut = false;
@@ -173,86 +201,41 @@
}
}
- // DEBUG: Cross-check event energy sum
- if (false) {
- List<ReconstructedParticle>recoParticles = event.get(ReconstructedParticle.class,"all particles (ron calib)");
- List<ReconstructedParticle>recoParticlesAdHoc = event.get(ReconstructedParticle.class,"all particles");
- List<ReconstructedParticle>perfParticles = event.get(ReconstructedParticle.class,"perfect particles");
- double recoSum = 0.0;
- double recoSumAdHoc = 0.0;
- double perfSum = 0.0;
- Hep3Vector perfSumVec = new BasicHep3Vector();
- Hep3Vector recoSumVec = new BasicHep3Vector();
- for (ReconstructedParticle part : recoParticles) {
- recoSum += part.getEnergy();
- recoSumVec = VecOp.add(recoSumVec, part.getMomentum());
- // Debug
- double partEnergy = part.getEnergy();
- Hep3Vector partMomVec = part.getMomentum();
- double partMomMag = partMomVec.magnitude();
- double partMass = Math.sqrt(partEnergy*partEnergy - partMomMag*partMomMag);
- System.out.println("DEBUG/RECO: Particle has Q="+part.getCharge()+", E="+partEnergy+", p="+partMomMag+" => m="+partMass+". SumE="+recoSum+", sumP="+recoSumVec.magnitude());
-
- }
- for (ReconstructedParticle part : recoParticlesAdHoc) { recoSumAdHoc += part.getEnergy(); }
- for (ReconstructedParticle part : perfParticles) {
- perfSum += part.getEnergy();
- perfSumVec = VecOp.add(perfSumVec, part.getMomentum());
- // Debug
- double partEnergy = part.getEnergy();
- Hep3Vector partMomVec = part.getMomentum();
- double partMomMag = partMomVec.magnitude();
- double partMass = Math.sqrt(partEnergy*partEnergy - partMomMag*partMomMag);
- System.out.println("DEBUG/PERFECT: Particle has Q="+part.getCharge()+", E="+partEnergy+", p="+partMomMag+" => m="+partMass+". SumE="+perfSum+", sumP="+perfSumVec.magnitude());
- }
- }
- // DEBUG: Check just photons:
- if (false) {
- List<ReconstructedParticle>recoParticles = event.get(ReconstructedParticle.class,"large photon particles");
- for (ReconstructedParticle part : recoParticles) {
- System.out.println("DEBUG/RECO: large photon particles: Particle has Q="+part.getCharge()+", E="+part.getEnergy()+", p="+part.getMomentum().magnitude());
- }
- recoParticles = event.get(ReconstructedParticle.class,"small photons (ron calib)");
- for (ReconstructedParticle part : recoParticles) {
- System.out.println("DEBUG/RECO: small photons: Particle has Q="+part.getCharge()+", E="+part.getEnergy()+", p="+part.getMomentum().magnitude());
- }
- recoParticles = event.get(ReconstructedParticle.class,"neutral hadron particles (ron calib)");
- for (ReconstructedParticle part : recoParticles) {
- System.out.println("DEBUG/RECO: neutral hadron particles (ron calib): Particle has Q="+part.getCharge()+", E="+part.getEnergy()+", p="+part.getMomentum().magnitude());
- }
- recoParticles = event.get(ReconstructedParticle.class,"charged hadron particles 2");
- for (ReconstructedParticle part : recoParticles) {
- System.out.println("DEBUG/RECO: charged particles: Particle has Q="+part.getCharge()+", E="+part.getEnergy()+", p="+part.getMomentum().magnitude());
- }
- }
-
-
List<ReconstructedParticle>jets = event.get(ReconstructedParticle.class,"jetOutput");
List<ReconstructedParticle>jetsAdHoc = event.get(ReconstructedParticle.class,"jetOutputAdHoc");
+ List<ReconstructedParticle>jetsNoMerge = event.get(ReconstructedParticle.class,"jetOutputNoMerge");
List<ReconstructedParticle>jetsPerfect = event.get(ReconstructedParticle.class,"jetOutputTrivial");
double esum = 0.;
double esumAdHoc = 0.;
+ double esumNoMerge = 0.;
double esumPerfect = 0.;
double pxsum = 0.;
double pxsumAdHoc = 0.;
+ double pxsumNoMerge = 0.;
double pxsumPerfect = 0.;
double pysum = 0.;
double pysumAdHoc = 0.;
+ double pysumNoMerge = 0.;
double pysumPerfect = 0.;
double pzsum = 0.;
double pzsumAdHoc = 0.;
+ double pzsumNoMerge = 0.;
double pzsumPerfect = 0.;
double chargesum = 0.;
double chargesumAdHoc = 0.;
+ double chargesumNoMerge = 0.;
double chargesumPerfect = 0.;
double[] je = new double[jets.size()];
double[] jeAdHoc = new double[jets.size()];
+ double[] jeNoMerge = new double[jets.size()];
double[] jePerfect = new double[jets.size()];
Hep3Vector[] jp = new Hep3Vector[jets.size()];
Hep3Vector[] jpAdHoc = new Hep3Vector[jets.size()];
+ Hep3Vector[] jpNoMerge = new Hep3Vector[jets.size()];
Hep3Vector[] jpPerfect = new Hep3Vector[jets.size()];
int ind = 0;
int indAdHoc = 0;
+ int indNoMerge = 0;
int indPerfect = 0;
// Loop over jets
@@ -278,6 +261,17 @@
chargesumAdHoc += jet.getCharge();
indAdHoc++;
}
+ for(ReconstructedParticle jet:jetsNoMerge) {
+ esumNoMerge += jet.getEnergy();
+ jeNoMerge[indNoMerge] = jet.getEnergy();
+ Hep3Vector v = jet.getMomentum();
+ jpNoMerge[indNoMerge] = v;
+ pxsumNoMerge += v.x();
+ pysumNoMerge += v.y();
+ pzsumNoMerge += v.z();
+ chargesumNoMerge += jet.getCharge();
+ indNoMerge++;
+ }
for(ReconstructedParticle jet:jetsPerfect) {
esumPerfect += jet.getEnergy();
jePerfect[indPerfect] = jet.getEnergy();
@@ -293,13 +287,15 @@
// Compute total momentum and hence event mass
double psum = Math.sqrt(pxsum*pxsum+pysum*pysum+pzsum*pzsum);
double psumAdHoc = Math.sqrt(pxsumAdHoc*pxsumAdHoc+pysumAdHoc*pysumAdHoc+pzsumAdHoc*pzsumAdHoc);
+ double psumNoMerge = Math.sqrt(pxsumNoMerge*pxsumNoMerge+pysumNoMerge*pysumNoMerge+pzsumNoMerge*pzsumNoMerge);
double psumPerfect = Math.sqrt(pxsumPerfect*pxsumPerfect+pysumPerfect*pysumPerfect+pzsumPerfect*pzsumPerfect);
double evtmass = Math.sqrt(esum*esum - psum*psum);
double evtmassAdHoc = Math.sqrt(esumAdHoc*esumAdHoc - psumAdHoc*psumAdHoc);
+ double evtmassNoMerge = Math.sqrt(esumNoMerge*esumNoMerge - psumNoMerge*psumNoMerge);
double evtmassPerfect = Math.sqrt(esumPerfect*esumPerfect - psumPerfect*psumPerfect);
- System.out.println("DEBUG: Energy sum = "+esum+" (reco-ron), "+esumAdHoc+" (reco-AH), "+esumPerfect+" (perf)");
- System.out.println("DEBUG: Momentum sum = "+psum+" (reco-ron), "+psumAdHoc+" (reco-AH), "+psumPerfect+" (perf)");
- System.out.println("DEBUG: => mass = "+evtmass+" (reco-ron), "+evtmassAdHoc+" (reco-AH), "+evtmassPerfect+" (perf)");
+ System.out.println("DEBUG: Energy sum = "+esum+" (reco-ron), "+esumAdHoc+" (reco-AH), "+esumPerfect+" (perf), "+esumNoMerge+" (reco-ron,nomerge)");
+ System.out.println("DEBUG: Momentum sum = "+psum+" (reco-ron), "+psumAdHoc+" (reco-AH), "+psumPerfect+" (perf), "+psumNoMerge+" (reco-ron,nomerge)");
+ System.out.println("DEBUG: => mass = "+evtmass+" (reco-ron), "+evtmassAdHoc+" (reco-AH), "+evtmassPerfect+" (perf), "+evtmassNoMerge+" (reco-ron,nomerge)");
if (partZToqq != null) {
System.out.println(" => Event residuals: mass = "+(evtmass-evtmassPerfect)+", event energy = "+(esum-esumPerfect)+", mom = "+(psum-psumPerfect)+", mass-truth="+(evtmass-partZToqq.getMass()));
} else {
@@ -315,12 +311,16 @@
m_h1resAdHoc.fill(evtmassAdHoc-evtmassPerfect);
m_h2resAdHoc.fill(esumAdHoc-esumPerfect);
m_h3resAdHoc.fill(psumAdHoc-psumPerfect);
+ m_h1resNoMerge.fill(evtmassNoMerge-evtmassPerfect);
+ m_h2resNoMerge.fill(esumNoMerge-esumPerfect);
+ m_h3resNoMerge.fill(psumNoMerge-psumPerfect);
m_h1perf.fill(evtmassPerfect);
m_h2perf.fill(esumPerfect);
m_h3perf.fill(psumPerfect);
if (partZToqq != null) {
m_h1trueRes.fill(evtmass-partZToqq.getMass());
m_h1trueResAdHoc.fill(evtmassAdHoc-partZToqq.getMass());
+ m_h1trueResNoMerge.fill(evtmassNoMerge-partZToqq.getMass());
}
// Look at individual jets
@@ -339,6 +339,7 @@
if (partZToqq != null) {
m_h1trueResPassingCut.fill(evtmass-partZToqq.getMass());
m_h1trueResPassingCutAdHoc.fill(evtmassAdHoc-partZToqq.getMass());
+ m_h1trueResPassingCutNoMerge.fill(evtmassNoMerge-partZToqq.getMass());
}
}
}
@@ -349,6 +350,7 @@
m_h1PassingCut2.fill(evtmass);
m_h1trueResPassingCut2.fill(evtmass-partZToqq.getMass());
m_h1trueResPassingCut2AdHoc.fill(evtmassAdHoc-partZToqq.getMass());
+ m_h1trueResPassingCut2NoMerge.fill(evtmassNoMerge-partZToqq.getMass());
}
m_eventCount++;
