lcsim/src/org/lcsim/contrib/SteveMagill
diff -u -r1.1 -r1.2
--- PPFAPlotDriver.java 23 Apr 2007 20:05:34 -0000 1.1
+++ PPFAPlotDriver.java 18 May 2007 19:58:28 -0000 1.2
@@ -32,6 +32,7 @@
private String _perfname;
private String _mcname;
+ private String _pfaname;
private String _perfjetname;
private String _pfajetname;
private AIDA aida = AIDA.defaultInstance();
@@ -39,6 +40,7 @@
public void process(EventHeader event)
{
+ List<ReconstructedParticle> pfarps = event.get(ReconstructedParticle.class,_pfaname);
List<ReconstructedParticle> perfrps = event.get(ReconstructedParticle.class,_perfname);
List<MCReconstructedParticle> mcrps = event.get(MCReconstructedParticle.class,_mcname);
List<ReconstructedParticle> perfjets = event.get(ReconstructedParticle.class,_perfjetname);
@@ -48,92 +50,52 @@
List<BasicCluster> phoeccls = event.get(BasicCluster.class,"PhotonECClusters");
List<ReconstructedParticle> ppjets = event.get(ReconstructedParticle.class,"PerfectJets");
-// List<CheatCluster> chHADclusters = event.get(CheatCluster.class,"HcalBarrDigiHitsCheatClusters");
-
- // here's some jet stuff for perfect jets
- double jESum = 0;
- double jPXSum = 0;
- double jPYSum = 0;
- double jPZSum = 0;
- double djm = 0;
- double[] jpz = new double[2];
- double[] jp = new double[2];
- int j = 0;
- if (perfPFAD) aida.cloud1D("Number of Perfect Jets").fill(perfjets.size());
- if (perfjets.size() == 2)
- {
- // calculate dijet mass - don't we have something that does this?
- for (ReconstructedParticle perfjet : perfjets)
- {
- if (perfPFAD) aida.cloud1D("Perfect Jet Energy").fill(perfjet.getEnergy());
- jESum += perfjet.getEnergy();
- Hep3Vector j3v = perfjet.getMomentum();
- jPXSum += j3v.x();
- jPYSum += j3v.y();
- jPZSum += j3v.z();
- jp[j] = Math.sqrt(j3v.x()*j3v.x()+j3v.y()*j3v.y()+j3v.z()*j3v.z());
- jpz[j] = j3v.z();
- j++;
- }
- double jE2 = jESum*jESum;
- double jP2 = jPXSum*jPXSum+jPYSum*jPYSum+jPZSum*jPZSum;
- djm += Math.sqrt(jE2-jP2);
- if (perfPFAD && jpz[0]/jp[0]<0.8 && jpz[1]/jp[1]<0.8) aida.cloud1D("Perf Barrel DiJet Mass").fill(djm);
- if (perfPFAD) aida.cloud1D("Perfect Dijet Mass").fill(djm);
- if (perfPFAD) aida.cloud1D("Diff Perfect Dijet Mass").fill((djm-91)/9.5);
- }
-
- // now real PFA jets
- double rjESum = 0;
- double rjPXSum = 0;
- double rjPYSum = 0;
- double rjPZSum = 0;
- double rdjm = 0;
- double[] jrpz = new double[2];
- double[] jrp = new double[2];
- int k = 0;
- if (perfPFAD) aida.cloud1D("Number of PFA Jets").fill(pfajets.size());
- if (pfajets.size() == 2)
- {
- // calculate dijet mass - don't we have something that does this?
- for (ReconstructedParticle pfajet : pfajets)
- {
- if (perfPFAD) aida.cloud1D("PFA Jet Energy").fill(pfajet.getEnergy());
- rjESum += pfajet.getEnergy();
- Hep3Vector rj3v = pfajet.getMomentum();
- rjPXSum += rj3v.x();
- rjPYSum += rj3v.y();
- rjPZSum += rj3v.z();
- jrp[k] = Math.sqrt(rj3v.x()*rj3v.x()+rj3v.y()*rj3v.y()+rj3v.z()*rj3v.z());
- jrpz[k] = rj3v.z();
- k++;
- }
- double rjE2 = rjESum*rjESum;
- double rjP2 = rjPXSum*rjPXSum+rjPYSum*rjPYSum+rjPZSum*rjPZSum;
- rdjm += Math.sqrt(rjE2-rjP2);
- if (perfPFAD) aida.cloud1D("PFA Dijet Mass").fill(rdjm);
- if (perfPFAD) aida.cloud1D("Diff PFA Dijet Mass").fill((rdjm-91)/9.5);
- if (perfPFAD && jrpz[0]/jrp[0]<0.8 && jrpz[1]/jrp[1]<0.8) aida.cloud1D("PFA Barrel DiJet Mass").fill(rdjm);
- if (perfPFAD && rdjm>0 && djm>0) aida.cloud1D("Difference PFAJet and PerfJet DiJet Mass").fill((rdjm-djm)/Math.sqrt(djm));
- }
-
- // Check some particle comparisons - photons first
- int NPFAPho = phobcls.size()+phoeccls.size();
- double PFAPhoClE = 0;
- for (BasicCluster phobcl : phobcls)
- {
- PFAPhoClE += phobcl.getEnergy()*1.015*1.013;
- }
- for (BasicCluster phoeccl : phoeccls)
- {
- PFAPhoClE += phoeccl.getEnergy()*1.015*1.013;
+ // Compare reconstructed particle lists
+ int npfapho = 0;
+ double rpphoE = 0;
+ int npfaneu = 0;
+ double rpneuE = 0;
+ int npfachh = 0;
+ double rpchhE = 0;
+ for (ReconstructedParticle pfarp : pfarps)
+ {
+ double rpch = pfarp.getCharge();
+ double rpmass = pfarp.getMass();
+ double rpE = pfarp.getEnergy();
+ if (rpch == 0)
+ {
+ if (rpmass == 0)
+ {
+ rpphoE += rpE;
+ npfapho++;
+ }else
+ {
+ rpneuE += rpE;
+ npfaneu++;
+ }
+ } else
+ {
+ rpchhE += rpE;
+ npfachh++;
+ }
}
+ aida.cloud1D("Number of PFA Photons").fill(npfapho);
+ aida.cloud1D("Number of PFA Neutral Hadrons").fill(npfaneu);
+ aida.cloud1D("Number of PFA Charged Hadrons").fill(npfachh);
+ aida.cloud1D("ESum of PFA Photons").fill(rpphoE);
+ aida.cloud1D("ESum of PFA Neutral Hadrons").fill(rpneuE);
+ aida.cloud1D("ESum of PFA Charged Hadrons").fill(rpchhE);
+ aida.cloud1D("Total Number of PFA Particles").fill(npfapho+npfaneu+npfachh);
+ aida.cloud1D("Total PFA Particle ESum").fill(rpphoE+rpneuE+rpchhE);
// now some stuff to compare perfect reco particles to true particles
- double perfpho = 0;
- double perfchpE = 0;
- double perfchpP = 0;
- double perfneu = 0;
+ double perfphoE = 0;
+ int nperfpho = 0;
+ int nperfchh = 0;
+ double perfchhE = 0;
+ double perfchhP = 0;
+ double perfneuE = 0;
+ int nperfneu = 0;
for (ReconstructedParticle perfrp : perfrps)
{
double perfpCh = perfrp.getCharge();
@@ -143,36 +105,94 @@
double perfpE = perfrp.getEnergy();
Hep3Vector perfpP = perfrp.getMomentum();
double pP = perfpP.magnitude();
+ double pmass = perfrp.getMass();
if (perfpCh == 0)
{
- if (perfpID == 22)
+ if (pmass == 0)
{
- perfpho += perfpE;
+ perfphoE += perfpE;
+ nperfpho++;
}else
{
- perfneu += perfpE;
+ perfneuE += perfpE;
+ nperfneu++;
}
} else
{
- perfchpE += perfpE;
- perfchpP += pP;
+ perfchhE += perfpE;
+ nperfchh++;
+ perfchhP += pP;
}
for (MCReconstructedParticle mcrp : mcrps)
{
MCParticle mcp = mcrp.getMCParticle();
}
}
+ aida.cloud1D("Number of Perfect Photons").fill(nperfpho);
+ aida.cloud1D("Number of Perfect Neutral Hadrons").fill(nperfneu);
+ aida.cloud1D("Number of Perfect Charged Hadrons").fill(nperfchh);
+ aida.cloud1D("ESum of Perfect Photons").fill(perfphoE);
+ aida.cloud1D("ESum of Perfect Neutral Hadrons").fill(perfneuE);
+ aida.cloud1D("ESum of Perfect Charged Hadrons").fill(perfchhE);
+ aida.cloud1D("Total Number of Perfect Particles").fill(nperfpho+nperfneu+nperfchh);
+ aida.cloud1D("Total Perfect Particle ESum").fill(perfphoE+perfneuE+perfchhE);
+
+ // compare some values perfect and real
+ if (perfPFAD)
+ {
+ aida.cloud1D("Difference NPFAPhotons NPerfPhotons").fill(npfapho-nperfpho);
+ aida.cloud1D("Difference PFAPhoE PerfPhoE").fill(rpphoE-perfphoE);
+ aida.cloud1D("Difference NPFANeutralH NPerfNeutralH").fill(npfaneu-nperfneu);
+ aida.cloud1D("Difference PFANeuE PerfNeuE").fill(rpneuE-perfneuE);
+ aida.cloud1D("Difference NPFA Particles NPerfect Particles").fill(npfapho+npfaneu+npfachh-nperfpho-nperfneu-nperfchh);
+ aida.cloud1D("Difference PFA ESum Perfect ESum").fill(rpphoE+rpneuE+rpchhE-perfphoE-perfneuE-perfchhE);
+ aida.cloud1D("Perfect E over P for Charge Hadrons").fill(perfchhP/perfchhE);
+ aida.cloud1D("Perf Pho ESum Ratio").fill(perfphoE/(perfphoE+perfneuE+perfchhE));
+ aida.cloud1D("Perf Ch H ESum Ratio").fill(perfchhE/(perfphoE+perfneuE+perfchhE));
+ aida.cloud1D("Perf Neu H ESum Ratio").fill(perfneuE/(perfphoE+perfneuE+perfchhE));
+ }
+
+ // Check some particle comparisons - photons first
+ int NPFAPho = phobcls.size()+phoeccls.size();
+ double PFAPhoClE = 0;
+ for (BasicCluster phobcl : phobcls)
+ {
+ PFAPhoClE += phobcl.getEnergy()*1.015*1.013;
+ }
+ for (BasicCluster phoeccl : phoeccls)
+ {
+ PFAPhoClE += phoeccl.getEnergy()*1.015*1.013;
+ }
+ aida.cloud1D("PFA Photon ESum").fill(PFAPhoClE);
+
+ // neutral hadrons
+ try
+ {
+ List<BasicCluster> neucls = event.get(BasicCluster.class,"NeuHClusters");
+ int NPFANeu = neucls.size();
+ double PFAneuClE = 0;
+ for (BasicCluster neucl : neucls)
+ {
+ PFAneuClE += neucl.getEnergy()*1.3;
+ }
+ aida.cloud1D("PFA Neutral ESum").fill(PFAneuClE);
+ aida.cloud1D("PFA Photon and NeutralH ESum").fill(PFAPhoClE+PFAneuClE);
+ }
+ catch(java.lang.IllegalArgumentException ex)
+ {
+ System.out.println("No Neutral Clusters");
+ }
// Check the cheated clusters also
double chclusE = 0;
- int NPerfPho = 0;
+ int NPerfPhoCl = 0;
double PerfPhoClE = 0;
for (CheatCluster chclus : chclusters)
{
chclusE += chclus.getEnergy();
if (chclus.getMCParticle().getCharge() == 0 && chclus.getMCParticle().getPDGID() == 22)
{
- NPerfPho++;
+ NPerfPhoCl++;
PerfPhoClE += chclus.getEnergy();
}
if (chclus.getMCParticle().getCharge() != 0)
@@ -183,20 +203,6 @@
}
}
if (perfPFAD) aida.cloud1D("Cheated CLuster E").fill(chclusE);
- if (perfPFAD) aida.cloud1D("Difference NPFAPhotons NPerfPhotons").fill(NPFAPho-NPerfPho);
- if (perfPFAD) aida.cloud1D("Difference PFAPhoE PerfPhoE").fill(PFAPhoClE-PerfPhoClE);
-
- double TPESum = perfpho+perfchpP+perfneu;
- if (perfPFAD) aida.cloud1D("Total Perfect ESum").fill(TPESum);
- if (perfPFAD) aida.cloud1D("Perfect Photon ESum").fill(perfpho);
- if (perfPFAD) aida.cloud1D("Perfect Charge P ESum").fill(perfchpE);
- if (perfPFAD) aida.cloud1D("Perfect Charge P PSum").fill(perfchpP);
- if (perfPFAD) aida.cloud1D("Perfect E over P for Charge P").fill(perfchpP/perfchpE);
- if (perfPFAD) aida.cloud1D("Perfect Neutral H ESum").fill(perfneu);
- if (perfPFAD) aida.cloud1D("Perfect Charge P + Photon ESum").fill(perfchpP+perfpho);
- if (perfPFAD) aida.cloud1D("Perf Pho ESum Ratio").fill(perfpho/TPESum);
- if (perfPFAD) aida.cloud1D("Perf Ch P ESum Ratio").fill(perfchpP/TPESum);
- if (perfPFAD) aida.cloud1D("Perf Neu H ESum Ratio").fill(perfneu/TPESum);
}
public void setPerfPartName(String perfname)
@@ -207,6 +213,10 @@
{
_mcname = mcname;
}
+ public void setPFAPartName(String pfaname)
+ {
+ _pfaname = pfaname;
+ }
public void setPerfJetName(String perfjetname)
{
_perfjetname = perfjetname;
