Commit in lcsim-contrib/src/main/java/org/lcsim/contrib/SteveMagill on MAIN | |||
PFADRSelect.java | +362 | added 1.1 |
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); + } + } +}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1