Commit in hps-java/src/main/java/org/lcsim/hps/util on MAIN | |||
CalculateAcceptanceFromMadGraph.java | +669 | added 1.1 |
main class for getting acceptance from radiative MadGraph events
diff -N CalculateAcceptanceFromMadGraph.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ CalculateAcceptanceFromMadGraph.java 5 Sep 2012 00:17:54 -0000 1.1 @@ -0,0 +1,669 @@
+/* + * To change this template, choose Tools | Templates + * and open the template in the editor. + */ +package org.lcsim.hps.util; + +/** + + @author mgraham + */ +import java.io.FileReader; +import java.io.IOException; +import java.io.StreamTokenizer; +import java.util.ArrayList; +import java.util.List; + +import hep.io.stdhep.StdhepBeginRun; +import hep.io.stdhep.StdhepEndRun; +import hep.io.stdhep.StdhepEvent; +import hep.io.stdhep.StdhepWriter; +import hep.physics.vec.BasicHep3Vector; +import hep.physics.vec.Hep3Vector; +import java.io.*; +import java.util.*; +import java.util.logging.Level; +import java.util.logging.Logger; +import org.apache.commons.cli.CommandLine; +import org.apache.commons.cli.CommandLineParser; +import org.apache.commons.cli.HelpFormatter; +import org.apache.commons.cli.Option; +import org.apache.commons.cli.Options; +import org.apache.commons.cli.ParseException; +import org.apache.commons.cli.PosixParser; +import org.lcsim.detector.IRotation3D; +import org.lcsim.detector.RotationGeant; +import org.lcsim.fit.helicaltrack.HelixParamCalculator; + +public class CalculateAcceptanceFromMadGraph { + + static boolean expDecay = false; //set a ctau decay length + static boolean flatDecay = false; //decay uniformily in some range + static boolean trident = false; //are these trident events or A' signal events + static double _declength = 0.0; //A' decay length (mm) + static double _xoff = 0.0; //set the x,y,z origin offsets... + static double _yoff = 0.0; + static double _zoff = 0.03; + static double aMass = 99; //Aprime mass (MeV) + static double sigx = 0.00001; // Luminous region size in x (mm) + static double sigy = 0.02; // Luminous region size in y (mm) + static double sigz = 0.02; // Luminous region size in z (mm) + static double maxLen = 200; // maximum decay length (mm) + static double _ecm = 2.2; //GeV + static int nInFiles = 1; //number of input files + static int nmax = 500000; //maximum number of events to write to 1 stdhep file (new one opens if n>nmax) + static String fileType = "lhe"; + static int nread = -1; // a running total of number of events read/written to stdhep files + static boolean _eventFilter = false; + static int _nEleRequired = 2; + static double bField = 0.5; + static double _pCut=0.1; + static double sensorWidth = 40; + static double sensorLength = 100; + static double gap = Math.sin(0.015); + static int nLayers = 7; + static double[] x = {50, 100, 200, 300, 500, 700, 900}; + static double[] ySize = {sensorLength / 2, sensorLength / 2, sensorLength / 2, sensorLength / 2, sensorLength / 2, sensorLength / 2, sensorLength / 2}; + static double[] zSize = {sensorWidth, sensorWidth, sensorWidth, sensorWidth, sensorWidth, sensorWidth, sensorWidth}; + static double[] zGap = {x[1] * gap, x[1] * gap, x[2] * gap, x[3] * gap, x[4] * gap, x[5] * gap, x[6] * gap}; + static int nGenerated = 0; + static int nPassTest = 0; + static int nPassTestLyr50 = 0; + static boolean elePass; + static boolean posPass; + static boolean recoilPass; + static int nbins = 160; + static double maxMass = 400.0; + static int[] genMass = new int[nbins]; + static int[] recoMassLyr1 = new int[nbins]; + static int[] recoMassLyr50 = new int[nbins]; + + public static void doAccounting(List<Boolean> passEle, List<Boolean> passPos, double invMass) { + boolean passLyr700 = false; + boolean passLyr100 = false; + boolean passLyr500 = false; + boolean passLyr50 = false; + //find the bin... + int bin = (int) (invMass / maxMass * nbins); +// System.out.println("invMass = "+invMass + "...goes in bin #"+bin); + genMass[bin]++; + //see if pos+ele is in layer 70cm layer (required form both test & full) + if (passEle.get(5) && passPos.get(5)) + passLyr700 = true; + //see if they were in layer 1 + if (passEle.get(1) && passPos.get(1)) + passLyr100 = true; + // for an event in layer 50, maybe only need for it to hit in layer @ 500 + if (passEle.get(4) && passPos.get(4)) + passLyr500 = true; + if (passEle.get(0) && passPos.get(0)) + passLyr50 = true; + + //ok...fill the histograms + if (passLyr700 && passLyr100) + recoMassLyr1[bin]++; + + if (passLyr700 && passLyr50) //require they go through 6 layers... + recoMassLyr50[bin]++; + + } + + private static Options createCommandLineOptions() { + Options options = new Options(); + + options.addOption(new Option("e", true, "Beam Energy (GeV)")); + options.addOption(new Option("n", true, "Number of files to run.")); + + return options; + } + + /** + @param args the command line arguments + @throws IOException + */ + public static void main(String[] args) throws IOException { + + + // Set up command line parsing. + Options options = createCommandLineOptions(); + + CommandLineParser parser = new PosixParser(); + + // Parse command line arguments. + CommandLine cl = null; + try { + cl = parser.parse(options, args); + System.out.println("Trying parser"); + } catch (ParseException e) { + throw new RuntimeException("Problem parsing command line options.", e); + } + + String ninString = String.valueOf(nInFiles); + String ecmString = String.valueOf(_ecm); + String eptString = String.valueOf(_ecm); + eptString = convertDecimal(eptString); + + + if (cl.hasOption("n")) { + ninString = cl.getOptionValue("n"); + nInFiles = Integer.valueOf(ninString); + System.out.println(ninString); + } + + if (cl.hasOption("e")) { + ecmString = cl.getOptionValue("e"); + _ecm = Double.valueOf(ninString); + System.out.println(ecmString); + eptString = convertDecimal(ecmString); + } + + + + + String fDir = "/nfs/slac/g/hps/mgraham/DarkPhoton/MadGraph/Events" + eptString + "Rad/"; +// String fDir = "./"; +// String inLabel = "W" + eptString + "GeV_Ap" + massString + "MeV_"; + String inLabel = "W" + eptString + "GeV_Rad_"; + + String inPost = "_unweighted_events.lhe"; + + +// String outDir = "/nfs/slac/g/hps/mgraham/DarkPhoton/SignalEvents/"; + String outDir = "./"; +// for (int i = 0; i < nInFiles; i++) { + for (int i = 933; i < 933 + nInFiles; i++) { + int fnum = i + 1; + String snum = "_" + fnum; + if (fnum < 10000) { + snum = "_0" + fnum; + } + if (fnum < 1000) { + snum = "_00" + fnum; + } + if (fnum < 100) { + snum = "_000" + fnum; + } + if (fnum < 10) { + snum = "_0000" + fnum; + } + + String infile = ""; + if (fileType.contains("dat")) { +// infile = fDir + fileLabel + snum + ".dat"; +// System.out.println(infile); + } else if (fileType.contains("lhe")) { + infile = fDir + inLabel + i + inPost; + System.out.println("Unzipping " + infile); + String cmd = "gunzip " + infile + ".gz"; + Process p = Runtime.getRuntime().exec(cmd); + try { + p.waitFor(); + } catch (InterruptedException ex) { + Logger.getLogger(CalculateAcceptanceFromMadGraph.class.getName()).log(Level.SEVERE, null, ex); + } + } + File f = new File(infile); + if (nread == -1 && f.exists()) { + System.out.println("==== processing " + infile + " ===="); + } + if (f.exists()) { + nread += process(infile); + } + if (fileType.contains("lhe")) { + String cmd = "gzip " + infile; + Process p = Runtime.getRuntime().exec(cmd); + try { + p.waitFor(); + } catch (InterruptedException ex) { + Logger.getLogger(CalculateAcceptanceFromMadGraph.class.getName()).log(Level.SEVERE, null, ex); + } + } + + } + //ok ... now spit out histograms + System.out.println("********* Histograms *********** "); + System.out.println("bin mass Gen Lyr1 Lyr50"); + for (int kk = 0; kk < nbins; kk++) { + double mass = ((double) kk / nbins) * maxMass; +// System.out.println(kk+"\t"+mass+"\t"+genMass[kk]+"\t"+recoMassLyr1[kk]+"\t"+recoMassLyr50[kk]); + System.out.printf("%d\t%4.1f\t%d\t%d\t%d\n", kk, mass, genMass[kk], recoMassLyr1[kk], recoMassLyr50[kk]); + } + + outputFile(outDir + "rad" + eptString + "_lyr1.dat", recoMassLyr1); + outputFile(outDir + "rad" + eptString + "_lyr50.dat", recoMassLyr50); + + } + + @SuppressWarnings("static-access") + private static int lineCounter(StreamTokenizer tok) throws IOException { + int lines = 0; + while (tok.nextToken() != tok.TT_EOF) { + if (tok.ttype == tok.TT_EOL) { + lines++; + } + if (tok.ttype == tok.TT_WORD && tok.sval.startsWith("nev")) { + return lines; + } + } + //shouldn't get here...but maybe + return lines; + } + + private static void outputFile(String fName, int[] hist) throws IOException { + FileOutputStream fos = new FileOutputStream(fName); + PrintWriter osw = new PrintWriter(fos); + for (int kk = 0; kk < nbins; kk++) { +// double mass = ((double)kk/nbins)*maxMass; +// System.out.println(kk+"\t"+mass+"\t"+genMass[kk]+"\t"+recoMassLyr1[kk]+"\t"+recoMassLyr50[kk]); +// System.out.printf("%d\t%4.1f\t%d\t%d\t%d\n",kk,mass,genMass[kk],recoMassLyr1[kk],recoMassLyr50[kk]); + osw.println(hist[kk]); + } + osw.close(); + fos.close(); + + } + + private static int getnevts(StreamTokenizer lctok) throws IOException { + int nevts = -1; + if (fileType.contains("dat")) { + return lineCounter(lctok); + } else if (fileType.contains("lhe")) { + while (nevts == -1) { + nevts = getNumberOfEvents(lctok); + } + return nevts; + } + return nevts; + } + + private static int process(String infile) throws IOException { + Random generator = new Random(); + + FileReader lc = new FileReader(infile); + StreamTokenizer lctok = new StreamTokenizer(lc); + lctok.resetSyntax(); + lctok.wordChars(33, 255); + lctok.parseNumbers(); + + lctok.whitespaceChars(0, ' '); + lctok.eolIsSignificant(true); + int nevts = getnevts(lctok); + lc.close(); + System.out.println("Found " + nevts + " events"); + + FileReader fr = new FileReader(infile); + + StreamTokenizer tok = new StreamTokenizer(fr); + + tok.resetSyntax(); + tok.wordChars(33, 255); + tok.parseNumbers(); + + tok.whitespaceChars(0, ' '); + tok.eolIsSignificant(true); + +// System.out.println("Found " + nevts + " events"); + int nreq = (int) nevts; + int ngen = (int) nevts; + int nwrit = (int) nevts; + float ecm = (float) _ecm; + float xsec = (float) 99999997952.; + double rn1 = 12345321; + double rn2 = 66666666; +// StdhepBeginRun sb = new StdhepBeginRun(nreq, ngen, nwrit, ecm, xsec, rn1, rn2); +// sw.writeRecord(sb); + + + tok.resetSyntax(); + tok.wordChars(33, 255); + tok.wordChars('0', '9'); // java.io.StreamTokenizer fails to parse + tok.wordChars('e', 'e'); // scientific notation like "1.09E-008". + tok.wordChars('E', 'E'); // The solution is to read and parse + tok.wordChars('.', '.'); // coordinates as "words". + tok.wordChars('+', '+'); // You run into trouble if the input file + tok.wordChars('-', '-'); // contains text with "e" or "E" which is + tok.whitespaceChars(0, ' '); + tok.eolIsSignificant(true); + + + + double[] beam = {0, 0, 0, 0}; + int nevhep = 0; + for (int icross = 0; icross < nwrit; icross++) { + Hep3Vector beamVec = + new BasicHep3Vector(sigx * generator.nextGaussian() + _xoff, + sigy * generator.nextGaussian() + _yoff, + sigz * generator.nextGaussian() + _zoff); + + beam[0] = beamVec.x(); + beam[1] = beamVec.y(); + beam[2] = beamVec.z(); + + double tmpDecLen = 0; + + + readLHEEvent(tok, beam, icross); + + + } + fr.close(); + + + return nwrit; + + + } + + @SuppressWarnings("static-access") + private static List<Double> getNumbersInLine(StreamTokenizer tok) throws IOException { + List<Double> nums = new ArrayList<Double>(); + while (tok.nextToken() != tok.TT_EOF) { + if (tok.ttype == tok.TT_EOL) { + break; + } + String tokVal = tok.sval; +// System.out.println(tokVal); + nums.add(Double.valueOf(tokVal).doubleValue()); + } + + return nums; + } + + @SuppressWarnings("static-access") + private static int getNumberOfEvents(StreamTokenizer tok) throws IOException { + boolean fndNumber = false; + boolean fndOf = false; + boolean fndEvents = false; + int evts = -1; + while (tok.nextToken() != tok.TT_EOF) { +// System.out.println(tok.toString()); + if (tok.ttype == tok.TT_EOL) { + break; + } + if (tok.ttype == tok.TT_WORD && tok.sval.contentEquals("Number")) { + // System.out.println(tok.toString()); + fndNumber = true; + } + if (tok.ttype == tok.TT_WORD && tok.sval.contentEquals("of")) { + fndOf = true; + } + if (tok.ttype == tok.TT_WORD && tok.sval.contentEquals("Events")) { + fndEvents = true; + } + if (tok.ttype == tok.TT_NUMBER && fndEvents && fndOf && fndNumber) { + evts = (int) tok.nval; + } + } + return evts; + } + + @SuppressWarnings("static-access") + private static void getToNextEvent(StreamTokenizer tok) throws IOException { + while (tok.nextToken() != tok.TT_EOF) // System.out.println(tok.toString()); + { + if (tok.ttype == tok.TT_WORD && tok.sval.contentEquals("<event>")) { + tok.nextToken();//get to the EOL + return; + } + } + } + + static private double expWeight(double x) { + return Math.exp(-x / _declength); + } + + static private double expWeight(double x, double gamma) { + return Math.exp(-x / (gamma * _declength)); + } + + static private double findMaxWeight() { + Random generator = new Random(); + int ntrials = 100000; + double maxlength = maxLen; + double maxWeight = 0; + for (int i = 0; i < ntrials; i++) { + double x = generator.nextDouble() * maxlength; + double wght = expWeight(x); + if (wght > maxWeight) { + maxWeight = wght; + } + } + + return maxWeight; + } + + static private double getDecayLength(double MaxWeight, double gamma) { + Random generator = new Random(); + double maxlength = maxLen; + double dl = 0; + double draw = generator.nextDouble(); + double tmpwght = 0; + while (tmpwght < draw) { + dl = generator.nextDouble() * maxlength; + tmpwght = expWeight(dl, gamma) / MaxWeight; + } + return dl; + } + + static private double getDecayLength(double MaxWeight) { + Random generator = new Random(); + double maxlength = maxLen; + double dl = 0; + double draw = generator.nextDouble(); + double tmpwght = 0; + while (tmpwght < draw) { + dl = generator.nextDouble() * maxlength; + tmpwght = expWeight(dl) / MaxWeight; + } + return dl; + } + + static private void readLHEEvent(StreamTokenizer tok, double[] beam, int nevhep) throws IOException { + Random generator = new Random(); + getToNextEvent(tok); + List<Double> nums = getNumbersInLine(tok); + + if (nums.size() != 6) { + throw new RuntimeException("Unexpected entry for number of particles"); + } + int nhep = nums.get(0).intValue(); +// System.out.println("Number of particles for event " + nevhep + ": " + nhep); + + double decLen = 0; + double maxWght = 0; + List<Boolean> passEle = new ArrayList<Boolean>(); + List<Boolean> passPos = new ArrayList<Boolean>(); + Hep3Vector pEle = new BasicHep3Vector(); + Hep3Vector pPos = new BasicHep3Vector(); + int i = 0; + int pdgid = 0; + double[] ApMom = {0, 0, 0}; + double ApMass = 0; + double mom[] = {0, 0, 0}; + int charge; + + for (int npart = 0; npart < nhep; npart++) { + List<Double> vals = getNumbersInLine(tok); + if (vals.size() != 13) { + throw new RuntimeException("Unexpected entry for a particle"); + } + if (vals.get(1).intValue() != 9) {//ignore the vertex for now +// int ip = vals.get(0).intValue(); +// if (ip != i + 1) { +// throw new RuntimeException("Particle numbering mismatch"); +// } + + pdgid = vals.get(0).intValue(); +// System.out.println(idhepTmp); + + if (Math.abs(pdgid) == 611) { + + for (int j = 0; j < 3; j++) + mom[j] = vals.get(j + 6); + + Hep3Vector p = rotate(mom[0], mom[1], mom[2]); + Hep3Vector o = rotate(beam[0], beam[1], beam[2]); + charge = 1; + if (pdgid == 611) + charge = -1; + HelixParamCalculator hpc = new HelixParamCalculator(p, o, charge, bField); + double d0 = hpc.getDCA(); + double phi0 = hpc.getPhi0(); + double z0 = hpc.getZ0(); + double slope = hpc.getSlopeSZPlane(); + double R = hpc.getRadius(); + double x0 = hpc.getX0(); + double y0 = hpc.getY0(); + double xc = getxc(R, d0, phi0); + double yc = getyc(R, d0, phi0);; + +// System.out.println(p.toString()); +// System.out.println("d0 = "+d0+"; phi0 = "+phi0+"; z0 = "+z0+"; slope = "+slope+"; R = "+R); +// System.out.println("x0 = "+x0+"; y0 = "+y0+"; xc = "+xc+"; yc = "+yc); + List<Boolean> passLayer = new ArrayList<Boolean>(); + for (int ii = 0; ii < nLayers; ii++) { + double pathL = PathToXPlane(x0, y0, xc, yc, R, x[ii]); +// System.out.println("path length "+pathL); + Hep3Vector posL = PointOnHelix(xc, yc, R, phi0, z0, slope, pathL); +// System.out.println("Position "+posL.toString()); + passLayer.add(inAcceptance(posL, ySize[ii], zSize[ii], zGap[ii])); + } + + if (pdgid == 611) { + passEle = passLayer; + pEle = p; + } else { + pPos = p; + passPos = passLayer; + } + } + + i++; + } + } + if (eventPass(pEle, pPos)) + doAccounting(passEle, passPos, 1000.0 * getInvMass(pEle, pPos)); + + } + + public static Hep3Vector rotate(double x, double y, double z) { + return new BasicHep3Vector(z, x, y); + } + + public static String convertDecimal(String num) { + if (num.contains(".")) { + num = num.replace(".", "pt"); + } + return num; + } + + public static boolean inAcceptance(Hep3Vector position, double yExt, double zExt, double zGap) { + double ypos = position.y(); + if (Math.abs(ypos) > yExt) + return false; + double zpos = position.z(); + if (Math.abs(zpos) < zGap) + return false; + if (Math.abs(zpos) > zGap + zExt) + return false; + + return true; + } + + public static boolean eventPass(Hep3Vector p1, Hep3Vector p2) { + + if (p1.magnitude() + p2.magnitude() < 0.8 * _ecm)//trigger requires 80% of beam energy + return false; + if (p1.magnitude() < _pCut) + return false; + if (p2.magnitude() < _pCut) + return false; + return true; + + } + + public static Double PathToXPlane(double x0, double y0, double xc, double yc, double RC, double x) { + // Create a list to hold the path lengths + Double path; + + double y = yc + Math.signum(RC) * Math.sqrt(RC * RC - Math.pow(x - xc, 2)); +// System.out.println("x = "+x+"; y = "+y); + double s = PathCalc(xc, yc, RC, x0, y0, x, y); + +// System.out.println("PathToXPlane : s = "+s+"; sFromClass = "+sFromClass); + + path = s; + + return path; + } + + private static double PathCalc(double xc, double yc, double RC, double x1, double y1, double x2, double y2) { + // Find the angle between these points measured wrt the circle center + double phi1 = Math.atan2(y1 - yc, x1 - xc); + double phi2 = Math.atan2(y2 - yc, x2 - xc); + double dphi = phi2 - phi1; + // Make sure dphi is in the valid range (-pi, pi) + if (dphi > Math.PI) + dphi -= 2. * Math.PI; + if (dphi < -Math.PI) + dphi += 2. * Math.PI; + // Return the arc length + return -RC * dphi; + } + + public static Hep3Vector PointOnHelix(double xc, double yc, double RC, double phi0, double z0, double slope, double s) { + // Find the azimuthal direction at this path length + double phi = phi0 - s / RC; + // Calculate the position on the helix at this path length + double x = xc - RC * Math.sin(phi); + double y = yc + RC * Math.cos(phi); + double z = z0 + s * slope; + // Return the position as a Hep3Vector + return new BasicHep3Vector(x, y, z); + } + + private static double getxc(double R, double d0, double phi0) { + return (R - d0) * Math.sin(phi0); + } + + private static double getyc(double R, double d0, double phi0) { + return -(R - d0) * Math.cos(phi0); + } + + public static double getInvMass(Hep3Vector p1, Hep3Vector p2) { + double esum = 0.; + double pxsum = 0.; + double pysum = 0.; + double pzsum = 0.; + double me = 0.000511; + // Loop over tracks + + + double p1x = p1.x(); + double p1y = p1.y(); + double p1z = p1.z(); + double p1mag2 = p1x * p1x + p1y * p1y + p1z * p1z; + double e1 = Math.sqrt(p1mag2 + me * me); + pxsum += p1x; + pysum += p1y; + pzsum += p1z; + esum += e1; + + double p2x = p2.x(); + double p2y = p2.y(); + double p2z = p2.z(); + double p2mag2 = p2x * p2x + p2y * p2y + p2z * p2z; + double e2 = Math.sqrt(p2mag2 + me * me); + pxsum += p2x; + pysum += p2y; + pzsum += p2z; + esum += e2; + double psum = Math.sqrt(pxsum * pxsum + pysum * pysum + pzsum * pzsum); + double evtmass = esum * esum - psum * psum; + + if (evtmass > 0) + return Math.sqrt(evtmass); + else + return -99; + } +}
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