Author: [log in to unmask]
Date: Thu Sep 1 04:26:03 2016
New Revision: 4477
Log: (empty)
Added:
java/trunk/users/src/main/java/org/hps/users/spaul/fee_old/
- copied from r4472, java/trunk/users/src/main/java/org/hps/users/spaul/feecc/
Removed:
java/trunk/users/src/main/java/org/hps/users/spaul/feecc/
Modified:
java/trunk/users/src/main/java/org/hps/users/spaul/fee/Cleanup.java
java/trunk/users/src/main/java/org/hps/users/spaul/fee/FeeHistogramDriver.java
java/trunk/users/src/main/java/org/hps/users/spaul/moller/MollerBeamtiltAnalysis.java
java/trunk/users/src/main/java/org/hps/users/spaul/trident/TridentHistogramDriver.java
Modified: java/trunk/users/src/main/java/org/hps/users/spaul/fee/Cleanup.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/fee/Cleanup.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/fee/Cleanup.java Thu Sep 1 04:26:03 2016
@@ -5,24 +5,34 @@
import java.util.ArrayList;
import java.util.List;
+import org.hps.recon.tracking.TrackUtils;
import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Track;
public class Cleanup {
- public static void purgeDuplicates(List<ReconstructedParticle> particles) {
+ public static void purgeParticlesWithSameCluster(List<ReconstructedParticle> particles) {
ArrayList<ReconstructedParticle> trashcan = new ArrayList();
outer : for(ReconstructedParticle p1 : particles){
if(trashcan.contains(p1))
continue;
+ if(p1.getClusters().size() == 0){
+ trashcan.add(p1);
+ continue;
+ }
for(ReconstructedParticle p2 : particles){
if(p1 == p2) continue;
if(trashcan.contains(p2))
continue;
+ if(p2.getClusters().size() == 0)
+ {
+ trashcan.add(p2);
+ continue;
+ }
if((p1.getType()>31)^(p2.getType()>31)){ //one is GBL, the other is seed.
continue;
}
-
- if(p1.getEnergy() == p2.getEnergy()){ //tracks matched to same cluster
+ if(p1.getClusters().get(0) == p2.getClusters().get(0)){ //tracks matched to same cluster
ReconstructedParticle loser = getWorseParticle(p1, p2);
trashcan.add(loser);
}
@@ -36,18 +46,41 @@
}
static ReconstructedParticle getWorseParticle(ReconstructedParticle p1, ReconstructedParticle p2){
- if(p1.getTracks().get(0).getChi2() < p2.getTracks().get(0).getChi2())
+ //first check if the tracks come from the target:
+
+ if(p1.getTracks().size() == 0)
+ return p1;
+ if(p2.getTracks().size() == 0)
+ return p2;
+
+ Track t1 = p1.getTracks().get(0);
+ Track t2 = p2.getTracks().get(0);
+
+
+ //first see if there is some clue that one of the tracks is better than the other:
+ if(Math.abs(TrackUtils.getDoca(t1))>1 && Math.abs(TrackUtils.getDoca(t2))<1)
+ return p2;
+ if(Math.abs(TrackUtils.getDoca(t2))>1 && Math.abs(TrackUtils.getDoca(t1))<1)
+ return p1;
+
+ //if they are both ok, use the one with the track that extrapolates closes to the cluster.
+ if(Math.hypot(FeeHistogramDriver.getDx(p1), FeeHistogramDriver.getDy(p1))
+ < Math.hypot(FeeHistogramDriver.getDx(p2), FeeHistogramDriver.getDy(p2)))
return p2;
return p1;
+
+
+ //if they are both ok, then just use whichever one has the smaller chi2
+ //if(p1.getTracks().get(0).getChi2() < p2.getTracks().get(0).getChi2())
+ //return p2;
+ //return p1;
+
}
- static double angle(Hep3Vector v1, Hep3Vector v2){
- return Math.acos(v1.x()*v2.x() + v1.y()*v2.y() + v1.z()*v2.z())
- /(v1.magnitude()*v2.magnitude());
- }
+
}
Modified: java/trunk/users/src/main/java/org/hps/users/spaul/fee/FeeHistogramDriver.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/fee/FeeHistogramDriver.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/fee/FeeHistogramDriver.java Thu Sep 1 04:26:03 2016
@@ -24,338 +24,399 @@
public class FeeHistogramDriver extends Driver{
-
+
double d0cut = 1;
double z0cut = .5;
-
+
//determined by beam energy
- double eMin, eMax, beamEnergy, seedEnergyCut;
- double pMin, pMax;
-
-
- double beamTiltX = .0295;
- double beamTiltY = -.0008;
-
- double maxChi2 = 50;
- //maximum difference between the reconstructed energy and momentum
- //double maxdE = .5;
- int nMin = 3;
-
- IHistogram1D theta;
- IHistogram1D thetaInRange;
- IHistogram2D thetaVsPhi;
- IHistogram2D thetaVsPhiInRange;
- IHistogram2D d0VsZ0, d0VsZ0_top, d0VsZ0_bottom;
- IHistogram1D d0;
- IHistogram1D z0;
- IHistogram1D timeHist;
- IHistogram1D chi2Hist;
-
-
- @Override
- public void detectorChanged(Detector det){
- BeamEnergyCollection beamEnergyCollection =
- this.getConditionsManager().getCachedConditions(BeamEnergyCollection.class, "beam_energies").getCachedData();
- beamEnergy = beamEnergyCollection.get(0).getBeamEnergy();
- eMin = .75*beamEnergy;
- eMax = 1.15*beamEnergy;
-
- pMin = .75*beamEnergy;
- pMax = 1.15*beamEnergy;
-
- seedEnergyCut = .38*beamEnergy;
- setupHistograms();
-
- }
-
- AIDA aida = AIDA.defaultInstance();
-
-
- private IHistogram1D thetaTopOnly;
-
-
- private IHistogram1D thetaTopOnlyInRange;
-
-
- private IHistogram1D thetaBottomOnlyInRange;
-
-
- private IHistogram1D energy, seedEnergy, energyTop, seedEnergyTop, energyBottom, seedEnergyBottom;
- private IHistogram2D clusterVsSeedEnergy, clusterVsSeedEnergyTop, clusterVsSeedEnergyBottom;
-
- private IHistogram1D thetaBottomOnly;
-
- private IHistogram1D charge;
-
-
- void setupHistograms(){
-
-
- energy = aida.histogram1D("cluster energy", 100, 0, 1.5*beamEnergy);
- energyTop = aida.histogram1D("cluster energy top", 100, 0, 1.5*beamEnergy);
- energyBottom = aida.histogram1D("cluster energy bottom", 100, 0, 1.5*beamEnergy);
-
- seedEnergy = aida.histogram1D("seed energy", 100, 0, beamEnergy);
- seedEnergyTop = aida.histogram1D("seed energy top", 100, 0, beamEnergy);
- seedEnergyBottom = aida.histogram1D("seed energy bottom", 100, 0, beamEnergy);
-
- clusterVsSeedEnergy = aida.histogram2D("cluster vs seed energy", 100, 0, 1.5*beamEnergy, 100, 0, beamEnergy);
- clusterVsSeedEnergyTop = aida.histogram2D("cluster vs seed energy top", 100, 0, 1.5*beamEnergy, 100, 0, beamEnergy);
- clusterVsSeedEnergyBottom = aida.histogram2D("cluster vs seed energy bottom", 100, 0, 1.5*beamEnergy, 100, 0, beamEnergy);
-
-
- int nBinsPhi = 400;
- thetaVsPhi = aida.histogram2D("theta vs phi", nBinsPhi, 0, .2, 628, -3.14, 3.14);
- thetaVsPhiInRange = aida.histogram2D("theta vs phi in range", nBinsPhi, 0, .2, 628, -3.14, 3.14);
-
- theta = aida.histogram1D("theta", nBinsPhi, 0, .2);
- thetaInRange = aida.histogram1D("theta in range", nBinsPhi, 0, .2);
-
- thetaTopOnly = aida.histogram1D("theta top", nBinsPhi, 0, .2);
- thetaTopOnlyInRange = aida.histogram1D("theta top in range", nBinsPhi, 0, .2);
-
- thetaBottomOnly = aida.histogram1D("theta bottom", nBinsPhi, 0, .2);
- thetaBottomOnlyInRange = aida.histogram1D("theta bottom in range", nBinsPhi, 0, .2);
- uxVsUy = aida.histogram2D("ux vs uy", 300, -.16, .24, 300, -.2, .2);
- uxVsUyInRange = aida.histogram2D("ux vs uy in range", 300, -.16, .24, 300, -.2, .2);
- d0VsZ0 = aida.histogram2D("d0 vs z0", 100, -5, 5, 100, -5, 5);
- d0 = aida.histogram1D("d0", 100, -5, 5);
- z0 = aida.histogram1D("z0", 100, -5, 5);
-
- d0VsZ0_top = aida.histogram2D("d0 vs z0 top", 100, -5, 5, 100, -5, 5);
- d0VsZ0_bottom = aida.histogram2D("d0 vs z0 bottom", 100, -5, 5, 100, -5, 5);
-
-
- pHist = aida.histogram1D("momentum", 100, 0, 1.5*beamEnergy);
-
- charge = aida.histogram1D("charge", 6, -3, 3);
-
-
- timeHist = aida.histogram1D("cluster time", 400, 0, 400);
- chi2Hist = aida.histogram1D("chi2", 200, 0, 200);
- clustSize = aida.histogram1D("cluster size", 20, 0, 20);
-
- z0VsTanLambda = aida.histogram2D("z0 vs tan lambda", 100, -5, 5, 100, -.1, .1);
- z0VsChi2 = aida.histogram2D("z0 vs chi2", 100, -5, 5, 100, 0, 100);
-
- }
-
-
- IHistogram2D uxVsUy, uxVsUyInRange;
-
- IHistogram1D pHist;
- IHistogram1D clustSize;
-
- @Override
- public void process(EventHeader event){
-
-
-
- //removed this criterion from this driver.
- //instead use a trigger filter driver preceding this driver in the steering file.
- //
-
- /*for (GenericObject gob : event.get(GenericObject.class,"TriggerBank"))
- {
- if (!(AbstractIntData.getTag(gob) == TIData.BANK_TAG)) continue;
- TIData tid = new TIData(gob);
- if (!tid.isSingle1Trigger())
- {
- return;
- }
- }*/
-
- List<ReconstructedParticle> particles = event.get(ReconstructedParticle.class, "FinalStateParticles");
- Cleanup.purgeDuplicates(particles);
- for(ReconstructedParticle p : particles){
- processParticle(p);
- }
- }
-
-
- private void processParticle(ReconstructedParticle part) {
-
- //reject duplicate particles that use seed tracks instead of GBL
- if(!TrackType.isGBL(part.getType()) && part.getTracks().size() >0)
- return;
-
- //reject particles that have no cluster in the Ecal
- if(part.getClusters().size() == 0)
- return;
- Cluster c = part.getClusters().get(0);
-
- //cluster time
- double time = c.getCalorimeterHits().get(0).getTime();
- timeHist.fill(time);
-
- if(!(time>minClustTime && time <maxClustTime))
- return;
-
-
- //find seed hit energy
- double seedEnergy = 0;
- for(CalorimeterHit hit : c.getCalorimeterHits()){
- if(hit.getCorrectedEnergy() > seedEnergy)
- seedEnergy = hit.getCorrectedEnergy();
- }
- double energy = c.getEnergy();
-
-
- //these are the histograms of seed and cluster energy before
- //the cuts on these variables
- this.seedEnergy.fill(seedEnergy);
- this.energy.fill(energy);
- this.clusterVsSeedEnergy.fill(energy, seedEnergy);
-
- if(part.getMomentum().y()> 0){
- this.seedEnergyTop.fill(seedEnergy);
- this.energyTop.fill(energy);
- this.clusterVsSeedEnergyTop.fill(energy, seedEnergy);
- } else {
- this.seedEnergyBottom.fill(seedEnergy);
- this.energyBottom.fill(energy);
- this.clusterVsSeedEnergyBottom.fill(energy, seedEnergy);
- }
-
-
- //seed energy cut
- if(seedEnergy < seedEnergyCut)
- return;
-
- //cluster energy cut
- if(energy > eMax || energy< eMin)
- return;
-
- clustSize.fill(c.getCalorimeterHits().size());
- if(c.getCalorimeterHits().size() < nMin)
- return;
-
-
-
- charge.fill(part.getCharge());
- if(part.getCharge() != -1)
- return;
-
- if(part.getTracks().size() == 0)
- return;
-
-
-
-
-
- Hep3Vector p = part.getMomentum();
-
- double pmag = part.getMomentum().magnitude();
- pHist.fill(pmag);
- if(pmag > pMax || pmag < pMin)
- return;
-
- double d = TrackUtils.getDoca(part.getTracks().get(0));
- d0.fill(d);
- double z = TrackUtils.getZ0(part.getTracks().get(0));
- z0.fill(z);
- d0VsZ0.fill(d, z);
-
- if(p.y() > 0)
- d0VsZ0_top.fill(d,z);
- else
- d0VsZ0_bottom.fill(d,z);
-
- z0VsTanLambda.fill(z, TrackUtils.getTanLambda(part.getTracks().get(0)));
- z0VsChi2.fill(z, part.getTracks().get(0).getChi2());
-
- if(abs(TrackUtils.getDoca(part.getTracks().get(0)))> d0cut)
- return;
- if(abs(TrackUtils.getZ0(part.getTracks().get(0)))> z0cut)
- return;
-
-
-
- chi2Hist.fill(part.getTracks().get(0).getChi2());
- if(part.getTracks().get(0).getChi2()>maxChi2)
- return;
-
-
-
-/*
+ double eMin, eMax, beamEnergy, seedEnergyCut;
+ public double getSeedEnergyCutFrac() {
+ return seedEnergyCutFrac;
+ }
+
+ public void setSeedEnergyCutFrac(double seedEnergyCut) {
+ this.seedEnergyCutFrac = seedEnergyCut;
+ }
+
+ double pMin, pMax;
+
+
+ double beamTiltX = .0295;
+ double beamTiltY = -.0008;
+
+ double maxChi2 = 50;
+ //maximum difference between the reconstructed energy and momentum
+ //double maxdE = .5;
+ int nMin = 3;
+
+ IHistogram1D theta;
+ IHistogram1D thetaInRange;
+ IHistogram2D[] thetaVsPhi;
+ IHistogram2D[] thetaVsPhiNoTrackQualityCuts;
+ IHistogram2D[] thetaVsPhiInRange;
+ IHistogram2D d0VsZ0, d0VsZ0_top, d0VsZ0_bottom;
+ IHistogram1D d0;
+ IHistogram1D z0;
+ IHistogram1D timeHist;
+ IHistogram1D chi2Hist;
+ private double seedEnergyCutFrac = .38;
+
+
+ @Override
+ public void detectorChanged(Detector det){
+ BeamEnergyCollection beamEnergyCollection =
+ this.getConditionsManager().getCachedConditions(BeamEnergyCollection.class, "beam_energies").getCachedData();
+ beamEnergy = beamEnergyCollection.get(0).getBeamEnergy();
+ eMin = .75*beamEnergy;
+ eMax = 1.15*beamEnergy;
+
+ pMin = .75*beamEnergy;
+ pMax = 1.15*beamEnergy;
+
+ seedEnergyCut = seedEnergyCutFrac *beamEnergy;
+ setupHistograms();
+
+ }
+
+ AIDA aida = AIDA.defaultInstance();
+
+
+ private IHistogram1D thetaTopOnly;
+
+
+ private IHistogram1D thetaTopOnlyInRange;
+
+
+ private IHistogram1D thetaBottomOnlyInRange;
+
+
+ private IHistogram1D energy, seedEnergy, energyTop, seedEnergyTop, energyBottom, seedEnergyBottom;
+ private IHistogram2D clusterVsSeedEnergy, clusterVsSeedEnergyTop, clusterVsSeedEnergyBottom;
+
+ private IHistogram1D thetaBottomOnly;
+
+ private IHistogram1D charge;
+ private IHistogram1D clustSizeGTP;
+
+
+ void setupHistograms(){
+
+
+ energy = aida.histogram1D("cluster energy", 100, 0, 1.5*beamEnergy);
+ energyTop = aida.histogram1D("cluster energy top", 100, 0, 1.5*beamEnergy);
+ energyBottom = aida.histogram1D("cluster energy bottom", 100, 0, 1.5*beamEnergy);
+
+ seedEnergy = aida.histogram1D("seed energy", 100, 0, beamEnergy);
+ seedEnergyTop = aida.histogram1D("seed energy top", 100, 0, beamEnergy);
+ seedEnergyBottom = aida.histogram1D("seed energy bottom", 100, 0, beamEnergy);
+
+ clusterVsSeedEnergy = aida.histogram2D("cluster vs seed energy", 100, 0, 1.5*beamEnergy, 100, 0, beamEnergy);
+ clusterVsSeedEnergyTop = aida.histogram2D("cluster vs seed energy top", 100, 0, 1.5*beamEnergy, 100, 0, beamEnergy);
+ clusterVsSeedEnergyBottom = aida.histogram2D("cluster vs seed energy bottom", 100, 0, 1.5*beamEnergy, 100, 0, beamEnergy);
+
+
+ int nBinsPhi = 400;
+
+ String[] regionNames = {"", " r1", " r2", " r3", " r4"};
+
+ thetaVsPhi = new IHistogram2D[regionNames.length];
+ thetaVsPhiInRange= new IHistogram2D[regionNames.length];
+ thetaVsPhiNoTrackQualityCuts = new IHistogram2D[regionNames.length];
+ thetaVsPhiChi2Cut= new IHistogram2D[regionNames.length];
+ thetaVsPhiTrackExtrapCut = new IHistogram2D[regionNames.length];
+ thetaVsPhiMomentumCut = new IHistogram2D[regionNames.length];
+ for(int i = 0; i< regionNames.length; i++){
+ thetaVsPhi[i] = aida.histogram2D("theta vs phi" + regionNames[i], nBinsPhi, 0, .2, 628, -3.14, 3.14);
+ thetaVsPhiInRange[i] = aida.histogram2D("theta vs phi in range" + regionNames[i], nBinsPhi, 0, .2, 628, -3.14, 3.14);
+ thetaVsPhiNoTrackQualityCuts[i] = aida.histogram2D("theta vs phi no track quality cuts" + regionNames[i], nBinsPhi, 0, .2, 628, -3.14, 3.14);
+ thetaVsPhiChi2Cut[i] = aida.histogram2D("theta vs phi chi2 cut" + regionNames[i], nBinsPhi, 0, .2, 628, -3.14, 3.14);
+ thetaVsPhiTrackExtrapCut[i] = aida.histogram2D("theta vs phi track extrap cut" + regionNames[i], nBinsPhi, 0, .2, 628, -3.14, 3.14);
+ thetaVsPhiMomentumCut[i] = aida.histogram2D("theta vs phi chi2 momentum cut" + regionNames[i], nBinsPhi, 0, .2, 628, -3.14, 3.14);
+ }
+ theta = aida.histogram1D("theta", nBinsPhi, 0, .2);
+ thetaInRange = aida.histogram1D("theta in range", nBinsPhi, 0, .2);
+
+ thetaTopOnly = aida.histogram1D("theta top", nBinsPhi, 0, .2);
+ thetaTopOnlyInRange = aida.histogram1D("theta top in range", nBinsPhi, 0, .2);
+
+ thetaBottomOnly = aida.histogram1D("theta bottom", nBinsPhi, 0, .2);
+ thetaBottomOnlyInRange = aida.histogram1D("theta bottom in range", nBinsPhi, 0, .2);
+ uxVsUy = aida.histogram2D("ux vs uy", 300, -.16, .24, 300, -.2, .2);
+ uxVsUyInRange = aida.histogram2D("ux vs uy in range", 300, -.16, .24, 300, -.2, .2);
+ d0VsZ0 = aida.histogram2D("d0 vs z0", 100, -5, 5, 100, -5, 5);
+ d0 = aida.histogram1D("d0", 100, -5, 5);
+ z0 = aida.histogram1D("z0", 100, -5, 5);
+
+ d0VsZ0_top = aida.histogram2D("d0 vs z0 top", 100, -5, 5, 100, -5, 5);
+ d0VsZ0_bottom = aida.histogram2D("d0 vs z0 bottom", 100, -5, 5, 100, -5, 5);
+
+
+ pHist = aida.histogram1D("momentum", 100, 0, 1.5*beamEnergy);
+
+ charge = aida.histogram1D("charge", 6, -3, 3);
+
+
+ timeHist = aida.histogram1D("cluster time", 400, 0, 400);
+ chi2Hist = aida.histogram1D("chi2", 200, 0, 200);
+ chi2RedHist = aida.histogram1D("chi2 red", 200, 0, 40);
+ clustSize = aida.histogram1D("cluster size", 20, 0, 20);
+ clustSizeGTP = aida.histogram1D("cluster size gtp (>55% beam energy)", 20, 0, 20);
+
+ z0VsTanLambda = aida.histogram2D("z0 vs tan lambda", 100, -5, 5, 100, -.1, .1);
+ z0VsChi2 = aida.histogram2D("z0 vs chi2", 100, -5, 5, 100, 0, 100);
+ nPassCutsPerEvent = aida.histogram1D("n pass cuts", 10, 0, 10);
+ nPassCutsPerEventFiducial = aida.histogram1D("n pass cuts fiducial", 10, 0, 10);
+ thetaVsMomentum = aida.histogram2D("theta vs energy", 100, 0, .2, 100, 0, 1.2*beamEnergy);
+ thetaVsMomentumFid = aida.histogram2D("theta vs energy fiducial", 100, 0, .2, 100, 0, 1.2*beamEnergy);
+ }
+
+
+ IHistogram2D uxVsUy, uxVsUyInRange;
+
+ IHistogram1D pHist;
+ IHistogram1D clustSize;
+
+ IHistogram1D nPassCutsPerEvent;
+ IHistogram1D nPassCutsPerEventFiducial;
+ private IHistogram1D chi2RedHist;
+ private IHistogram2D thetaVsPhiChi2Cut[];
+ private IHistogram2D thetaVsPhiMomentumCut[];
+ private IHistogram2D thetaVsPhiTrackExtrapCut[];
+ @Override
+ public void process(EventHeader event){
+
+ fillHistogramsGTP(event);
+
+ //int col = getColumn(event);
+
+ List<ReconstructedParticle> particles = event.get(ReconstructedParticle.class, "FinalStateParticles");
+ Cleanup.purgeParticlesWithSameCluster(particles);
+ int nPassFidCuts = 0, nPassCuts = 0;
+ for(ReconstructedParticle part : particles){
+
+
+ //reject duplicate particles that use seed tracks instead of GBL
+ if(!TrackType.isGBL(part.getType()) && part.getTracks().size() >0)
+ continue;
+
+ //reject particles that have no cluster in the Ecal
+ if(part.getClusters().size() == 0)
+ continue;
+ Cluster c = part.getClusters().get(0);
+
+ //cluster time
+ double time = c.getCalorimeterHits().get(0).getTime();
+ timeHist.fill(time);
+
+ if(!(time>minClustTime && time <maxClustTime))
+ continue;
+
+
+ //find seed hit energy
+ double seedEnergy = 0;
+ int col = 0;
+ for(CalorimeterHit hit : c.getCalorimeterHits()){
+ if(hit.getCorrectedEnergy() > seedEnergy){
+ seedEnergy = hit.getCorrectedEnergy();
+ col = hit.getIdentifierFieldValue("ix");
+ }
+ }
+ double energy = c.getEnergy();
+
+
+ //these are the histograms of seed and cluster energy before
+ //the cuts on these variables
+ this.seedEnergy.fill(seedEnergy);
+ this.energy.fill(energy);
+ this.clusterVsSeedEnergy.fill(energy, seedEnergy);
+
+ if(part.getMomentum().y()> 0){
+ this.seedEnergyTop.fill(seedEnergy);
+ this.energyTop.fill(energy);
+ this.clusterVsSeedEnergyTop.fill(energy, seedEnergy);
+ } else {
+ this.seedEnergyBottom.fill(seedEnergy);
+ this.energyBottom.fill(energy);
+ this.clusterVsSeedEnergyBottom.fill(energy, seedEnergy);
+ }
+
+
+ //seed energy cut
+ if(seedEnergy < seedEnergyCut)
+ continue;
+
+ //cluster energy cut
+ if(energy > eMax || energy< eMin)
+ continue;
+
+ clustSize.fill(c.getCalorimeterHits().size());
+ if(c.getCalorimeterHits().size() < nMin)
+ continue;
+
+
+
+ charge.fill(part.getCharge());
+ if(part.getCharge() != -1)
+ continue;
+
+ if(part.getTracks().size() == 0)
+ continue;
+
+
+ Hep3Vector p = part.getMomentum();
+ double px = p.x(), py = p.y(), pz = p.z();
+
+ double cx = cos(beamTiltX);
+ double sy = sin(beamTiltY);
+ double cy = cos(beamTiltY);
+ double sx = sin(beamTiltX);
+
+ double pxtilt = px*cx -pz*sx;
+ double pytilt = -py*sx*sy + py*cy -pz*sy*cx;
+ double pztilt = px*cy*sx + py*sy +pz*cy*cx;
+
+
+ double theta = atan(hypot(pxtilt, pytilt)/pztilt);
+ double phi =atan2(pytilt, pxtilt);
+ fill(thetaVsPhiNoTrackQualityCuts,theta, phi, col);
+
+ if(abs(getDx(part))>trackExtrapCut || abs(getDy(part))> trackExtrapCut)
+ continue;
+ fill(thetaVsPhiTrackExtrapCut, theta, phi, col);
+
+ chi2Hist.fill(part.getTracks().get(0).getChi2());
+ chi2RedHist.fill(part.getTracks().get(0).getChi2()/(2*part.getTracks().get(0).getTrackerHits().size()-5));
+
+ if(part.getTracks().get(0).getChi2()>maxChi2)
+ continue;
+ fill(thetaVsPhiChi2Cut, theta, phi, col);
+
+ double pmag = part.getMomentum().magnitude();
+ pHist.fill(pmag);
+ if(pmag > pMax || pmag < pMin)
+ continue;
+ fill(thetaVsPhiMomentumCut, theta, phi, col);
+
+ double d = TrackUtils.getDoca(part.getTracks().get(0));
+ d0.fill(d);
+ double z = TrackUtils.getZ0(part.getTracks().get(0));
+ z0.fill(z);
+ d0VsZ0.fill(d, z);
+
+ if(p.y() > 0)
+ d0VsZ0_top.fill(d,z);
+ else
+ d0VsZ0_bottom.fill(d,z);
+
+ z0VsTanLambda.fill(z, TrackUtils.getTanLambda(part.getTracks().get(0)));
+ z0VsChi2.fill(z, part.getTracks().get(0).getChi2());
+
+ if(abs(TrackUtils.getDoca(part.getTracks().get(0)))> d0cut)
+ continue;
+ if(abs(TrackUtils.getZ0(part.getTracks().get(0)))> z0cut)
+ continue;
+
+
+
+
+
+
+
+ /*
double px = p.x(), pz = p.z();
double pxtilt = px*cos(beamTiltX)-pz*sin(beamTiltX);
double py = p.y();
double pztilt = pz*cos(beamTiltX)+px*sin(beamTiltX);
-
+
double pytilt = py*cos(beamTiltY)-pztilt*sin(beamTiltY);
pztilt = pz*cos(beamTiltY) + pytilt*sin(beamTiltY);
-*/
- double px = p.x(), py = p.y(), pz = p.z();
-
- double cx = cos(beamTiltX);
- double sy = sin(beamTiltY);
- double cy = cos(beamTiltY);
- double sx = sin(beamTiltX);
-
- double pxtilt = px*cx -pz*sx;
- double pytilt = -py*sx*sy + py*cy -pz*sy*cx;
- double pztilt = px*cy*sx + py*sy +pz*cy*cx;
-
-
- double theta = atan(hypot(pxtilt, pytilt)/pztilt);
- double phi =atan2(pytilt, pxtilt);
-
- uxVsUy.fill(pxtilt/pztilt, pytilt/pztilt);
- thetaVsPhi.fill(theta, phi);
- this.theta.fill(theta);
- if(phi > 0)
- thetaTopOnly.fill(theta);
- else
- thetaBottomOnly.fill(theta);
- if(cb.inRange(theta, phi)){
- thetaVsPhiInRange.fill(theta, phi);
- thetaInRange.fill(theta);
- if(phi > 0)
- thetaTopOnlyInRange.fill(theta);
- else
- thetaBottomOnlyInRange.fill(theta);
- uxVsUyInRange.fill(pxtilt/pztilt, pytilt/pztilt);
- }
- }
-
- public void setBinning(String binning){
- System.out.println("binning scheme:\n" + binning);
- this.cb = new CustomBinning(binning);
- }
- /*
- @Override
- protected void endOfData(){
- createRenormalizationHistogram();
- }
-
- private void createRenormalizationHistogram(){
-
- //the purpose of this histogram is so that there is a simple histogram
- //that each bin in the theta histograms can be normalized to the
- //angular range in phi that I am using for the cut inside of that theta bin.
-
- int xbins = thetaVsPhi.xAxis().bins();
- double xmin = thetaVsPhi.xAxis().lowerEdge();
- double xmax = thetaVsPhi.xAxis().upperEdge();
- double ymin = thetaVsPhi.xAxis().lowerEdge();
- double ymax = thetaVsPhi.xAxis().upperEdge();
-
- IHistogram1D phiHist = aida.histogram1D("phiWidth", xbins, xmin, xmax);
-
- for(int i = 0; i< cb.nTheta; i++){
- double phiWidth = cb.getPhiWidth(i);
- for(int j = 0; j<phiWidth*1000; j++){
- phiHist.fill(cb.thetaMin[i]+.001);
- }
- }
-
-
- }*/
-
- CustomBinning cb;
- double minClustTime = 40;
- double maxClustTime = 50;
+ */
+
+
+ uxVsUy.fill(pxtilt/pztilt, pytilt/pztilt);
+ fill(thetaVsPhi, theta, phi, col);
+ this.theta.fill(theta);
+ if(phi > 0)
+ thetaTopOnly.fill(theta);
+ else
+ thetaBottomOnly.fill(theta);
+ nPassCuts++;
+ thetaVsMomentum.fill(theta, pz);
+ if(cb.inRange(theta, phi)){
+ fill(thetaVsPhiInRange, theta, phi, col);
+ thetaInRange.fill(theta);
+ if(phi > 0)
+ thetaTopOnlyInRange.fill(theta);
+ else
+ thetaBottomOnlyInRange.fill(theta);
+ uxVsUyInRange.fill(pxtilt/pztilt, pytilt/pztilt);
+ nPassFidCuts++;
+ thetaVsMomentumFid.fill(theta, pz);
+ }
+ }
+ nPassCutsPerEvent.fill(nPassCuts);
+ nPassCutsPerEventFiducial.fill(nPassFidCuts);
+ }
+ //returns the column of the seed hit of the trigger cluster in the GTP
+ /*private int getColumn(EventHeader event) {
+ for (GenericObject gob : event.get(GenericObject.class,"TriggerBank"))
+ {
+ if (!(AbstractIntData.getTag(gob) == TIData.BANK_TAG)) continue;
+ TIData tid = new TIData(gob);
+
+ }
+
+
+ List<Cluster> l = event.get(Cluster.class, "EcalClusters");
+ double maxClusterEnergy = 0;
+ int col = 0;
+ for(Cluster c : l){
+ if(c.getEnergy()>maxClusterEnergy){
+ maxClusterEnergy = c.getEnergy();
+ double maxHitEnergy = 0;
+ for(CalorimeterHit hit : c.getCalorimeterHits()){
+ if(hit.getCorrectedEnergy() > maxHitEnergy){
+ col = hit.getIdentifierFieldValue("ix");
+ maxHitEnergy = hit.getCorrectedEnergy();
+ }
+ }
+ }
+ }
+
+ return col;
+ }*/
+
+ private void fill(IHistogram2D[] thetaVsPhiHist, double theta,
+ double phi, int col) {
+ thetaVsPhiHist[0].fill(theta, phi);
+ if(col <= -13 || col>=6)
+ thetaVsPhiHist[1].fill(theta, phi);
+ else if (col <= -9 || col >= 2)
+ thetaVsPhiHist[2].fill(theta, phi);
+ else if (col <= -7 || col >= -2)
+ thetaVsPhiHist[3].fill(theta, phi);
+ else
+ thetaVsPhiHist[4].fill(theta, phi);
+
+ }
+
+ private void fillHistogramsGTP(EventHeader event) {
+ if(event.hasCollection(Cluster.class, "EcalClusters"))
+ for(Cluster c : event.get(Cluster.class, "EcalClusters")){
+ double energy = c.getEnergy();
+ if(energy>beamEnergy*.55)
+ clustSizeGTP.fill(c.getSize());
+ }
+ }
+
+ public void setBinning(String binning){
+ System.out.println("binning scheme:\n" + binning);
+ this.cb = new CustomBinning(binning);
+ }
+
+ CustomBinning cb;
+ double minClustTime = 40;
+ double maxClustTime = 50;
public double getMinClustTime() {
@@ -379,8 +440,22 @@
public void setMaxClustTime(double maxClustTime) {
this.maxClustTime = maxClustTime;
}
-
+
IHistogram2D z0VsTanLambda;
IHistogram2D z0VsChi2;
+ IHistogram2D thetaVsMomentumFid;
+ IHistogram2D thetaVsMomentum;
+ public double trackExtrapCut = 10;
+ public static double getDx(ReconstructedParticle p){
+ double xc = p.getClusters().get(0).getPosition()[0];
+ double xt = TrackUtils.getTrackPositionAtEcal(p.getTracks().get(0)).x();
+ return xt-xc;
+ }
+
+ public static double getDy(ReconstructedParticle p) {
+ double yc = p.getClusters().get(0).getPosition()[1];
+ double yt = TrackUtils.getTrackPositionAtEcal(p.getTracks().get(0)).y();
+ return yt-yc;
+ }
}
Modified: java/trunk/users/src/main/java/org/hps/users/spaul/moller/MollerBeamtiltAnalysis.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/moller/MollerBeamtiltAnalysis.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/moller/MollerBeamtiltAnalysis.java Thu Sep 1 04:26:03 2016
@@ -1,14 +1,18 @@
package org.hps.users.spaul.moller;
+
+import static java.lang.Math.atan;
+import static java.lang.Math.atan2;
+import static java.lang.Math.cos;
+import static java.lang.Math.hypot;
+import static java.lang.Math.sin;
import java.util.List;
import hep.aida.IHistogram1D;
import hep.aida.IHistogram2D;
-import hep.aida.IPlotter;
-import hep.aida.IProfile1D;
+import hep.physics.vec.Hep3Vector;
import org.hps.recon.tracking.TrackType;
-import org.hps.users.spaul.StyleUtil;
import org.lcsim.event.EventHeader;
import org.lcsim.event.ReconstructedParticle;
import org.lcsim.event.Vertex;
@@ -28,9 +32,12 @@
public void process(EventHeader event){
for(int i = 0; i< mollerCollections.length; i++){
List<Vertex> mollers = event.get(Vertex.class, mollerCollections[i]);
+ Cleanup.purgeDuplicates(mollers);
+ int nPassCuts = 0;
for(Vertex v : mollers){
- if(!passesCuts(v))
+ if(!passesBradleysCuts(v))
continue;
+ nPassCuts ++;
ReconstructedParticle m = v.getAssociatedParticle();
ReconstructedParticle top;
ReconstructedParticle bottom;
@@ -55,75 +62,81 @@
double sum = m.getMomentum().z();//top.getMomentum().z()+bottom.getMomentum().z();
double mass = m.getMass();
- if(diff > -.05 && diff < .05){
- hpypz_mid[i].fill(pypz);
- hpxpz_mid[i].fill(pxpz);
- }
-
- if(diff > .2 && diff < .3){
- hpypz_topHighE[i].fill(pypz);
- hpxpz_topHighE[i].fill(pxpz);
- }
-
- if(diff > -.3 && diff < -.2){
- hpypz_botHighE[i].fill(pypz);
- hpxpz_botHighE[i].fill(pxpz);
- }
+
this.diff[i].fill(diff);
this.sum[i].fill(sum);
this.mass[i].fill(mass);
- pypz_vs_diff[i].fill(diff,pypz );
- pxpz_vs_diff[i].fill(diff, pxpz );
-
-
-
- pxpz_vs_sum[i].fill(sum, pxpz );
- pypz_vs_sum[i].fill(sum, pypz );
-
- pxpz_vs_mass[i].fill(mass, pxpz );
- pypz_vs_mass[i].fill(mass, pypz );
- timediff[i].fill(top.getClusters().get(0).getCalorimeterHits().get(0).getTime()
+ pypz_vs_diff[i].fill(pypz, diff );
+ pxpz_vs_diff[i].fill(pxpz, diff );
+
+
+
+ pxpz_vs_sum[i].fill(pxpz, sum );
+ pypz_vs_sum[i].fill(pypz, sum );
+
+ pxpz_vs_mass[i].fill(pxpz, mass);
+ pypz_vs_mass[i].fill(pypz, mass);
+ //only do this if both tracks have associated clusters:
+ if(top.getClusters().size()*bottom.getClusters().size()>0){
+ timediff[i].fill(top.getClusters().get(0).getCalorimeterHits().get(0).getTime()
-bottom.getClusters().get(0).getCalorimeterHits().get(0).getTime());
-
-
- double theta, phi;
- double pxpz_top = top.getMomentum().x()/top.getMomentum().z();
- double pypz_top = top.getMomentum().y()/top.getMomentum().z();
-
- theta = Math.hypot(pypz_top-(-.0008), pxpz_top-.0295);
- phi = Math.atan2(pypz_top-(-.0008), pxpz_top-.0295);
- thetaPhi[i].fill(theta, phi);
-
- double pxpz_bot = bottom.getMomentum().x()/bottom.getMomentum().z();
- double pypz_bot = bottom.getMomentum().y()/bottom.getMomentum().z();
-
- theta = Math.hypot(pypz_bot-(-.0008), pxpz_bot-.0295);
- phi = Math.atan2(pypz_bot-(-.0008), pxpz_bot-.0295);
- thetaPhi[i].fill(theta, phi);
-
- /*if(moreEnergetic.getMomentum().y() > 0)
- {
- pypz_tophighE.fill(pypz);
- pxpz_tophighE.fill(pxpz);
+ meanClusterTime[i].fill((top.getClusters().get(0).getCalorimeterHits().get(0).getTime()
+ +bottom.getClusters().get(0).getCalorimeterHits().get(0).getTime())/2.);
+ }
+
+ double [] res = getThetaAndPhi(top.getMomentum());
+ double theta1 = res[0];
+ double phi1 = res[1];
+ thetaPhi[i].fill(theta1, phi1);
+
+ res = getThetaAndPhi(bottom.getMomentum());
+ double theta2 = res[0];
+ double phi2 = res[1];
+ thetaPhi[i].fill(theta2, phi2);
+
+ vtxchi2[i].fill(v.getChi2());
+ trkchi2[i].fill(top.getTracks().get(0).getChi2());
+ trkchi2[i].fill(bottom.getTracks().get(0).getChi2());
+
+ vtx_xy[i].fill(v.getPosition().x(), v.getPosition().y());
+
+ double dPhi = phi1-phi2 + Math.PI;
+ while(dPhi > Math.PI)
+ dPhi-= 2*Math.PI;
+ while(dPhi< -Math.PI)
+ dPhi += 2*Math.PI;
+ dphi[i].fill(dPhi*180/Math.PI);
+ dphi_vs_dtheta[i].fill(dPhi*180/Math.PI, theta1-theta2);
+
+ theta_theta[i].fill(sin(theta1/2)*sin(theta2/2));
+
}
- if(moreEnergetic.getMomentum().y() < 0)
- {
- pypz_bottomhighE.fill(pypz);
- pxpz_bottomhighE.fill(pxpz);
- }*/
- }
+ this.nPassCuts[i].fill(nPassCuts);
}
}
- private double _maxVtxChi2 = 15;
- private double _maxTrkChi2 = 30;
+ //private double _maxVtxChi2 = 15;
+ //private double _maxTrkChi2 = 30;
+ private double _maxVtxChi2 = 30;
+ private double _maxTrkChi2 = 50;
private double _maxMass = .037;
private double _minMass = .030;
private double _minPz = 1.0;
private double _maxPz = 1.1;
- private boolean passesCuts(Vertex vertex){
+ private double _maxTrkPz = .8;
+ public double get_maxTrkPz() {
+ return _maxTrkPz;
+ }
+
+
+ public void set_maxTrkPz(double _maxTrkPz) {
+ this._maxTrkPz = _maxTrkPz;
+ }
+
+
+ private boolean passesCuts(Vertex vertex){
ReconstructedParticle m = vertex.getAssociatedParticle();
if(!TrackType.isGBL(m.getType()))
return false;
@@ -132,6 +145,7 @@
if(m.getMass() > _maxMass || m.getMass() < _minMass)
return false;
+
if(m.getParticles().get(0).getCharge() != -1
|| m.getParticles().get(1).getCharge() != -1 )
return false;
@@ -149,15 +163,61 @@
return false;
if(m.getParticles().get(1).getTracks().get(0).getChi2() > _maxTrkChi2)
return false;
+
+ if(m.getParticles().get(0).getMomentum().z()>_maxTrkPz )
+ return false;
+ if(m.getParticles().get(1).getMomentum().z()>_maxTrkPz )
+ return false;
+
+ //very loose cut around the nominal beam direction
+ if(Math.abs(m.getMomentum().x()/m.getMomentum().z()-.030)> .01)
+ return false;
+ if(Math.abs(m.getMomentum().y()/m.getMomentum().z()-.000)> .01)
+ return false;
+
return true;
}
-
- private IHistogram1D hpypz[], hpxpz[], diff[], sum[], mass[],
- hpypz_topHighE[], hpxpz_topHighE[],
- hpypz_botHighE[], hpxpz_botHighE[],
- hpypz_mid[], hpxpz_mid[];
-
- private IHistogram2D[] thetaPhi;
+
+ private boolean passesBradleysCuts(Vertex vertex){
+ ReconstructedParticle m = vertex.getAssociatedParticle();
+ if(!TrackType.isGBL(m.getType()))
+ return false;
+ if(m.getMomentum().z() < .95 || m.getMomentum().z() > 1.15)
+ return false;
+
+ //both particles must have clusters
+ if(m.getParticles().get(0).getClusters().size()==0 || m.getParticles().get(1).getClusters().size()==0)
+ return false;
+
+
+ ReconstructedParticle p1 = m.getParticles().get(0);
+ ReconstructedParticle p2 = m.getParticles().get(1);
+
+ //timing difference cut (2 ns)
+ if(Math.abs(p1.getClusters().get(0).getCalorimeterHits().get(0).getTime()-
+ p2.getClusters().get(0).getCalorimeterHits().get(0).getTime())>2)
+ return false;
+
+
+ double phi1 = getThetaAndPhi(p1.getMomentum())[1];
+ double phi2 = getThetaAndPhi(p2.getMomentum())[1];
+
+ double dPhi = phi1-phi2 + Math.PI;
+ while(dPhi>Math.PI)
+ dPhi-= 2*Math.PI;
+ while(dPhi<-Math.PI)
+ dPhi+= 2*Math.PI;
+
+ if(Math.abs(dPhi)>.35) //about 20 degrees
+ return false;
+ return true;
+ }
+
+ private IHistogram1D hpypz[], hpxpz[], diff[], sum[], mass[], trkchi2[], vtxchi2[], dphi[];
+
+ private IHistogram2D[] thetaPhi, dphi_vs_dtheta;
+
+ private IHistogram1D[] meanClusterTime;
public double getMaxVtxChi2() {
return _maxVtxChi2;
@@ -172,6 +232,8 @@
public double getMaxTrkChi2() {
return _maxTrkChi2;
}
+
+
public void setMaxTrkChi2(double _maxTrkChi2) {
@@ -219,9 +281,12 @@
}
- private IHistogram1D vtx_x[], vtx_y[], timediff[];
-
- private IProfile1D pxpz_vs_diff[], pypz_vs_diff[], pxpz_vs_sum[], pypz_vs_sum[],
+ private IHistogram2D vtx_xy[];
+ private IHistogram1D timediff[];
+ private IHistogram1D nPassCuts[], theta_theta[];
+
+
+ private IHistogram2D pxpz_vs_diff[], pypz_vs_diff[], pxpz_vs_sum[], pypz_vs_sum[],
pxpz_vs_mass[], pypz_vs_mass[];
@@ -237,31 +302,36 @@
hpypz = new IHistogram1D[3];
hpxpz = new IHistogram1D[3];
- hpypz_mid = new IHistogram1D[3];
- hpxpz_mid = new IHistogram1D[3];
- hpypz_topHighE = new IHistogram1D[3];
- hpxpz_topHighE = new IHistogram1D[3];
- hpypz_botHighE = new IHistogram1D[3];
- hpxpz_botHighE = new IHistogram1D[3];
-
-
- pxpz_vs_diff= new IProfile1D[3];
- pypz_vs_diff= new IProfile1D[3];
+
+
+ pxpz_vs_diff= new IHistogram2D[3];
+ pypz_vs_diff= new IHistogram2D[3];
diff= new IHistogram1D[3];
sum= new IHistogram1D[3];
- pxpz_vs_sum= new IProfile1D[3];
- pypz_vs_sum= new IProfile1D[3];
-
- pxpz_vs_mass= new IProfile1D[3];
- pypz_vs_mass= new IProfile1D[3];
+ pxpz_vs_sum= new IHistogram2D[3];
+ pypz_vs_sum= new IHistogram2D[3];
+
+ pxpz_vs_mass= new IHistogram2D[3];
+ pypz_vs_mass= new IHistogram2D[3];
//vtx_x= new IHistogram1D[3];
//vtx_y= new IHistogram1D[3];
mass= new IHistogram1D[3];
timediff= new IHistogram1D[3];
+ nPassCuts = new IHistogram1D[3];
+ vtx_xy = new IHistogram2D[3];
+ trkchi2 = new IHistogram1D[3];
+ vtxchi2 = new IHistogram1D[3];
+ meanClusterTime = new IHistogram1D[3];
+
+ dphi = new IHistogram1D[3];
+
+ theta_theta = new IHistogram1D[3];
+
+ dphi_vs_dtheta = new IHistogram2D[3];
for(int i = 0; i< 3; i++){
@@ -270,94 +340,57 @@
hpypz[i] = aida.histogram1D(mollerCollections[i]+"/"+"pypz", 60, -.005,.005);
hpxpz[i] = aida.histogram1D(mollerCollections[i]+"/"+"pxpz", 60, .025,.035);
-
- hpypz_mid[i] = aida.histogram1D(mollerCollections[i]+"/"+"pypz mid", 60, -.005,.005);
- hpxpz_mid[i] = aida.histogram1D(mollerCollections[i]+"/"+"pxpz mid", 60, .025,.035);
-
- hpypz_topHighE[i] = aida.histogram1D(mollerCollections[i]+"/"+"pypz top", 30, -.005,.005);
- hpxpz_topHighE[i] = aida.histogram1D(mollerCollections[i]+"/"+"pxpz top", 30, .025,.035);
-
- hpypz_botHighE[i] = aida.histogram1D(mollerCollections[i]+"/"+"pypz bot", 30, -.005,.005);
- hpxpz_botHighE[i] = aida.histogram1D(mollerCollections[i]+"/"+"pxpz bot", 30, .025,.035);
-
-
- pxpz_vs_diff[i] = aida.profile1D(mollerCollections[i]+"/"+"pxpz vs diff", 25, -.60, .60);
- pypz_vs_diff[i] = aida.profile1D(mollerCollections[i]+"/"+"pypz vs diff", 25, -.60, .60);
+ pxpz_vs_diff[i] = aida.histogram2D(mollerCollections[i]+"/"+"pxpz vs diff", 50, .025, .035, 50, -.60, .60);
+ pypz_vs_diff[i] = aida.histogram2D(mollerCollections[i]+"/"+"pypz vs diff", 50, -.005, .005, 50, -.60, .60);
diff[i] = aida.histogram1D(mollerCollections[i]+"/"+"diff", 50, -.60, .60);
- sum[i] = aida.histogram1D(mollerCollections[i]+"/"+"sum", 50, 1.0, 1.1);
-
- pxpz_vs_sum[i] = aida.profile1D(mollerCollections[i]+"/"+"pxpz vs sum", 25, 1.0, 1.1);
- pypz_vs_sum[i] = aida.profile1D(mollerCollections[i]+"/"+"pypz vs sum", 25, 1.0, 1.1);
-
- pxpz_vs_mass[i] = aida.profile1D(mollerCollections[i]+"/"+"pxpz vs mass", 25, .03, .037);
- pypz_vs_mass[i] = aida.profile1D(mollerCollections[i]+"/"+"pypz vs mass", 25, .03, .037);
-
- //vtx_x[i] = aida.histogram1D(mollerCollections[i]+"/"+"vtx x", 60, -1, 1);
- //vtx_y[i] = aida.histogram1D(mollerCollections[i]+"/"+"vtx y", 60, -1, 1);
- mass[i] = aida.histogram1D(mollerCollections[i]+"/"+"mass", 60, .030, .037);
+ sum[i] = aida.histogram1D(mollerCollections[i]+"/"+"sum", 50, .95, 1.15);
+
+ pxpz_vs_sum[i] = aida.histogram2D(mollerCollections[i]+"/"+"pxpz vs sum", 50, .025, .035, 50, .95, 1.15);
+ pypz_vs_sum[i] = aida.histogram2D(mollerCollections[i]+"/"+"pypz vs sum", 50, -.005, .005, 50, .95, 1.15);
+
+ pxpz_vs_mass[i] = aida.histogram2D(mollerCollections[i]+"/"+"pxpz vs mass", 50, .025, .035, 50,.03, .037);
+ pypz_vs_mass[i] = aida.histogram2D(mollerCollections[i]+"/"+"pypz vs mass", 50, -.005, .005, 50, .03, .037);
+
+ vtx_xy[i] = aida.histogram2D(mollerCollections[i]+"/"+"vtx xy", 50, -5, 5, 50, -5, 5);
+
+ mass[i] = aida.histogram1D(mollerCollections[i]+"/"+"mass", 300, 0, .060);
timediff[i] = aida.histogram1D(mollerCollections[i]+"/"+"time diff", 60, -6, 6);
+
+ trkchi2[i] = aida.histogram1D(mollerCollections[i]+"/"+"trk chi2", 100, 0, 100);
+ vtxchi2[i] = aida.histogram1D(mollerCollections[i]+"/"+"vtx chi2", 100, 0, 100);
+
+ nPassCuts[i] = aida.histogram1D(mollerCollections[i]+"/"+"moller events pass cut per event", 20, 0, 20);
+ meanClusterTime[i] = aida.histogram1D(mollerCollections[i]+"/"+"mean cluster time", 200, 0, 200);
+
+ dphi[i] = aida.histogram1D(mollerCollections[i]+"/"+"dphi (degrees)", 90, -180, 180);
+ theta_theta[i] = aida.histogram1D(mollerCollections[i]+"/"+"sin(theta1_2)sin(theta2_2)", 100, 0, .002);
+ dphi_vs_dtheta[i] = aida.histogram2D(mollerCollections[i]+"/"+"dphi vs dtheta", 180, -180, 180, 100, 0, .05);
}
- /*pypz_tophighE = aida.histogram1D("topHighE pypz", 60, -.005,.005);
- pxpz_tophighE = aida.histogram1D("topHighE pxpz", 60, .025,.035);
- pypz_bottomhighE = aida.histogram1D("bottomHighE pypz", 60, -.005,.005);
- pxpz_bottomhighE = aida.histogram1D("bottomHighE pxpz", 60, .025,.035);*/
- /*if(display){
- IPlotter p = aida.analysisFactory().createPlotterFactory().create();
- StyleUtil.setSize(p, 1300, 900);
- //p.createRegions(3, 2);
- p.createRegions(4, 3);
-
- p.region(0).plot(hpypz);
- p.region(1).plot(hpxpz);
- p.region(2).plot(timediff);
- p.region(3).plot(pypz_vs_diff);
- p.region(4).plot(pxpz_vs_diff);
- p.region(5).plot(diff);
- p.region(6).plot(pypz_vs_sum);
- p.region(7).plot(pxpz_vs_sum);
- p.region(8).plot(sum);
-
- p.region(9).plot(pypz_vs_mass);
- p.region(10).plot(pxpz_vs_mass);
- p.region(11).plot(mass);
- StyleUtil.stylize(p.region(0),"py/pz", "py/pz", "#");
- StyleUtil.stylize(p.region(1),"px/pz", "px/pz", "#");
- StyleUtil.stylize(p.region(2),"time diff (t-b)", "diff (ns)", "#");
- StyleUtil.stylize(p.region(3),"py/pz vs diff", "diff (GeV)", "py/pz");
- StyleUtil.stylize(p.region(4),"px/pz vs diff", "diff (GeV)", "px/pz");
- StyleUtil.stylize(p.region(5),"diff", "diff (GeV)", "#");
-
- StyleUtil.stylize(p.region(6),"py/pz vs sum", "sum (GeV)", "py/pz");
- StyleUtil.stylize(p.region(7),"px/pz vs sum", "sum (GeV)", "px/pz");
- StyleUtil.stylize(p.region(8),"sum", "sum (GeV)", "#");
-
- StyleUtil.stylize(p.region(9),"py/pz vs mass", "mass (GeV)", "py/pz");
- StyleUtil.stylize(p.region(10),"px/pz vs mass", "mass (GeV)", "px/pz");
- StyleUtil.stylize(p.region(11),"mass", "mass (GeV)", "#");
-
- p.show();
-
- IPlotter p2 = aida.analysisFactory().createPlotterFactory().create();
-
- p2.createRegions(2, 1);
-
-
- p2.region(0).plot(hpypz_botHighE);
- p2.region(1).plot(hpxpz_botHighE);
- p2.region(0).plot(hpypz_mid);
- p2.region(1).plot(hpxpz_mid);
- p2.region(0).plot(hpypz_topHighE);
- p2.region(1).plot(hpxpz_topHighE);
-
- StyleUtil.stylize(p2.region(0),"py/pz", "py/pz", "#");
- StyleUtil.stylize(p2.region(1),"px/pz", "py/pz", "#");
- StyleUtil.noFillHistogramBars(p2.region(0));
- StyleUtil.noFillHistogramBars(p2.region(1));
- p2.show();
- }
- */
+
+ }
+
+
+ double beamTiltX = .0295;
+ double beamTiltY = -.0008;
+
+ double[] getThetaAndPhi(Hep3Vector p){
+ double px = p.x(), py = p.y(), pz = p.z();
+
+ double cx = cos(beamTiltX);
+ double sy = sin(beamTiltY);
+ double cy = cos(beamTiltY);
+ double sx = sin(beamTiltX);
+
+ double pxtilt = px*cx -pz*sx;
+ double pytilt = -py*sx*sy + py*cy -pz*sy*cx;
+ double pztilt = px*cy*sx + py*sy +pz*cy*cx;
+
+
+ double theta = atan(hypot(pxtilt, pytilt)/pztilt);
+ double phi = atan2(pytilt, pxtilt);
+ return new double[]{theta, phi};
}
}
Modified: java/trunk/users/src/main/java/org/hps/users/spaul/trident/TridentHistogramDriver.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/trident/TridentHistogramDriver.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/trident/TridentHistogramDriver.java Thu Sep 1 04:26:03 2016
@@ -9,6 +9,7 @@
import org.hps.recon.tracking.TrackType;
import org.hps.users.spaul.StyleUtil;
+import org.lcsim.event.Cluster;
import org.lcsim.event.EventHeader;
import org.lcsim.event.ReconstructedParticle;
import org.lcsim.event.Vertex;
@@ -28,22 +29,25 @@
public void process(EventHeader event){
for(int i = 0; i< tridentCollections.length; i++){
List<Vertex> tridents = event.get(Vertex.class, tridentCollections[i]);
+
for(Vertex v : tridents){
+
+
if(!passesCuts(v))
continue;
- ReconstructedParticle m = v.getAssociatedParticle();
- ReconstructedParticle top;
- ReconstructedParticle bottom;
- if(m.getParticles().get(0).getMomentum().y()>0){
- top = m.getParticles().get(0);
- bottom = m.getParticles().get(1);
+ ReconstructedParticle part = v.getAssociatedParticle();
+ ReconstructedParticle plus;
+ ReconstructedParticle minus;
+ if(part.getParticles().get(0).getCharge()>0){
+ plus = part.getParticles().get(0);
+ minus = part.getParticles().get(1);
}else{
- top = m.getParticles().get(1);
- bottom = m.getParticles().get(0);
+ plus = part.getParticles().get(1);
+ minus = part.getParticles().get(0);
}
- double pypz = m.getMomentum().y()/m.getMomentum().z();
- double pxpz = m.getMomentum().x()/m.getMomentum().z();
+ double pypz = part.getMomentum().y()/part.getMomentum().z();
+ double pxpz = part.getMomentum().x()/part.getMomentum().z();
//double pypz = (top.getMomentum().y()+bottom.getMomentum().y())/(top.getMomentum().z()+bottom.getMomentum().z());
//double pxpz = (top.getMomentum().x()+bottom.getMomentum().x())/(top.getMomentum().z()+bottom.getMomentum().z());
@@ -51,71 +55,104 @@
hpxpz[i].fill(pxpz);
- double diff = top.getMomentum().z()-bottom.getMomentum().z();
- double sum = m.getMomentum().z();//top.getMomentum().z()+bottom.getMomentum().z();
- double mass = m.getMass();
-
-
-
- this.diff[i].fill(diff);
- this.sum[i].fill(sum);
+ double diff = plus.getMomentum().z()-minus.getMomentum().z();
+ double sum = part.getMomentum().z();
+ double mass = part.getMass();
+
+
+
+ this.pz_diff[i].fill(diff);
+ this.pz_sum[i].fill(sum);
this.mass[i].fill(mass);
- pypz_vs_diff[i].fill(diff,pypz );
- pxpz_vs_diff[i].fill(diff, pxpz );
-
-
-
- pxpz_vs_sum[i].fill(sum, pxpz );
- pypz_vs_sum[i].fill(sum, pypz );
-
- pxpz_vs_mass[i].fill(mass, pxpz );
- pypz_vs_mass[i].fill(mass, pypz );
- timediff[i].fill(top.getClusters().get(0).getCalorimeterHits().get(0).getTime()
- -bottom.getClusters().get(0).getCalorimeterHits().get(0).getTime());
- /*if(moreEnergetic.getMomentum().y() > 0)
- {
- pypz_tophighE.fill(pypz);
- pxpz_tophighE.fill(pxpz);
- }
- if(moreEnergetic.getMomentum().y() < 0)
- {
- pypz_bottomhighE.fill(pypz);
- pxpz_bottomhighE.fill(pxpz);
- }*/
+
+ double energy = v.getAssociatedParticle().getEnergy();
+ this.mass_vs_energy[i].fill(mass, energy);
+ this.mass_vs_pz[i].fill(mass, part.getMomentum().z());
+
+ timediff[i].fill(plus.getClusters().get(0).getCalorimeterHits().get(0).getTime()
+ -minus.getClusters().get(0).getCalorimeterHits().get(0).getTime());
+ this.chi2Vtx[i].fill(v.getChi2());
+ this.chi2TrkPlus[i].fill(plus.getTracks().get(0).getChi2());
+ this.chi2TrkMinus[i].fill(minus.getTracks().get(0).getChi2());
+
+ this.vtx_x[i].fill(v.getPosition().x());
+ this.vtx_y[i].fill(v.getPosition().y());
+ this.vtx_xy[i].fill(v.getPosition().x(), v.getPosition().y());
+
+ this.pz_vs_pz[i].fill(minus.getMomentum().z(), plus.getMomentum().z());
}
}
}
- static double BIG_NUMBER = Double.POSITIVE_INFINITY;
- double _maxVtxChi2 = 10;
- double _maxTrkChi2 = 15;
- double _maxMass = 1.0;
+
+ double _maxVtxChi2 = 100;
+ double _maxTrkChi2 = 150;
+ double _maxMass = .2;
double _minMass = 0;
- double _minVtxPz = 2.0;
- double _maxVtxPz = 4.0;
- double _maxTrkPz = 0;
+ double _minVtxPz = 0.0;
+ double _maxVtxPz = 3.0;
+
+ double _maxTrkPz = 1.9;
double _maxPzDiff = 4.0;
+ double _timeDiffCut = 1.2;
boolean passesCuts(Vertex vertex){
ReconstructedParticle part = vertex.getAssociatedParticle();
+
+ ReconstructedParticle p1 = part.getParticles().get(0);
+ ReconstructedParticle p2 = part.getParticles().get(1);
+
+
+ //first make sure the track type is GBL
if(!TrackType.isGBL(part.getType()))
return false;
+ // make sure both tracks are matched to clusters
+ if(p1.getClusters().size() == 0)
+ return false;
+ if(p2.getClusters().size() == 0)
+ return false;
+
+ Cluster c1 = p1.getClusters().get(0);
+ Cluster c2 = p2.getClusters().get(0);
+
+
+ //make sure the clusters are on opposite sides of detector.
+ if(c1.getPosition()[1]*c2.getPosition()[1] >0)
+ return false;
+
+ //plot the time difference (top minus bottom) versus the energy sum
+ double dt = c1.getCalorimeterHits().get(0).getTime()
+ -c2.getCalorimeterHits().get(0).getTime();
+ timediff_vs_esum[0].fill(Math.signum(c1.getPosition()[1])*dt,
+ part.getEnergy());
+
+ //and that they are from the same beam bunch
+
+ if(Math.abs(dt)>_timeDiffCut)
+ return false;
+
+
+ //make sure the total momentum is a reasonable range.
if(part.getMomentum().z() > _maxVtxPz || part.getMomentum().z() < _minVtxPz)
return false;
+
+ // mass within a proper window.
if(part.getMass() > _maxMass || part.getMass() < _minMass)
return false;
-
- if(part.getParticles().get(0).getCharge()*part.getParticles().get(1).getCharge() != -1 )
- return false;
-
+
+ //fee momentum cut
+ if(p1.getMomentum().z() > 1.9 || p2.getMomentum().z() > 1.9)
+ return false;
+
+ //oppositely charged particles.
+ if(part.getParticles().get(0).getCharge()
+ + part.getParticles().get(1).getCharge() != 0)
+ return false;
+
+ // make sure the chi^2 of the vertex fit is reasonable
if(vertex.getChi2() > _maxVtxChi2)
return false;
-
-
- if(part.getParticles().get(0).getClusters().size() == 0)
- return false;
- if(part.getParticles().get(1).getClusters().size() == 0)
- return false;
-
+
+ //and also the chi^2 of the individual tracks are reasonable as well.
if(part.getParticles().get(0).getTracks().get(0).getChi2() > _maxTrkChi2)
return false;
if(part.getParticles().get(1).getTracks().get(0).getChi2() > _maxTrkChi2)
@@ -123,13 +160,20 @@
return true;
}
- IHistogram1D hpypz[], hpxpz[], diff[], sum[], mass[];
+ IHistogram1D hpypz[], hpxpz[], pz_diff[], pz_sum[], mass[];
IHistogram1D chi2TrkPlus[];
IHistogram1D chi2TrkMinus[];
IHistogram1D chi2Vtx[];
-
+ IHistogram1D vtx_x[], vtx_y[], timediff[];
+
+
+ private IHistogram2D[] mass_vs_energy;
+ private IHistogram2D[] vtx_xy;
+ private IHistogram2D[] pz_vs_pz;
+ private IHistogram2D[] mass_vs_pz;
+ private IHistogram2D[] timediff_vs_esum;
@@ -191,10 +235,7 @@
this._maxVtxPz = _maxPz;
}
- IHistogram1D vtx_x[], vtx_y[], timediff[];
-
- IProfile1D pxpz_vs_diff[], pypz_vs_diff[], pxpz_vs_sum[], pypz_vs_sum[],
- pxpz_vs_mass[], pypz_vs_mass[];
+
//IHistogram1D pypz_tophighE, pxpz_tophighE;
@@ -206,51 +247,59 @@
hpxpz = new IHistogram1D[3];
- pxpz_vs_diff= new IProfile1D[3];
- pypz_vs_diff= new IProfile1D[3];
-
- diff= new IHistogram1D[3];
-
- sum= new IHistogram1D[3];
-
- pxpz_vs_sum= new IProfile1D[3];
- pypz_vs_sum= new IProfile1D[3];
-
- pxpz_vs_mass= new IProfile1D[3];
- pypz_vs_mass= new IProfile1D[3];
-
- //vtx_x= new IHistogram1D[3];
- //vtx_y= new IHistogram1D[3];
+ pz_diff= new IHistogram1D[3];
+
+ pz_sum= new IHistogram1D[3];
+
+
+
+
+ vtx_x= new IHistogram1D[3];
+ vtx_y= new IHistogram1D[3];
+
+ vtx_xy = new IHistogram2D[3];
+
mass= new IHistogram1D[3];
timediff= new IHistogram1D[3];
-
+ timediff_vs_esum = new IHistogram2D[3];
+ mass_vs_energy = new IHistogram2D[3];
+ mass_vs_pz = new IHistogram2D[3];
+ chi2Vtx = new IHistogram1D[3];
+ chi2TrkPlus = new IHistogram1D[3];
+ chi2TrkMinus = new IHistogram1D[3];
+ pz_vs_pz = new IHistogram2D[3];
for(int i = 0; i< 3; i++){
- hpypz[i] = aida.histogram1D(tridentCollections[i]+"/"+"pypz", 60, -.005,.005);
- hpxpz[i] = aida.histogram1D(tridentCollections[i]+"/"+"pxpz", 60, .025,.035);
-
-
+ hpypz[i] = aida.histogram1D(tridentCollections[i]+"/"+"pypz", 80, -.020,.020);
+ hpxpz[i] = aida.histogram1D(tridentCollections[i]+"/"+"pxpz", 80, .010,.050);
+
+ vtx_x[i] = aida.histogram1D(tridentCollections[i]+"/"+"vtx x", 200, -2, 2);
+ vtx_y[i] = aida.histogram1D(tridentCollections[i]+"/"+"vtx y", 200, -2, 2);
+
+
+ vtx_xy[i] = aida.histogram2D(tridentCollections[i]+"/"+"vtx xy", 200, -2, 2, 200, -2, 2);
+
+ pz_diff[i] = aida.histogram1D(tridentCollections[i]+"/"+"pz pos - ele",120, -_maxPzDiff, _maxPzDiff);
+
+ pz_sum[i] = aida.histogram1D(tridentCollections[i]+"/"+"pz sum", 100, _minVtxPz, _maxVtxPz);
-
- pxpz_vs_diff[i] = aida.profile1D(tridentCollections[i]+"/"+"pxpz vs diff", 25, -.60, .60);
- pypz_vs_diff[i] = aida.profile1D(tridentCollections[i]+"/"+"pypz vs diff", 25, -.60, .60);
-
- diff[i] = aida.histogram1D(tridentCollections[i]+"/"+"diff", 50, -.60, .60);
-
- sum[i] = aida.histogram1D(tridentCollections[i]+"/"+"sum", 50, 1.0, 1.1);
-
- pxpz_vs_sum[i] = aida.profile1D(tridentCollections[i]+"/"+"pxpz vs sum", 25, 1.0, 1.1);
- pypz_vs_sum[i] = aida.profile1D(tridentCollections[i]+"/"+"pypz vs sum", 25, 1.0, 1.1);
-
- pxpz_vs_mass[i] = aida.profile1D(tridentCollections[i]+"/"+"pxpz vs mass", 25, .03, .037);
- pypz_vs_mass[i] = aida.profile1D(tridentCollections[i]+"/"+"pypz vs mass", 25, .03, .037);
-
//vtx_x[i] = aida.histogram1D(tridentCollections[i]+"/"+"vtx x", 60, -1, 1);
//vtx_y[i] = aida.histogram1D(tridentCollections[i]+"/"+"vtx y", 60, -1, 1);
- mass[i] = aida.histogram1D(tridentCollections[i]+"/"+"mass", 60, .030, .037);
+
+ //for 0-200 MeV, this yields 0.1 MeV mass bins
+ mass[i] = aida.histogram1D(tridentCollections[i]+"/"+"mass", 2000, _minMass, _maxMass);
+ //
+ mass_vs_energy[i] = aida.histogram2D(tridentCollections[i]+"/"+"mass vs energy sum", 2000, _minMass, _maxMass, 100, 0, 5);
+ mass_vs_pz[i] = aida.histogram2D(tridentCollections[i]+"/"+"mass vs pz sum", 2000, _minMass, _maxMass, 100, 0, 5);
+
+
timediff[i] = aida.histogram1D(tridentCollections[i]+"/"+"time diff", 60, -6, 6);
+ timediff_vs_esum[i] = aida.histogram2D(tridentCollections[i]+"/"+"time diff vs energy sum", 60, -6, 6, 100, 0, 3);
+
+ chi2Vtx[i] = aida.histogram1D(tridentCollections[i]+"/"+"chi2 vertex", 100, 0, 100);
+ chi2TrkPlus[i] = aida.histogram1D(tridentCollections[i]+"/"+"chi2 track plus", 100, 0, 100);
+ chi2TrkMinus[i] = aida.histogram1D(tridentCollections[i]+"/"+"chi2 track minus", 100, 0, 100);
+ pz_vs_pz[i] = aida.histogram2D(tridentCollections[i]+"/"+"pz ele vs pos", 100, 0, 4, 100, 0, 4);
}
-
-
}
}
