Author: [log in to unmask]
Date: Mon May 16 12:45:06 2016
New Revision: 4353
Log: (empty)
Added:
java/trunk/users/src/main/java/org/hps/users/spaul/fee/
java/trunk/users/src/main/java/org/hps/users/spaul/fee/Cleanup.java
java/trunk/users/src/main/java/org/hps/users/spaul/fee/CustomBinning.java
java/trunk/users/src/main/java/org/hps/users/spaul/fee/FeeHistogramDriver.java
java/trunk/users/src/main/java/org/hps/users/spaul/fee/FormFactor.java
java/trunk/users/src/main/java/org/hps/users/spaul/fee/MottIntegral.java
java/trunk/users/src/main/java/org/hps/users/spaul/trident/
java/trunk/users/src/main/java/org/hps/users/spaul/trident/TridentHistogramDriver.java
Modified:
java/trunk/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java
java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java
java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinningXY.java
Added: 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 (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/fee/Cleanup.java Mon May 16 12:45:06 2016
@@ -0,0 +1,53 @@
+package org.hps.users.spaul.fee;
+
+import hep.physics.vec.Hep3Vector;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.lcsim.event.ReconstructedParticle;
+
+public class Cleanup {
+
+ public static void purgeDuplicates(List<ReconstructedParticle> particles) {
+ ArrayList<ReconstructedParticle> trashcan = new ArrayList();
+ outer : for(ReconstructedParticle p1 : particles){
+ if(trashcan.contains(p1))
+ continue;
+ for(ReconstructedParticle p2 : particles){
+ if(p1 == p2) continue;
+ if(trashcan.contains(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
+ ReconstructedParticle loser = getWorseParticle(p1, p2);
+ trashcan.add(loser);
+ }
+ if(trashcan.contains(p1))
+ continue outer;
+
+ }
+
+ }
+ particles.removeAll(trashcan);
+ }
+
+ static ReconstructedParticle getWorseParticle(ReconstructedParticle p1, ReconstructedParticle p2){
+ 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());
+ }
+
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/fee/CustomBinning.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/fee/CustomBinning.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/fee/CustomBinning.java Mon May 16 12:45:06 2016
@@ -0,0 +1,100 @@
+package org.hps.users.spaul.fee;
+
+import java.io.File;
+import java.io.FileNotFoundException;
+import java.util.Scanner;
+
+public class CustomBinning {
+
+ public CustomBinning(String s){
+ this(new Scanner(s));
+ }
+
+
+ public CustomBinning(File f) throws FileNotFoundException{
+ this(new Scanner(f));
+ }
+
+ private CustomBinning(Scanner s){
+
+
+
+ nTheta = s.nextInt(); //number of bins in theta;
+ thetaMax = new double[nTheta];
+ thetaMin = new double[nTheta];
+
+ phiMax = new double[nTheta][];
+ phiMin = new double[nTheta][];
+ int i = 0;
+ while(s.hasNext()){ //new row
+ int nPhi = s.nextInt();
+ thetaMin[i] = s.nextDouble();
+ thetaMax[i] = s.nextDouble();
+
+ phiMax[i] = new double[nPhi];
+ phiMin[i] = new double[nPhi];
+ for(int j = 0; j<nPhi; j++){
+ phiMin[i][j] = s.nextDouble();
+ phiMax[i][j] = s.nextDouble();
+ }
+ i++;
+ }
+ }
+ double[][] phiMax;
+ double[][] phiMin;
+ public double thetaMax[], thetaMin[];
+ public int nTheta;
+
+ double getSteradians(int binNumber){
+ double t1 = thetaMin[binNumber];
+ double t2 = thetaMax[binNumber];
+ double dCos = Math.cos(t1)-Math.cos(t2);
+ double dPhiTot = 0;
+ for(int i = 0; i< phiMax[binNumber].length; i++){
+ dPhiTot += phiMax[binNumber][i]-phiMin[binNumber][i];
+ }
+ return 2*dPhiTot*dCos; //factor of two because top and bottom
+ }
+
+ double getPhiWidth(int thetaBin){
+ double dPhiTot = 0;
+ for(int i = 0; i< phiMax[thetaBin].length; i++){
+ dPhiTot += phiMax[thetaBin][i]-phiMin[thetaBin][i];
+ }
+ return 2*dPhiTot; //factor of two because top and bottom
+
+ }
+ boolean inRange(double theta, double phi){
+ phi = Math.abs(phi);
+ /*int i =(int) Math.floor((theta-theta0)/deltaTheta);
+ if(i>= nTheta || i<0)
+ return false;*/
+ if(theta > thetaMax[nTheta-1] || theta < thetaMin[0])
+ return false;
+ int i;
+ boolean found = false;
+ for(i = 0; i< nTheta; i++){
+ if(theta > thetaMin[i] && theta < thetaMax[i]){
+ found = true;
+ break;
+ }
+ }
+ if(!found)
+ return false;
+
+ for(int j = 0; j<phiMax[i].length; j++){
+ if(phi>phiMin[i][j] && phi< phiMax[i][j])
+ return true;
+ }
+ return false;
+
+ }
+ public double getTotSteradians() {
+ double tot = 0;
+ for(int i = 0; i<nTheta; i++){
+ tot += getSteradians(i);
+ }
+ return tot;
+ }
+
+}
Added: 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 (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/fee/FeeHistogramDriver.java Mon May 16 12:45:06 2016
@@ -0,0 +1,386 @@
+package org.hps.users.spaul.fee;
+
+import java.util.List;
+
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.physics.vec.Hep3Vector;
+
+import org.hps.conditions.beam.BeamEnergy.BeamEnergyCollection;
+import org.hps.recon.tracking.TrackType;
+import org.hps.recon.tracking.TrackUtils;
+import org.hps.record.triggerbank.AbstractIntData;
+import org.hps.record.triggerbank.TIData;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.GenericObject;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.geometry.Detector;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+import static java.lang.Math.*;
+
+
+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 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;
+
+
+ public double getMinClustTime() {
+ return minClustTime;
+ }
+
+
+
+ public void setMinClustTime(double minClustTime) {
+ this.minClustTime = minClustTime;
+ }
+
+
+
+ public double getMaxClustTime() {
+ return maxClustTime;
+ }
+
+
+
+ public void setMaxClustTime(double maxClustTime) {
+ this.maxClustTime = maxClustTime;
+ }
+
+ IHistogram2D z0VsTanLambda;
+ IHistogram2D z0VsChi2;
+
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/fee/FormFactor.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/fee/FormFactor.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/fee/FormFactor.java Mon May 16 12:45:06 2016
@@ -0,0 +1,43 @@
+package org.hps.users.spaul.fee;
+
+public abstract class FormFactor {
+
+ abstract double getFormFactor(double q2);
+
+ double getFormFactorSquared(double q2){
+ return Math.pow(getFormFactor(q2), 2);
+ }
+ static double hbarc = .197;
+ static FormFactor carbon = new FormFactor(){
+
+ @Override
+ public double getFormFactor(double Q2) {
+ double Z= 6;
+ double rp = 0.8786;
+ double a = 1.64;
+ double b = Math.sqrt(a*a*(1-1/12.)+rp*rp);
+ return (1-(Z-2)/(6*Z)*a*a*Q2/(hbarc*hbarc))*Math.exp(-1/4.*b*b*Q2/(hbarc*hbarc));
+ }
+
+ };
+ static FormFactor tungsten = new FormFactor(){
+
+ @Override
+ public double getFormFactor(double Q2) {
+ double Z= 74;
+ double r = 6.87;
+ double x = Math.sqrt(Q2*r*r/(hbarc*hbarc));
+ double F_calc = 3/(x*x*x)*(Math.sin(x)-x*Math.cos(x));
+ return F_calc;
+ }
+
+ };
+ static FormFactor get(int Z){
+ if(Z == 6)
+ return carbon;
+ if(Z == 74)
+ return tungsten;
+ System.err.println("Form Factor for " + Z + " not implemented");
+ return null;
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/fee/MottIntegral.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/fee/MottIntegral.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/fee/MottIntegral.java Mon May 16 12:45:06 2016
@@ -0,0 +1,66 @@
+package org.hps.users.spaul.fee;
+
+public class MottIntegral {
+ /**
+ * @param bin
+ * @param scale = Z^2alpha^2/(4*E^2)
+ * @param recoil = 2E/M
+ * @return the integral of 1/sin^4(th/2)*cos^2(th/2)/(1+a*sin^2(th/2)) times dPhi*sin(theta),
+ * which appears in the integral of mott scattering.
+ *
+ * NOTE the sin(theta) is because dOmega = dTheta*dPhi*sin(theta)
+ */
+ static double mottIntegral(double recoil, double scale, int bin, CustomBinning cb){
+ double dPhi = 0;
+ for(int i = 0; i< cb.phiMax[bin].length; i++)
+ dPhi += 2*(cb.phiMax[bin][i] - cb.phiMin[bin][i]); //factor of 2 from top and bottom
+
+ double Imax = integral(recoil, cb.thetaMax[bin]);
+
+ double Imin = integral(recoil, cb.thetaMin[bin]);
+ double dI = Imax-Imin;
+
+
+ double retval = scale*dPhi*dI;
+ return retval;
+ }
+
+ static double mottIntegral(double recoil, double scale, double thetaMin, double thetaMax){
+ double Imax = integral(recoil, thetaMax);
+
+ double Imin = integral(recoil, thetaMin);
+ double dI = Imax-Imin;
+
+ double dPhi = 2*Math.PI; //full range in phi
+ double retval = scale*dPhi*dI;
+ return retval;
+ }
+
+ static double integral(double a, double th){
+ double sinth22 = Math.pow(Math.sin(th/2), 2);
+ return 2*(-1/sinth22+(1+a)*(Math.log(2/sinth22+2*a)));
+ }
+
+ static double mottIntegralWithFormFactor(double recoil, double scale, double thetaMin, double thetaMax, FormFactor ff, int nsteps, double E){
+ double sum = 0;
+ double prev = -1;
+
+ double dTheta = (thetaMax-thetaMin)/nsteps;
+ for(int i = 0; i<nsteps; i++){
+ double theta = i*(thetaMax-thetaMin)/nsteps+thetaMin;
+ double I = integral(recoil, theta);
+ if(i != 0)
+ {
+ double ts = Math.pow(Math.sin(theta/2),2);
+ double q2 = 4*E*E*ts/(1+recoil*ts);
+ double f2 = ff.getFormFactorSquared(q2);
+ sum+= (I-prev)*f2;
+ }
+ prev = I;
+
+ }
+ double dPhi = 2*Math.PI; //full range in phi
+ double retval = scale*dPhi*sum;
+ return retval;
+ }
+}
Modified: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java Mon May 16 12:45:06 2016
@@ -21,7 +21,7 @@
}
pw.println();
}
- ShowCustomBinning.main(new String[]{"generatedbins.txt"});
+ ShowCustomBinning.main(new String[]{"generatedbins.txt", "."});
}
@@ -60,10 +60,11 @@
// make the angular cuts on the tracks such that the particles that go into that cut
// are expected to be within 4 mm (~= 2 times the angular resolution of 1.5 mrad) of
// the ecal cuts.
- double d = 4;
+ double d = 0;
boolean inRange = EcalUtil.fid_ECal_spherical_more_strict(thetaMin, phi, d) && EcalUtil.fid_ECal_spherical_more_strict(thetaMax, phi, d)
- && EcalUtil.fid_ECal_spherical_more_strict(thetaMin, -phi, d) && EcalUtil.fid_ECal_spherical_more_strict(thetaMax, -phi, d);
+ && EcalUtil.fid_ECal_spherical_more_strict(thetaMin, -phi, d) && EcalUtil.fid_ECal_spherical_more_strict(thetaMax, -phi, d)
+ && svt_fiducial(thetaMin, phi)&& svt_fiducial(thetaMin, -phi)&& svt_fiducial(thetaMax, phi) && svt_fiducial(thetaMax, -phi);
if(inRange && !prevInRange)
phiBins[edgeNumber++] = phi;
if(prevInRange && !inRange)
@@ -87,4 +88,15 @@
}
+ private static boolean svt_fiducial(double theta, double phi) {
+ double ux = Math.sin(theta)*Math.cos(phi);
+ double uy = Math.sin(theta)*Math.sin(phi);
+ double uz = Math.cos(theta);
+
+ ux/= uz;
+ uy/= uz;
+
+ return Math.abs(uy)> 0.008 && Math.abs(uy)<0.059 && ux < 0.1600 && ux > -.119;
+ }
+
}
Modified: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java Mon May 16 12:45:06 2016
@@ -70,13 +70,63 @@
drawEcalOutline(g);
drawFidEcalOutline(g);
+ drawSVTOutline5(g);
drawCustomBinRectangles(g);
g.setColor(Color.BLACK);
drawXAxis(g);
drawYAxis(g);
}
- void drawFidEcalOutline(Graphics g){
+ private void drawSVTOutline5(Graphics g) {
+ g.setColor(Color.RED);
+ double x_edge_high = 0.160;
+ double y_edge_low = 0.008;
+ double y_edge_high = 0.059;
+ double x_edge_low = -.119;
+
+ double ux1,uy1, ux2, uy2;
+ double nPoints = 200;
+ for(int i = 0; i< nPoints-1; i++){
+ ux1 = x_edge_high;
+ ux2 = x_edge_high;
+ uy1 = y_edge_low + i*(y_edge_high-y_edge_low)/nPoints;
+ uy2 = y_edge_low+ (i+1)*(y_edge_high-y_edge_low)/nPoints;
+ drawLine(g, ux1, uy1, ux2, uy2);
+ drawLine(g, ux1, -uy1, ux2, -uy2);
+
+ ux1 = x_edge_low;
+ ux2 = x_edge_low;
+ uy1 = y_edge_low + i*(y_edge_high-y_edge_low)/nPoints;
+ uy2 = y_edge_low+ (i+1)*(y_edge_high-y_edge_low)/nPoints;
+ drawLine(g, ux1, uy1, ux2, uy2);
+ drawLine(g, ux1, -uy1, ux2, -uy2);
+
+ ux1 = x_edge_low + i*(x_edge_high-x_edge_low)/nPoints;
+ ux2 = x_edge_low + (i+1)*(x_edge_high-x_edge_low)/nPoints;
+ uy1 = y_edge_low;
+ uy2 = y_edge_low;
+ drawLine(g, ux1, uy1, ux2, uy2);
+ drawLine(g, ux1, -uy1, ux2, -uy2);
+
+ ux1 = x_edge_low + i*(x_edge_high-x_edge_low)/nPoints;
+ ux2 = x_edge_low + (i+1)*(x_edge_high-x_edge_low)/nPoints;
+ uy1 = y_edge_high;
+ uy2 = y_edge_high;
+ drawLine(g, ux1, uy1, ux2, uy2);
+ drawLine(g, ux1, -uy1, ux2, -uy2);
+ }
+ }
+
+ protected void drawLine(Graphics g, double ux1, double uy1, double ux2, double uy2){
+ double theta1 = Math.atan(Math.hypot(ux1, uy1));
+ double theta2 = Math.atan(Math.hypot(ux2, uy2));
+
+ double phi1 = Math.atan2(uy1, ux1);
+ double phi2 = Math.atan2(uy2, ux2);
+
+ g.drawLine(getX(theta1), getY(phi1), getX(theta2), getY(phi2));
+ }
+ void drawFidEcalOutline(Graphics g){
g.setColor(Color.GRAY);
double x_edge_low = -262.74;
double x_edge_high = 347.7;
Modified: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinningXY.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinningXY.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinningXY.java Mon May 16 12:45:06 2016
@@ -45,6 +45,15 @@
}
}
}
+
+ protected void drawLine(Graphics g, double ux1, double uy1, double ux2, double uy2){
+
+ ux1+=.0305;
+ ux2+=.0305;
+
+ g.drawLine(getX(ux1), getY(uy1), getX(ux2), getY(uy2));
+ }
+
private void drawLineFromPolar(double theta1, double phi1, double theta2,
double phi2, Graphics g, Polygon p) {
double[] xy1 = toXY(theta1, phi1);
Added: 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 (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/trident/TridentHistogramDriver.java Mon May 16 12:45:06 2016
@@ -0,0 +1,256 @@
+package org.hps.users.spaul.trident;
+
+import java.util.List;
+
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.aida.IPlotter;
+import hep.aida.IProfile1D;
+
+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;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+public class TridentHistogramDriver extends Driver{
+
+ private String[] tridentCollections = new String[]{
+ "TargetConstrainedV0Vertices",
+ "UnconstrainedV0Vertices",
+ "BeamspotConstrainedV0Vertices",
+ };
+
+
+ @Override
+ 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);
+ }else{
+ top = m.getParticles().get(1);
+ bottom = m.getParticles().get(0);
+ }
+
+ double pypz = m.getMomentum().y()/m.getMomentum().z();
+ double pxpz = m.getMomentum().x()/m.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());
+
+ hpypz[i].fill(pypz);
+ 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);
+ 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);
+ }*/
+ }
+ }
+ }
+
+ static double BIG_NUMBER = Double.POSITIVE_INFINITY;
+ double _maxVtxChi2 = 10;
+ double _maxTrkChi2 = 15;
+ double _maxMass = 1.0;
+ double _minMass = 0;
+ double _minVtxPz = 2.0;
+ double _maxVtxPz = 4.0;
+ double _maxTrkPz = 0;
+ double _maxPzDiff = 4.0;
+ boolean passesCuts(Vertex vertex){
+ ReconstructedParticle part = vertex.getAssociatedParticle();
+ if(!TrackType.isGBL(part.getType()))
+ return false;
+ if(part.getMomentum().z() > _maxVtxPz || part.getMomentum().z() < _minVtxPz)
+ return false;
+ if(part.getMass() > _maxMass || part.getMass() < _minMass)
+ return false;
+
+ if(part.getParticles().get(0).getCharge()*part.getParticles().get(1).getCharge() != -1 )
+ return false;
+
+ 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;
+
+ if(part.getParticles().get(0).getTracks().get(0).getChi2() > _maxTrkChi2)
+ return false;
+ if(part.getParticles().get(1).getTracks().get(0).getChi2() > _maxTrkChi2)
+ return false;
+ return true;
+ }
+
+ IHistogram1D hpypz[], hpxpz[], diff[], sum[], mass[];
+
+ IHistogram1D chi2TrkPlus[];
+ IHistogram1D chi2TrkMinus[];
+ IHistogram1D chi2Vtx[];
+
+
+
+
+
+ public double getMaxVtxChi2() {
+ return _maxVtxChi2;
+ }
+
+
+ public void setMaxVtxChi2(double _maxVtxChi2) {
+ this._maxVtxChi2 = _maxVtxChi2;
+ }
+
+
+ public double getMaxTrkChi2() {
+ return _maxTrkChi2;
+ }
+
+
+ public void setMaxTrkChi2(double _maxTrkChi2) {
+ this._maxTrkChi2 = _maxTrkChi2;
+ }
+
+
+ public double getMaxMass() {
+ return _maxMass;
+ }
+
+
+ public void setMaxMass(double _maxMass) {
+ this._maxMass = _maxMass;
+ }
+
+
+ public double getMinMass() {
+ return _minMass;
+ }
+
+
+ public void setMinMass(double _minMass) {
+ this._minMass = _minMass;
+ }
+
+
+ public double getMinVtxPz() {
+ return _minVtxPz;
+ }
+
+
+ public void setMinVtxPz(double _minPz) {
+ this._minVtxPz = _minPz;
+ }
+
+
+ public double getMaxVtxPz() {
+ return _maxVtxPz;
+ }
+
+ public void setMaxVtxPz(double _maxPz) {
+ 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;
+ //IHistogram1D pypz_bottomhighE, pxpz_bottomhighE;
+ @Override
+ public void startOfData(){
+ AIDA aida = AIDA.defaultInstance();
+ hpypz = new IHistogram1D[3];
+ 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];
+ mass= new IHistogram1D[3];
+ timediff= new IHistogram1D[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);
+
+
+
+
+ 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);
+ timediff[i] = aida.histogram1D(tridentCollections[i]+"/"+"time diff", 60, -6, 6);
+ }
+
+
+ }
+}
