diff -u -r1.3 -r1.4
--- ConvertToStdhep.java	5 Dec 2011 19:54:13 -0000	1.3
+++ ConvertToStdhep.java	6 Mar 2012 05:12:15 -0000	1.4
@@ -28,7 +28,7 @@

     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 _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;

@@ -200,8 +200,7 @@

         tok.whitespaceChars(0, ' ');
         tok.eolIsSignificant(true);
 

-        double maxWght = 0;
-        if (expDecay) maxWght = findMaxWeight();

+      

 
         double[] beam = {0, 0, 0, 0};
         int nevhep = 0;

@@ -210,11 +209,9 @@

             beam[1] = sigy * generator.nextGaussian() + _yoff;
             beam[2] = sigz * generator.nextGaussian() + _zoff;
             double tmpDecLen = 0;

-            if (expDecay) tmpDecLen = getDecayLength(maxWght);
-            if (flatDecay) tmpDecLen = generator.nextDouble() * maxLen;
-

+           

             if (fileType.contains("lhe"))

-                writeLHEEvent(tok, beam, tmpDecLen, icross);

+                writeLHEEvent(tok, beam,  icross);

             else if (fileType.contains("dat"))
                 writeDATEvent(tok, beam, tmpDecLen, icross);
 

@@ -290,6 +287,10 @@

     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();

@@ -305,6 +306,19 @@

         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;

@@ -389,7 +403,7 @@

                 if (!expDecay && !flatDecay)
                     vhep[4 * i + 0] = beam[0] + _declength;
                 else {

-                    double totApMom = Math.sqrt(ApMom[0] * ApMom[0] + ApMom[1] * ApMom[1] + ApMom[2] * ApMom[2]);

+                    double totApMom = Math.sqrt(ApMom[0] * ApMom[0] + ApMom[1] * ApMom[1] + ApMom[2] * ApMom[2]);                  

 //                        System.out.println("Decay at : " + tmpDecLen);
                     vhep[4 * i + 0] = beam[0] + decLen * ApMom[0] / totApMom;
                     vhep[4 * i + 1] = beam[1] + decLen * ApMom[1] / totApMom;

@@ -419,7 +433,8 @@

         sw.writeRecord(ev);
     }
 

-    static private void writeLHEEvent(StreamTokenizer tok, double[] beam, double decLen, int nevhep) throws IOException {

+    static private void writeLHEEvent(StreamTokenizer tok, double[] beam,  int nevhep) throws IOException {
+         Random generator = new Random();

         getToNextEvent(tok);
         List<Double> nums = getNumbersInLine(tok);
 

@@ -428,7 +443,10 @@

         int nhep = nums.get(0).intValue();
 //        System.out.println("Number of particles for event " + nevhep + ": " + nhep);
 

-

+        double decLen = 0;
+        double maxWght = 0;
+          
+        if (expDecay) maxWght = findMaxWeight();

         int isthep[] = new int[nhep];
         int idhep[] = new int[nhep];
         int jmohep[] = new int[2 * nhep];

@@ -438,12 +456,13 @@

         int i = 0;
         int idhepTmp = 0;
         double[] ApMom = {0, 0, 0};

+        double ApMass=0;
+        double ApEnergy=0;

         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");

@@ -452,8 +471,6 @@

 
                 if (vals.get(1).intValue() == -1)
                     isthep[i] = 3;

-
-

                 idhepTmp = vals.get(0).intValue();
 //                System.out.println(idhepTmp);
                 idhep[i] = idhepTmp;

@@ -471,21 +488,26 @@

                     ApMom[0] = phep[5 * i + 2];
                     ApMom[1] = phep[5 * i + 1];
                     ApMom[2] = phep[5 * i + 0];

+                    ApMass=vals.get(11);
+                    ApEnergy=vals.get(10);

                 }
                 for (int j = 0; j < 4; j++)
                     vhep[4 * i + j] = beam[j];
 

-                if (!trident && (Math.abs(idhepTmp) == 611 || Math.abs(idhepTmp) == 622))

+                if (!trident && (Math.abs(idhepTmp) == 611 || Math.abs(idhepTmp) == 622)){
+                    double gamma = ApEnergy/ApMass;
+                    if (expDecay) decLen = getDecayLength(maxWght, gamma);
+                     if (flatDecay) decLen = generator.nextDouble() * maxLen;

                     if (!expDecay && !flatDecay)
                         vhep[4 * i + 0] = beam[0] + decLen;

-                    else {

+                    else {                    

                         double totApMom = Math.sqrt(ApMom[0] * ApMom[0] + ApMom[1] * ApMom[1] + ApMom[2] * ApMom[2]);
                         System.out.println("Decay at : " + decLen);
                         vhep[4 * i + 0] = beam[0] + decLen * ApMom[0] / totApMom;
                         vhep[4 * i + 1] = beam[1] + decLen * ApMom[1] / totApMom;
                         vhep[4 * i + 2] = beam[2] + decLen * ApMom[2] / totApMom;
                     }

-

+                }

                 // swap x and z axes...
                 /*Don't do this anymore!  We do stuff in JLAB frame now...
                 double px = phep[5 * i + 0];