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