 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/util on MAIN | |||
| CalculateAcceptanceFromMadGraph.java | +196 | -99 | 1.2 -> 1.3 |
fix a dumb mistake in accounting
diff -u -r1.2 -r1.3 --- CalculateAcceptanceFromMadGraph.java 5 Sep 2012 19:41:04 -0000 1.2 +++ CalculateAcceptanceFromMadGraph.java 8 Jan 2013 21:04:14 -0000 1.3 @@ -5,34 +5,23 @@
package org.lcsim.hps.util; /**
- - @author mgraham
+ * + * @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.BasicHepLorentzVector;
import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.HepLorentzVector;
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 {
@@ -49,83 +38,121 @@
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 = 1000; //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 boolean _isMuon=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[] ySize = {sensorLength / 2, sensorLength / 2, sensorLength / 2, sensorLength / 2, sensorLength / 2, sensorLength / 2, sensorLength / 2};
+ static double[] ySize = {sensorLength / 2, sensorLength / 2, sensorLength / 2, sensorLength/2 , sensorLength / 2, sensorLength , sensorLength };
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 double[] zGap = {x[0] * gap, x[1] * gap, x[2] * gap, x[3] * gap, x[4] * gap, x[5] * gap, x[6] * gap};
+ static double[] ySizeFull = {sensorWidth, sensorWidth, sensorWidth, 3 * sensorWidth / 2, 2 * sensorWidth, 7 * sensorWidth / 2, 4 * sensorLength,};
+ static double[] zSizeFull = {sensorWidth, sensorWidth, 2 * sensorWidth, 2 * sensorWidth, sensorLength, sensorLength, sensorLength};
+ static double[] zGapFull = {x[0] * 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 nPassFull = 0;
static int nPassTestLyr50 = 0;
- static boolean elePass; - static boolean posPass; - static boolean recoilPass; - static int nbins = 160; - static double maxMass = 400.0;
+ /* + static double maxMass = 1000.0; + static double bField = 1.0; + static double _pCut = 0.25; + static double _ecm = 6.6; //GeV + */ + static double maxMass = 1000.0; + static double bField = 1.0; + static double _pCut = 0.05; + static double _ecm = 4.4; //GeV + static double binSize = 1.0; + static int nbins = (int) (maxMass / binSize);
static int[] genMass = new int[nbins];
static int[] recoMassLyr1 = new int[nbins];
+ static int[] recoMassLyr1Full = new int[nbins];
static int[] recoMassLyr50 = new int[nbins];
- public static void doAccounting(List<Boolean> passEle, List<Boolean> passPos, boolean passEvt, double invMass) {
+ public static void doAccounting(List<Boolean> passEle, List<Boolean> passPos, List<Boolean> passEleFull, List<Boolean> passPosFull, boolean passEvt, double invMass) {
+
+ boolean passLyr700Full = false;
+ boolean passLyr100Full = false;
+
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);
- if(bin<nbins){
- 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]++;
- }else{
- System.out.println("Mass out of range! "+invMass);
+// System.out.println("invMass = "+invMass + "...goes in bin #"+bin);
+ if (bin < nbins) {
+ genMass[bin]++;
+ nGenerated++;
+ //see if pos+ele is in layer 70cm layer (required for both test & full)
+ if (passEle.get(5) && passPos.get(5)) {
+ passLyr700 = true;
+ }
+ if (passEleFull.get(5) && passPosFull.get(5)) {
+ passLyr700Full = true;
+ }
+
+ //see if they were in layer 1
+ if (passEle.get(1) && passPos.get(1)) {
+ passLyr100 = true;
+ }
+ //see if they were in layer 1
+ if (passEleFull.get(1) && passPosFull.get(1)) {
+ passLyr100Full = 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 (passLyr700Full && passLyr100Full && passEvt) {
+ recoMassLyr1Full[bin]++;
+ nPassFull++;
+ }
+ //ok...fill the histograms
+ if (passLyr700 && passLyr100 && passEvt) {
+ recoMassLyr1[bin]++;
+ nPassTest++;
+ }
+ if (passLyr700 && passLyr50 && passEvt) { //require they go through 6 layers...
+ recoMassLyr50[bin]++;
+ nPassTestLyr50++;
+ }
+
+ } else {
+// System.out.println("Mass out of range! "+invMass);
}
-
+
}
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."));
options.addOption(new Option("b", true, "B-Field"));
+ options.addOption(new Option("t", true, "Rad, BH, or FullRadBH"));
+ options.addOption(new Option("u", false, "Is muon decay?"));
return options;
}
/**
- @param args the command line arguments - @throws IOException
+ * @param args the command line arguments + * @throws IOException
*/
public static void main(String[] args) throws IOException {
@@ -143,18 +170,22 @@
} catch (ParseException e) {
throw new RuntimeException("Problem parsing command line options.", e);
}
-
+
String ninString = String.valueOf(nInFiles);
String ecmString = String.valueOf(_ecm);
+ String bString = "0.5";
String eptString = String.valueOf(_ecm);
+ String typeString = "Rad";
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");
@@ -162,20 +193,31 @@
System.out.println(ecmString);
eptString = convertDecimal(ecmString);
}
+ if (cl.hasOption("b")) {
+ bString = cl.getOptionValue("b");
+ bField = Double.valueOf(bString);
+ System.out.println(bString);
+ }
- -
+ if (cl.hasOption("t")) {
+ typeString = cl.getOptionValue("t");
+ System.out.println(typeString);
+ }
+
+ if (cl.hasOption("u")) {
+ _isMuon=true;
+ }
- String fDir = "/nfs/slac/g/hps/mgraham/DarkPhoton/MadGraph/Events" + eptString + "Rad/";
+ String fDir = "/nfs/slac/g/hps/mgraham/DarkPhoton/MadGraph/Events" + eptString + typeString + "/";
// String fDir = "./"; // String inLabel = "W" + eptString + "GeV_Ap" + massString + "MeV_";
- String inLabel = "W" + eptString + "GeV_Rad_";
+ String inLabel = "W" + eptString + "GeV_" + typeString + "_";
String inPost = "_unweighted_events.lhe";
// String outDir = "/nfs/slac/g/hps/mgraham/DarkPhoton/SignalEvents/";
- String outDir = "./";
+ String outDir = "./Acceptance/";
for (int i = 0; i < nInFiles; i++) {
int fnum = i + 1;
String snum = "_" + fnum;
@@ -226,16 +268,29 @@
}
//ok ... now spit out histograms
+ System.out.println("nGenerated = " + nGenerated);
+ System.out.println("nPass(Test) = " + nPassTest);
+ System.out.println("nPass(Test Lyr50)= " + nPassTestLyr50);
+ System.out.println("nPass(Full) = " + nPassFull);
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]);
+ System.out.printf("%d\t%4.1f\t%d\t%d\t%d\t%d\n", kk, mass, genMass[kk], recoMassLyr1[kk], recoMassLyr50[kk], recoMassLyr1Full[kk]);
}
+ /* + outputFile(outDir + typeString + eptString + "_Test_1T.dat", recoMassLyr1); + outputFile(outDir + typeString + eptString + "_Full_1T.dat", recoMassLyr1Full); + outputFile(outDir + typeString + eptString + "_lyr50_1T.dat", recoMassLyr50); + outputFile(outDir + typeString + eptString + "_generated_1T.dat", genMass); + */ + + outputFile(outDir + typeString + eptString + "_Test.dat", recoMassLyr1); + outputFile(outDir + typeString + eptString + "_Full.dat", recoMassLyr1Full); + outputFile(outDir + typeString + eptString + "_lyr50.dat", recoMassLyr50); + outputFile(outDir + typeString + eptString + "_generated.dat", genMass);
- outputFile(outDir + "rad" + eptString + "_lyr1.dat", recoMassLyr1); - outputFile(outDir + "rad" + eptString + "_lyr50.dat", recoMassLyr50);
}
@@ -474,19 +529,25 @@
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();
+ List<Boolean> passEleFull = new ArrayList<Boolean>(); + List<Boolean> passPosFull = new ArrayList<Boolean>(); + List<Boolean> passRecoil = new ArrayList<Boolean>(); + List<Boolean> passRecoilFull = new ArrayList<Boolean>(); + //Hep3Vector pEle = new BasicHep3Vector(); + //Hep3Vector pPos = new BasicHep3Vector(); + //Hep3Vector pRecoil = new BasicHep3Vector(); + HepLorentzVector pEle = new BasicHepLorentzVector(); + HepLorentzVector pPos = new BasicHepLorentzVector(); + HepLorentzVector pRecoil = new BasicHepLorentzVector();
int i = 0;
int pdgid = 0;
double[] ApMom = {0, 0, 0};
- double ApMass = 0;
double mom[] = {0, 0, 0};
+ double ene=0;
int charge;
-
+ boolean foundRecoil = false;
for (int npart = 0; npart < nhep; npart++) {
List<Double> vals = getNumbersInLine(tok);
if (vals.size() != 13) {
@@ -500,17 +561,22 @@
pdgid = vals.get(0).intValue();
// System.out.println(idhepTmp);
-
- if (Math.abs(pdgid) == 611) {
-
- for (int j = 0; j < 3; j++)
+// System.out.println(pdgid+" "+vals.get(1).intValue());
+ if (Math.abs(pdgid) == 611 || (Math.abs(pdgid) == 11 && vals.get(1).intValue() == 1)) {
+// if (Math.abs(pdgid) == 611) {
+ // System.out.println("Ok...getting info for this particle");
+ for (int j = 0; j < 3; j++) {
mom[j] = vals.get(j + 6);
-
+ } + ene= vals.get(10);
Hep3Vector p = rotate(mom[0], mom[1], mom[2]);
+ // Hep3Vector p = rotate(mom[1], mom[0], mom[2]); //flip x,y because my trident files have cut on thetaX that may bias things
Hep3Vector o = rotate(beam[0], beam[1], beam[2]);
charge = 1;
- if (pdgid == 611)
+ if (pdgid == 611) {
charge = -1;
+ } + HepLorentzVector pl=new BasicHepLorentzVector(ene,p);
HelixParamCalculator hpc = new HelixParamCalculator(p, o, charge, bField);
double d0 = hpc.getDCA();
double phi0 = hpc.getPhi0();
@@ -522,33 +588,50 @@
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);
+// 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>();
+ List<Boolean> passLayerFull = 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]));
+ passLayerFull.add(inAcceptance(posL, ySizeFull[ii], zSizeFull[ii], zGapFull[ii]));
}
- if (pdgid == 611) {
+ if (pdgid == 611) { //electron from A'
passEle = passLayer;
- pEle = p;
- } else {
- pPos = p;
+ passEleFull = passLayerFull;
+ pEle = pl;
+ } else if (pdgid == -611||pdgid==-11) {//positron
+ pPos = pl;
passPos = passLayer;
+ passPosFull = passLayerFull;
+ } else if (pdgid == 11) {
+ if (!foundRecoil) {
+ foundRecoil=true;
+ passRecoil = passLayer;
+ passRecoilFull = passLayerFull;
+ pRecoil = pl;
+ }
+// else{
+// passEle = passLayer;
+// passEleFull = passLayerFull;
+// pEle = p;
+// }
}
+
}
i++;
}
}
- - doAccounting(passEle, passPos, eventPass(pEle,pPos),1000.0 * getInvMass(pEle, pPos)); -
+ double invMass=getInvMass(pEle.v3(), pPos.v3()); + doAccounting(passEle, passPos, passEleFull, passPosFull, eventPass(pEle.v3(), pPos.v3()), 1000.0 * invMass); + // doAccounting(passRecoil, passPos, passRecoilFull, passPosFull, eventPass(pRecoil, pPos), 1000.0 * getInvMass(pRecoil, pPos));
}
public static Hep3Vector rotate(double x, double y, double z) {
@@ -564,29 +647,38 @@
public static boolean inAcceptance(Hep3Vector position, double yExt, double zExt, double zGap) {
double ypos = position.y();
- if (Math.abs(ypos) > yExt)
+ if (Math.abs(ypos) > yExt) {
return false;
+ }
double zpos = position.z();
- if (Math.abs(zpos) < zGap)
+ if (Math.abs(zpos) < zGap) {
return false;
- if (Math.abs(zpos) > zGap + zExt)
+ }
+ if (Math.abs(zpos) > zGap + zExt) {
return false;
+ }
return true;
}
public static boolean eventPass(Hep3Vector p1, Hep3Vector p2) {
- // System.out.println("p1.magnitude = "+p1.magnitude()+"; p2.magnitude = "+p2.magnitude()+"; 0.8*_ecm = "+0.8*_ecm);
-
+ // System.out.println("p1.magnitude = "+p1.magnitude()+"; p2.magnitude = "+p2.magnitude()+"; 0.8*_ecm = "+0.8*_ecm);
+
if (p1.magnitude() + p2.magnitude() < 0.8 * _ecm)//trigger requires 80% of beam energy
+ {
return false;
+ }
// System.out.println("Passed totenergy");
- if (p1.magnitude() < _pCut)
+ if (p1.magnitude() < _pCut) {
return false;
- if (p2.magnitude() < _pCut)
+ }
+ if (p2.magnitude() < _pCut) {
return false;
- if(p2.z()*p1.z()>0) // this is basically the opposite quadrant cut in the trigger (B-field makes them opposite in y)
+ }
+ if (p2.z() * p1.z() > 0) // this is basically the opposite quadrant cut in the trigger (B-field makes them opposite in y)
+ {
return false;
+ }
// System.out.println("Event is good!!!!!");
return true;
@@ -613,10 +705,12 @@
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)
+ if (dphi > Math.PI) {
dphi -= 2. * Math.PI;
- if (dphi < -Math.PI)
+ }
+ if (dphi < -Math.PI) {
dphi += 2. * Math.PI;
+ }
// Return the arc length
return -RC * dphi;
}
@@ -646,6 +740,8 @@
double pysum = 0.;
double pzsum = 0.;
double me = 0.000511;
+ if(_isMuon) + me=0.1057;
// Loop over tracks
@@ -671,9 +767,10 @@
double psum = Math.sqrt(pxsum * pxsum + pysum * pysum + pzsum * pzsum);
double evtmass = esum * esum - psum * psum;
- if (evtmass > 0)
+ if (evtmass > 0) {
return Math.sqrt(evtmass);
- else
+ } 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