Commit in hps-java/src/main/java/org/lcsim/hps/util on MAIN | |||
CalculateAcceptanceFromMadGraph.java | +402 | -447 | 1.3 -> 1.4 |
Remove test-run geometry; add high acceptance geometry.
diff -u -r1.3 -r1.4 --- CalculateAcceptanceFromMadGraph.java 8 Jan 2013 21:04:14 -0000 1.3 +++ CalculateAcceptanceFromMadGraph.java 2 May 2013 22:24:12 -0000 1.4 @@ -5,8 +5,8 @@
package org.lcsim.hps.util; /**
- * - * @author mgraham
+ + @author mgraham
*/ import hep.physics.vec.BasicHep3Vector; import hep.physics.vec.BasicHepLorentzVector;
@@ -24,122 +24,109 @@
import org.apache.commons.cli.PosixParser; import org.lcsim.fit.helicaltrack.HelixParamCalculator;
-public class CalculateAcceptanceFromMadGraph {
+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 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 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 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 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 , sensorLength }; - static double[] zSize = {sensorWidth, sensorWidth, sensorWidth, sensorWidth, sensorWidth, sensorWidth, sensorWidth}; - 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 int _nEleRequired=2; + 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, sensorLength}; + static double[] zSize={sensorWidth, sensorWidth, sensorWidth, sensorWidth, sensorWidth, sensorWidth, sensorWidth}; + 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 nPassNoMax=0; + static int nPassFull=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, 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;
+ 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[] recoMassLyr1NoMax=new int[nbins]; + static int[] recoMassLyr50Full=new int[nbins]; + + public static void doAccounting(List<Boolean> passEleNoMax, List<Boolean> passPosNoMax, List<Boolean> passEleFull, List<Boolean> passPosFull, boolean passEvt, double invMass){ + + boolean passLyr700Full=false; + boolean passLyr100Full=false; + boolean passLyr50Full=false; + boolean passLyr700NoMax=false; + boolean passLyr100NoMax=false; +
//find the bin...
- int bin = (int) (invMass / maxMass * nbins);
+ int bin=(int) (invMass/maxMass*nbins);
// System.out.println("invMass = "+invMass + "...goes in bin #"+bin);
- if (bin < nbins) {
+ 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; - }
+ + if(passEleFull.get(5)&&passPosFull.get(5)) + passLyr700Full=true; + if(passEleNoMax.get(5)&&passPosNoMax.get(5)) + passLyr700NoMax=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; - }
+ if(passEleFull.get(1)&&passPosFull.get(1)) + passLyr100Full=true; + + if(passEleNoMax.get(1)&&passPosNoMax.get(1)) + passLyr100NoMax=true; + + if(passEleFull.get(0)&&passPosFull.get(0)) + passLyr50Full=true;
//ok...fill the histograms
- if (passLyr700Full && passLyr100Full && passEvt) {
+ 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++;
+ if(passLyr700NoMax&&passLyr100NoMax&&passEvt){ + recoMassLyr1NoMax[bin]++; + nPassNoMax++;
}
+ +
- } else {
+ } else{
// System.out.println("Mass out of range! "+invMass); } }
- private static Options createCommandLineOptions() { - Options options = new Options();
+ 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."));
@@ -151,133 +138,128 @@
} /**
- * @param args the command line arguments - * @throws IOException
+ @param args the command line arguments + @throws IOException
*/
- public static void main(String[] args) throws IOException {
+ public static void main(String[] args) throws IOException{
// Set up command line parsing.
- Options options = createCommandLineOptions();
+ Options options=createCommandLineOptions();
- CommandLineParser parser = new PosixParser();
+ CommandLineParser parser=new PosixParser();
// Parse command line arguments.
- CommandLine cl = null; - try { - cl = parser.parse(options, args);
+ CommandLine cl=null; + try{ + cl=parser.parse(options, args);
System.out.println("Trying parser");
- } catch (ParseException e) {
+ } 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);
+ 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);
+ 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(ecmString);
+ + + if(cl.hasOption("e")){ + ecmString=cl.getOptionValue("e"); + _ecm=Double.valueOf(ecmString);
System.out.println(ecmString);
- eptString = convertDecimal(ecmString);
+ eptString=convertDecimal(ecmString);
}
- if (cl.hasOption("b")) { - bString = cl.getOptionValue("b"); - bField = Double.valueOf(bString); - System.out.println(bString);
+ if(cl.hasOption("b")){ + bString=cl.getOptionValue("b"); + bField=Double.valueOf(bString); + System.out.println(bString);
}
- if (cl.hasOption("t")) { - typeString = cl.getOptionValue("t");
+ if(cl.hasOption("t")){ + typeString=cl.getOptionValue("t");
System.out.println(typeString); }
- - if (cl.hasOption("u")) {
+ + if(cl.hasOption("u"))
_isMuon=true;
- }
- String fDir = "/nfs/slac/g/hps/mgraham/DarkPhoton/MadGraph/Events" + eptString + typeString + "/";
+ 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_" + typeString + "_";
+ String inLabel="W"+eptString+"GeV_"+typeString+"_";
- String inPost = "_unweighted_events.lhe";
+ String inPost="_unweighted_events.lhe";
// String outDir = "/nfs/slac/g/hps/mgraham/DarkPhoton/SignalEvents/";
- String outDir = "./Acceptance/"; - for (int i = 0; i < 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 outDir="./Acceptance/"; + for(int i=0; i<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")) {
+ 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 {
+ } 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) {
+ } 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 {
+ 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) {
+ } catch(InterruptedException ex){
Logger.getLogger(CalculateAcceptanceFromMadGraph.class.getName()).log(Level.SEVERE, null, ex); } } } //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("nGenerated = "+nGenerated); +// System.out.println("nPass(Test) = " + nPassTest); +// System.out.println("nPass(Test Lyr50)= " + nPassTestLyr50); + System.out.println("nPass(Full) = "+nPassFull); + System.out.println("nPass(NoMax) = "+nPassNoMax);
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;
+ 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\t%d\n", kk, mass, genMass[kk], recoMassLyr1[kk], recoMassLyr50[kk], recoMassLyr1Full[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]);
+ System.out.printf("%d\t%4.1f\t%d\t%d\t%d\t%d\n", kk, mass, genMass[kk], recoMassLyr1Full[kk], recoMassLyr50Full[kk], recoMassLyr1Full[kk]);
} /* outputFile(outDir + typeString + eptString + "_Test_1T.dat", recoMassLyr1);
@@ -286,74 +268,71 @@
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 + typeString + eptString + "_Test.dat", recoMassLyr1); + outputFile(outDir+typeString+eptString+"_Full.dat", recoMassLyr1Full); + outputFile(outDir+typeString+eptString+"_NoMax.dat", recoMassLyr1NoMax); + outputFile(outDir+typeString+eptString+"_lyr50.dat", recoMassLyr50Full); + outputFile(outDir+typeString+eptString+"_generated.dat", genMass);
} @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) {
+ 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")) {
+ 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++) {
+ 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")) {
+ 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); - }
+ 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();
+ private static int process(String infile) throws IOException{ + Random generator=new Random();
- FileReader lc = new FileReader(infile); - StreamTokenizer lctok = new StreamTokenizer(lc);
+ 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);
+ int nevts=getnevts(lctok);
lc.close();
- System.out.println("Found " + nevts + " events");
+ System.out.println("Found "+nevts+" events");
- FileReader fr = new FileReader(infile);
+ FileReader fr=new FileReader(infile);
- StreamTokenizer tok = new StreamTokenizer(fr);
+ StreamTokenizer tok=new StreamTokenizer(fr);
tok.resetSyntax(); tok.wordChars(33, 255);
@@ -363,13 +342,13 @@
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;
+ 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);
@@ -387,19 +366,19 @@
- 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[] 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;
+ double tmpDecLen=0;
readLHEEvent(tok, beam, icross);
@@ -415,13 +394,12 @@
} @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) {
+ 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;
+ String tokVal=tok.sval;
// System.out.println(tokVal); nums.add(Double.valueOf(tokVal).doubleValue()); }
@@ -430,347 +408,324 @@
} @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) {
+ 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) {
+ if(tok.ttype==tok.TT_EOL)
break;
- } - if (tok.ttype == tok.TT_WORD && tok.sval.contentEquals("Number")) {
+ 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; - }
+ 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>")) {
+ 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){ + return Math.exp(-x/_declength);
}
- static private double expWeight(double x, double gamma) { - return Math.exp(-x / (gamma * _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; - }
+ 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;
+ 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;
+ 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();
+ static private void readLHEEvent(StreamTokenizer tok, double[] beam, int nevhep) throws IOException{ + Random generator=new Random();
getToNextEvent(tok);
- List<Double> nums = getNumbersInLine(tok);
+ List<Double> nums=getNumbersInLine(tok);
- if (nums.size() != 6) {
+ if(nums.size()!=6)
throw new RuntimeException("Unexpected entry for number of particles");
- } - int nhep = nums.get(0).intValue();
+ int nhep=nums.get(0).intValue();
// System.out.println("Number of particles for event " + nevhep + ": " + nhep);
- List<Boolean> passEle = new ArrayList<Boolean>(); - List<Boolean> passPos = new ArrayList<Boolean>(); - 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();
+ List<Boolean> passEleNoMax=new ArrayList<Boolean>(); + List<Boolean> passPosNoMax=new ArrayList<Boolean>(); + List<Boolean> passEleFull=new ArrayList<Boolean>(); + List<Boolean> passPosFull=new ArrayList<Boolean>(); + List<Boolean> passRecoilNoMax=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 mom[] = {0, 0, 0};
+ HepLorentzVector pEle=new BasicHepLorentzVector(); + HepLorentzVector pPos=new BasicHepLorentzVector(); + HepLorentzVector pRecoil=new BasicHepLorentzVector(); + int i=0; + int pdgid=0; + double[] ApMom={0, 0, 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) {
+ boolean foundRecoil=false; + 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
+ 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();
+ pdgid=vals.get(0).intValue();
// System.out.println(idhepTmp); // 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||(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) { - charge = -1; - } - HepLorentzVector pl=new BasicHepLorentzVector(ene,p); - 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("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) + charge=-1; + HepLorentzVector pl=new BasicHepLorentzVector(ene, p); + 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>(); - List<Boolean> passLayerFull = new ArrayList<Boolean>(); - for (int ii = 0; ii < nLayers; ii++) { - double pathL = PathToXPlane(x0, y0, xc, yc, R, x[ii]);
+ List<Boolean> passLayerNoMax=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);
+ 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]));
+ // passLayer.add(inAcceptance(posL, ySize[ii], zSize[ii], zGap[ii])); + passLayerFull.add(inAcceptance(posL, ySize[ii], zSize[ii], zGap[ii])); + passLayerNoMax.add(inAcceptance(posL, 9999999, 9999999, zGap[ii]));
}
- if (pdgid == 611) { //electron from A' - passEle = passLayer; - 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; -// } - } -
+ if(pdgid==611){ //electron from A' + passEleNoMax=passLayerNoMax; + passEleFull=passLayerFull; + pEle=pl; + } else if(pdgid==-611||pdgid==-11){//positron + pPos=pl; + passPosNoMax=passLayerNoMax; + passPosFull=passLayerFull; + } else if(pdgid==11) + if(!foundRecoil){ + foundRecoil=true; + passRecoilNoMax=passLayerNoMax; + passRecoilFull=passLayerFull; + pRecoil=pl; + } // else{ + // passEle = passLayer; + // passEleFull = passLayerFull; + // pEle = p; + // } +
} i++; } } 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));
+ doAccounting(passEleNoMax, passPosNoMax, 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) {
+ 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"); - }
+ 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) {
+ 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) {
+ double zpos=position.z(); + 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) {
+ 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);
- if (p1.magnitude() + p2.magnitude() < 0.8 * _ecm)//trigger requires 80% of beam energy - {
+ 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; }[truncated at 1000 lines; 142 more skipped]
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