lcsim-contrib/src/main/java/org/lcsim/contrib/SteveMagill
diff -N PFADRSelect.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PFADRSelect.java 14 Feb 2012 17:56:52 -0000 1.1
@@ -0,0 +1,362 @@
+package org.lcsim.contrib.SteveMagill;
+
+import java.io.File;
+import java.util.*;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.util.JetDriver;
+import org.lcsim.util.Driver;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.jet.*;
+import org.lcsim.util.loop.LCIODriver;
+
+/*
+ * PFASelect includes fiducial cut to select region for PFA analysis
+ */
+
+public class PFADRSelect extends Driver
+{
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ public PFADRSelect()
+ {
+ // run analysis program to make reconstructed particles
+ add(new PFACSClusdMJetDriver());
+
+// find jets with Reco Particles from analysis
+ // find jets with all tracks as input
+ JetDriver trjet = new JetDriver();
+ trjet.setInputCollectionName("AllTrackRecoParticles");
+ trjet.setOutputCollectionName("TrRPJets");
+// JetFinder tworealjets = new FixNumberOfJetsFinder(2); // uses Jade by default
+// mcjet.setFinder(tworealjets);
+ DurhamJetFinder kttrjet = new DurhamJetFinder(0.02);
+ JetFinder twokttrjets = new FixNumberOfJetsFinder(2,kttrjet);
+ trjet.setFinder(twokttrjets);
+// JetFinder sixrealjets = new FixNumberOfJetsFinder(6);
+// PFAjet.setFinder(sixrealjets);
+ // for Durham jet, ycut = 0.25 at ZPole, 0.01 at 500 GeV ZZ
+// JetFinder rktjet = new DurhamJetFinder(0.1025);
+// mcjet.setFinder(rktjet);
+ add(trjet);
+
+ // jets from Corrected Scin Cluster Reco Particles
+ JetDriver cojet = new JetDriver();
+ cojet.setInputCollectionName("CorrClusRecoParticles");
+ cojet.setOutputCollectionName("CorrRPJets");
+// JetFinder tworealjets = new FixNumberOfJetsFinder(2); // uses Jade by default
+// cojet.setFinder(tworealjets);
+ DurhamJetFinder ktcojet = new DurhamJetFinder(0.02);
+ JetFinder twoktcojets = new FixNumberOfJetsFinder(2,ktcojet);
+ cojet.setFinder(twoktcojets);
+// JetFinder sixktcojets = new FixNumberOfJetsFinder(6,ktcojet);
+// cojet.setFinder(sixktcojets);
+ // for Durham jet, ycut = 0.25 at ZPole, 0.01 at 500 GeV ZZ
+// JetFinder rktjet = new DurhamJetFinder(0.1025);
+// cojet.setFinder(rktjet);
+ add(cojet);
+
+ // find jets with MC Reco Particles
+ JetDriver mcjet = new JetDriver();
+ mcjet.setInputCollectionName("MyCheatPerfectReconParticles");
+// mcjet.setInputCollectionName("MCFSRecoParticles");
+ mcjet.setOutputCollectionName("MCRPJets");
+// JetFinder tworealjets = new FixNumberOfJetsFinder(2); // uses Jade by default
+// mcjet.setFinder(tworealjets);
+ DurhamJetFinder ktjet = new DurhamJetFinder(0.02);
+ JetFinder twoktjets = new FixNumberOfJetsFinder(2,ktjet);
+ mcjet.setFinder(twoktjets);
+// JetFinder sixktjets = new FixNumberOfJetsFinder(6,ktjet);
+// mcjet.setFinder(sixktjets);
+ // for Durham jet, ycut = 0.25 at ZPole, 0.01 at 500 GeV ZZ
+// JetFinder rktjet = new DurhamJetFinder(0.1025);
+// mcjet.setFinder(rktjet);
+ add(mcjet);
+
+ // jets from final reco particle list using PFA
+ JetDriver pfajet = new JetDriver();
+ pfajet.setInputCollectionName("CSPFARecoParticles");
+ pfajet.setOutputCollectionName("CSPFARPJets");
+// JetFinder tworealjets = new FixNumberOfJetsFinder(2); // uses Jade by default
+// pfajet.setFinder(tworealjets);
+ DurhamJetFinder ktpfajet = new DurhamJetFinder(0.02);
+ JetFinder twoktpfajets = new FixNumberOfJetsFinder(2,ktpfajet);
+ pfajet.setFinder(twoktpfajets);
+// JetFinder sixktpfajets = new FixNumberOfJetsFinder(6,ktpfajet);
+// pfajet.setFinder(sixktpfajets);
+ // for Durham jet, ycut = 0.25 at ZPole, 0.01 at 500 GeV ZZ
+// JetFinder rktjet = new DurhamJetFinder(0.1025);
+// pfajet.setFinder(rktjet);
+ add(pfajet);
+
+// do deltaM corrections to jets
+ // this is using all tracks with endpoints as determined by PFA Mip finder
+ TrBMCorr tbmdriver = new TrBMCorr();
+ tbmdriver.setInputTrackList("ReconTracks");
+ tbmdriver.setTrackSpacepointMap("TrackILPosMap");
+ tbmdriver.setDeltaMObjectName("TrDelM");
+ add(tbmdriver);
+
+ // gets DeltaM from tracks for each track jet and makes map linking
+ // jet to its deltaM using all tracks except those used in PFA
+ PFATJetBMCorr ptjbmdriver = new PFATJetBMCorr();
+ ptjbmdriver.setInputJetList("TrRPJets");
+ ptjbmdriver.setTrCClusMapName("TrackCalClusMap");
+ ptjbmdriver.setTrackSpacepointMap("TrackILPosMap");
+ ptjbmdriver.setDeltaMObjectName("PTJetDelM");
+ ptjbmdriver.setTrackDelMMap("PTJetDMMap");
+ add(ptjbmdriver);
+
+ // gets DeltaM from tracks for each track jet and makes map linking
+ // jet to its deltaM using all tracks
+ TJetBMCorr tjbmdriver = new TJetBMCorr();
+ tjbmdriver.setInputJetList("TrRPJets");
+ tjbmdriver.setTrCClusMapName("TrackCalClusMap");
+ tjbmdriver.setTrackSpacepointMap("TrackILPosMap");
+ tjbmdriver.setDeltaMObjectName("TJetDelM");
+ tjbmdriver.setTrackDelMMap("TJetDMMap");
+ add(tjbmdriver);
+
+ // links corr scin jets and track jets
+ double ctdist = 0.7;
+ STrJetLinkDriver cstdriver = new STrJetLinkDriver(ctdist);
+ String cstlinkmap = ("STrJetMap");
+ cstdriver.setInputScinJetName("CorrRPJets");
+ cstdriver.setInputTrackJetName("TrRPJets");
+ cstdriver.setOutputLinkMapName(cstlinkmap);
+ add(cstdriver);
+
+ // links PFA jets and track jets
+ double pctdist = 0.7;
+ STrJetLinkDriver pcstdriver = new STrJetLinkDriver(pctdist);
+ String pcstlinkmap = ("PSTrJetMap");
+ pcstdriver.setInputScinJetName("CSPFARPJets");
+ pcstdriver.setInputTrackJetName("TrRPJets");
+ pcstdriver.setOutputLinkMapName(pcstlinkmap);
+ add(pcstdriver);
+
+// calculate dijet masses
+ // from MC Reco Particle jets
+ MCRPdijetm djm = new MCRPdijetm();
+ djm.setInputJets("MCRPJets");
+ add(djm);
+
+ // from CS Cluster jets with delM from all tracks per jet
+ CorrRPdijetm codjm = new CorrRPdijetm();
+ codjm.setInputScinJets("CorrRPJets");
+ codjm.setSTJetMap(cstlinkmap);
+ codjm.setTJetDMMap("TJetDMMap");
+ add(codjm);
+
+ // from CS + PFA jets with delM from unassoc tracks per jet
+ PFARPdijetm pfadjm = new PFARPdijetm();
+ pfadjm.setInputScinJets("CSPFARPJets");
+ pfadjm.setSTJetMap(pcstlinkmap);
+ pfadjm.setTJetDMMap("PTJetDMMap");
+ add(pfadjm);
+
+ // write out some events maybe
+// LCIODriver ldr = new LCIODriver();
+// String[] keepCollections = {"CSPFARPJets","TrRPJets","CorrRPJets"};
+// ldr.setWriteOnlyCollections(keepCollections);
+// ldr.setOutputFilePath("C:/Documents and Settings/adminsrm/My Documents/MYoutput.slcio");
+// File output = new File(System.getProperty("Desktop"),"MYoutput.slcio");
+// add(ldr);
+
+ }
+
+ protected void process(EventHeader event)
+ {
+
+// Only process events with both jets costheta < costhetaMAX
+ List<MCParticle> mcl = event.get(MCParticle.class,"MCParticle");
+ List<Double> Ddm = new ArrayList<Double>();
+ List<Double> DdE = new ArrayList<Double>();
+ List<Double> DdpT = new ArrayList<Double>();
+ boolean keepit = true;
+ double Zmass = 0.;
+ double qqbarE = 0.;
+ double qqbarET = 0.;
+ int ievt = 0;
+ double costmax = 0.9; // 0.8 is barrel, max ~.95 or so to avoid beampipe
+ for(MCParticle p:mcl)
+ {
+ int id = Math.abs(p.getPDGID());
+ if( (id == 1 )||(id == 2 )||(id == 3 ) )
+ {
+ if(p.getParents().get(0).getPDGID() == 23 )
+ {
+ Zmass = p.getParents().get(0).getMass();
+ qqbarE = p.getParents().get(0).getEnergy();
+ Ddm.add(Zmass);
+ DdE.add(qqbarE);
+// System.out.println("Z Mass " + Zmass);
+ Hep3Vector pp = p.getMomentum();
+ double cost = Math.abs(pp.z()/pp.magnitude());
+ if (cost>0. && cost<1.) qqbarET = qqbarE*cost;
+ DdpT.add(qqbarET);
+ if(cost > costmax)
+ {
+ keepit = false;
+ continue;
+ }
+ }
+ }
+ }
+
+ if(!keepit)
+ {
+ ievt++;
+ return;
+ }
+
+ event.put("DJmass", Ddm);
+ event.put("DJEnergy", DdE);
+ event.put("DJpT", DdpT);
+
+ super.process(event);
+
+ // add some plots here - get jet stuff
+ List<ReconstructedParticle> perfjets = event.get(ReconstructedParticle.class,"MCRPJets");
+ List<ReconstructedParticle> pfajets = event.get(ReconstructedParticle.class,"CSPFARPJets");
+ List<ReconstructedParticle> csscinjets = event.get(ReconstructedParticle.class,"CorrRPJets");
+ List<ReconstructedParticle> rpList = event.get(ReconstructedParticle.class,"CSPFARecoParticles");
+
+ // check RP stuff
+ double rpESum = 0.;
+ double rpPxSum = 0.;
+ double rpPySum = 0.;
+ double rpPzSum = 0.;
+ for(ReconstructedParticle rp : rpList)
+ {
+ aida.cloud1D("PFA RP Energy per RP").fill(rp.getEnergy());
+ rpESum += rp.getEnergy();
+ Hep3Vector rpmom = rp.getMomentum();
+ rpPxSum += rpmom.x();
+ rpPySum += rpmom.y();
+ rpPzSum += rpmom.z();
+ }
+ double rpPSum = Math.sqrt(rpPxSum*rpPxSum+rpPySum*rpPySum+rpPzSum*rpPzSum);
+ aida.cloud1D("PFA RP Diff ESum").fill(rpESum-qqbarE);
+// aida.cloud1D("PFA RP PSum").fill(rpPSum);
+ aida.cloud1D("PFA RP Diff Mass").fill(Math.sqrt(rpESum*rpESum-rpPSum*rpPSum)-Zmass);
+
+ // 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;
+ 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)
+ {
+ aida.cloud1D("Perfect Jet Energy").fill(perfjet.getEnergy());
+ aida.cloud1D("Perfect Jet Mass").fill(perfjet.getMass());
+// aida.cloud1D("Number of particles in Jet").fill(perfjet.getParticles().size());
+ 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());
+ // force jet mass 0
+// jESum += jp[j];
+// aida.cloud1D("Perfect Jet Energy zero mass").fill(jp[j]);
+ jpz[j] = j3v.z();
+ j++;
+ }
+ double jE2 = jESum*jESum;
+ double jP2 = jPXSum*jPXSum+jPYSum*jPYSum+jPZSum*jPZSum;
+ djm = Math.sqrt(jE2-jP2);
+ aida.cloud1D("Z Mass").fill(Zmass);
+ aida.cloud1D("Perfect Dijet Mass").fill(djm);
+ aida.cloud1D("Perfect Dijet Energy").fill(Math.sqrt(jE2));
+ aida.cloud1D("Diff Perfect Dijet Mass").fill(djm-Zmass);
+ }
+
+ // 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;
+ 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)
+ {
+ aida.cloud1D("PFA Jet Energy").fill(pfajet.getEnergy());
+ aida.cloud1D("PFA Jet Mass").fill(pfajet.getMass());
+ 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());
+ // force jet mass 0
+// rjESum += jrp[k];
+// aida.cloud1D("PFA Jet Energy zero mass").fill(jrp[k]);
+ jrpz[k] = rj3v.z();
+ k++;
+ }
+ double rjE2 = rjESum*rjESum;
+ double rjP2 = rjPXSum*rjPXSum+rjPYSum*rjPYSum+rjPZSum*rjPZSum;
+ rdjm = Math.sqrt(rjE2-rjP2);
+ aida.cloud1D("PFA Dijet Mass").fill(rdjm);
+ aida.cloud1D("PFA Dijet Energy").fill(Math.sqrt(rjE2));
+ aida.cloud1D("Diff PFA Dijet Mass").fill(rdjm-Zmass);
+ aida.cloud1D("Difference PFAJet and PerfJet DiJet Mass").fill(rdjm-djm);
+ }
+
+ // now CS Scin cluster jets
+ double sjESum = 0;
+ double sjPXSum = 0;
+ double sjPYSum = 0;
+ double sjPZSum = 0;
+ double sdjm = 0;
+ double[] jsrpz = new double[2];
+ double[] jsrp = new double[2];
+ int ks = 0;
+ aida.cloud1D("Number of CS Scin Jets").fill(csscinjets.size());
+ if (csscinjets.size() == 2)
+ {
+ // calculate dijet mass - don't we have something that does this?
+ for (ReconstructedParticle cssjet : csscinjets)
+ {
+ aida.cloud1D("CS S Jet Energy").fill(cssjet.getEnergy());
+ aida.cloud1D("CS S Jet Mass").fill(cssjet.getMass());
+ sjESum += cssjet.getEnergy();
+ Hep3Vector sj3v = cssjet.getMomentum();
+ sjPXSum += sj3v.x();
+ sjPYSum += sj3v.y();
+ sjPZSum += sj3v.z();
+ jsrp[ks] = Math.sqrt(sj3v.x()*sj3v.x()+sj3v.y()*sj3v.y()+sj3v.z()*sj3v.z());
+ // force jet mass 0
+// sjESum += jsrp[ks];
+// aida.cloud1D("CS S Jet Energy zero mass").fill(jrp[ks]);
+ jsrpz[ks] = sj3v.z();
+ ks++;
+ }
+ double sjE2 = sjESum*sjESum;
+ double sjP2 = sjPXSum*sjPXSum+sjPYSum*sjPYSum+sjPZSum*sjPZSum;
+ sdjm = Math.sqrt(sjE2-sjP2);
+ aida.cloud1D("CS S Dijet Mass").fill(sdjm);
+ aida.cloud1D("CS S Dijet Energy").fill(Math.sqrt(sjE2));
+ aida.cloud1D("Diff CS S Dijet Mass").fill(sdjm-Zmass);
+ aida.cloud1D("Difference CSSJet and PerfJet DiJet Mass").fill(sdjm-djm);
+ }
+ }
+}