hps-java/src/main/java/org/lcsim/hps/util
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;
+ }
+}