 |
 |
| Commit in mcd-analysis/src/main/java/org/lcsim/mcd/evtgen on MAIN | |||
| MuonColliderFlukaToStdhepConverter.java | +379 | added 1.1 | |
Class to convert text files output by Takashi from Fluka into stdhep files. One particle per event. Particle / event weights are not currently handled.
diff -N MuonColliderFlukaToStdhepConverter.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MuonColliderFlukaToStdhepConverter.java 23 Sep 2013 23:08:16 -0000 1.1 @@ -0,0 +1,379 @@
+package org.lcsim.mcd.evtgen;
+
+import hep.io.stdhep.StdhepEvent;
+import hep.physics.particle.properties.ParticleType;
+import hep.physics.particle.properties.ParticlePropertyManager;
+import hep.io.stdhep.StdhepEndRun;
+import hep.io.stdhep.StdhepBeginRun;
+import hep.io.stdhep.StdhepWriter;
+import java.io.BufferedReader;
+import java.io.EOFException;
+import java.io.File;
+import java.io.FileInputStream;
+import java.io.IOException;
+import java.io.InputStreamReader;
+
+import hep.physics.particle.properties.ParticlePropertyProvider;
+
+import java.util.Arrays;
+import static java.lang.Math.sqrt;
+
+/**
+ *
+ * @author Norman A. Graf
+ *
+ * @version $Id: MuonColliderFlukaToStdhepConverter.java,v 1.1 2013/09/23 23:08:16 ngraf Exp $
+ */
+public class MuonColliderFlukaToStdhepConverter
+{
+
+ public static void main(String args[]) throws IOException
+ {
+ boolean debug = false;
+
+ // remind user of correct usage
+ if (args.length == 0)
+ {
+ usage();
+ }
+ if (args.length == 1 && args[0].equals("-h"))
+ {
+ usage();
+ }
+
+ File outputDir = new File(".");
+ if (args.length > 1)
+ {
+ outputDir = new File(args[1]);
+ }
+ // check if output directory exists
+ if (!outputDir.exists())
+ {
+ System.out.println("\n\n Directory " + outputDir + " does not exist!");
+ System.exit(1);
+ }
+
+ System.out.println("processing: " + args[0]);
+ File f = new File(args[0]);
+ if (!f.exists())
+ {
+ System.out.println("\n\n File " + f + " does not exist!");
+ System.exit(1);
+ }
+ // file exists, get its name...
+ String stdhepFileName = f.getName();
+ int outPutFileNum = 0;
+ int numEventsInFile = 0;
+ int maxnumEventsInFile = 10000;
+ stdhepFileName += "_" + outPutFileNum + "_" + maxnumEventsInFile + ".stdhep";
+
+ String outputFile = outputDir + "/" + stdhepFileName;
+
+ System.out.println(outputFile);
+
+ /*--------------------------------------------------------*/
+ /* NEVHEP - event number (or some special meaning*/
+ /* (see documentation for details) */
+ /* NHEP - actual number of entries in current */
+ /* event. */
+ /* ISTHEP[IHEP] - status code for IHEP'th entry - see */
+ /* documentation for details */
+ /* IDHEP [IHEP] - IHEP'th particle identifier according*/
+ /* to PDG. */
+ /* JMOHEP[IHEP][0] - pointer to position of 1st mother */
+ /* JMOHEP[IHEP][1] - pointer to position of 2nd mother */
+ /* JDAHEP[IHEP][0] - pointer to position of 1st daughter */
+ /* JDAHEP[IHEP][1] - pointer to position of 2nd daughter */
+ /* PHEP [IHEP][0] - X momentum [Gev/c] */
+ /* PHEP [IHEP][1] - Y momentum [Gev/c] */
+ /* PHEP [IHEP][2] - Z momentum [Gev/c] */
+ /* PHEP [IHEP][3] - Energy [Gev] */
+ /* PHEP [IHEP][4] - Mass[Gev/c^2] */
+ /* VHEP [IHEP][0] - X vertex [mm] */
+ /* VHEP [IHEP][1] - Y vertex [mm] */
+ /* VHEP [IHEP][2] - Z vertex [mm] */
+ /* VHEP [IHEP][3] - production time [mm/c] */
+ /*========================================================*/
+
+ //IstHep convention:
+ // 0 - final state particle if JdaHEP=0
+ // intermediate particle if JdaHEP>0
+ // (NEUGEN extension; was "null")
+ // 1 - final state particle
+ // 2 - intermediate state
+ // 3 - documentation line
+ // 4-10 - reserved for future
+ // 11-200 - reserved for specific model use
+ // 201+ - reserved for users
+
+ int _nmax = 500000;
+ int _eventnum = 0;
+ int _nhep = 0;
+ int[] _isthep = new int[_nmax];
+ int[] _idhep = new int[_nmax];
+ int[] _jmohep = new int[2 * _nmax];
+ int[] _jdahep = new int[2 * _nmax];
+ double[] _phep = new double[5 * _nmax];
+ double[] _vhep = new double[4 * _nmax];
+
+ int _nevents = 0;
+
+ //Dummy values...
+ int nevtreq = 1;
+ int nevtgen = 1;
+ int nevtwrt = 1;
+ float stdecom = 2.F;
+ float stdxsec = 2.F;
+ double stdseed1 = 137.;
+ double stdseed2 = 138.;
+
+ double[] mom = new double[3];
+ double[] pos = new double[3];
+ int[] id2pdg = new int[26];
+ double[] idmasses = new double[26];
+ String[] idnames = new String[26];
+ setupFlukaParticleTable(id2pdg, idmasses, idnames);
+ StdhepWriter w = null;
+ try
+ {
+ w = new StdhepWriter(outputFile, "Imported Stdhep Events v1.0", "From file " + f.getName(), 10);
+ w.setCompatibilityMode(false);
+ } catch (java.io.IOException ex)
+ {
+ System.err.println("Error opening file: " + outputFile);
+ ex.printStackTrace();
+ System.exit(1);
+ }
+ // write a begin run record
+ double energy = 0.;
+ w.writeRecord(new StdhepBeginRun(nevtreq, nevtgen, nevtwrt, stdecom, stdxsec, stdseed1, stdseed2));
+ // now loop over contents of this file...
+ FileInputStream fin2 = new FileInputStream(f);
+ try
+ {
+ BufferedReader myInput = new BufferedReader(new InputStreamReader(fin2));
+// double[] values = new double[8];
+ String thisLine;
+ _nhep = 0; // number of particles in this event
+ // read file line by line, each line represents one particle
+
+ while ((thisLine = myInput.readLine()) != null)
+ {
+ if (debug)
+ {
+ System.out.println(thisLine);
+ }
+ String[] nums = thisLine.trim().split("\\s+");
+ //should be 11 entries per line
+//Every row consists of
+// 0 iregion
+// 1 itype
+// 2 weight
+// 3 energy(GeV)
+// 4 x
+// 5 y
+// 6 z (cm)
+// 7 u
+// 8 v
+// 9 w (direction cosines)
+// 10 time (ns)
+
+ int id = Integer.parseInt(nums[1]);
+ if (id > 25)
+ {
+ System.out.println("Bad particle ID!!!");
+ }
+ int pdgId = id2pdg[id];
+ double mass = idmasses[id];
+ //Fluka is cm, stdhep is mm
+ pos[0] = Double.parseDouble(nums[4]) * 10.;
+ pos[1] = Double.parseDouble(nums[5]) * 10.;
+ pos[2] = Double.parseDouble(nums[6]) * 10.;
+ //Fluka and stdhep in GeV
+ energy = Double.parseDouble(nums[3]);
+ double momentum = sqrt(energy * energy - mass * mass);
+ mom[0] = Double.parseDouble(nums[7]) * momentum;
+ mom[1] = Double.parseDouble(nums[8]) * momentum;
+ mom[2] = Double.parseDouble(nums[9]) * momentum;
+ //Fluka time in ns, stdhep in mm/c
+ // .299792458 m 1 1000 mm
+ // ns * -------------- * --- * --------
+ // ns c m
+ double time = Double.parseDouble(nums[10]) * 299.792458;
+
+ // now populate the HEPEVT "common block"
+ _isthep[_nhep] = 1; // final state particle
+ _idhep[_nhep] = pdgId;
+ _phep[0 + 5 * _nhep] = mom[0]; //px
+ _phep[1 + 5 * _nhep] = mom[1]; //py
+ _phep[2 + 5 * _nhep] = mom[2]; //pz
+ _phep[3 + 5 * _nhep] = energy; //E
+ _phep[4 + 5 * _nhep] = mass; // mass
+ _vhep[0 + 4 * _nhep] = pos[0]; //x
+ _vhep[1 + 4 * _nhep] = pos[1]; //y
+ _vhep[2 + 4 * _nhep] = pos[2]; //z
+ _vhep[3 + 4 * _nhep] = time; //t
+// // increment the number of particles in this event
+ _nhep++;
+ if (_nhep >= _nmax)
+ {
+ throw new RuntimeException("\nAre you sure you want to create an event with more than " + _nmax + " particles?\nIf so, please recompile.");
+ }
+
+ //
+ // Create an StdhepEvent for this event and write it out...
+ //
+ StdhepEvent event = new StdhepEvent(_eventnum++, _nhep, _isthep, _idhep, _jmohep, _jdahep, _phep, _vhep);
+ w.writeRecord(event);
+ if (debug)
+ {
+ System.out.println("writing event " + _eventnum + " with " + energy + " GeV");
+ System.out.println("pos: " + Arrays.toString(pos));
+ System.out.println("mom: " + Arrays.toString(mom));
+ System.out.println("pdg: " + pdgId + " " + idnames[id]);
+ }
+ _nevents++;
+ }
+ // write an end run record
+ w.writeRecord(new StdhepEndRun(nevtreq, nevtgen, nevtwrt, stdecom, stdxsec, stdseed1, stdseed2));
+
+ // close the file
+ try
+ {
+ System.out.println(" Closing file: " + outputFile + " with " + _nevents + " events");
+ w.close();
+ } catch (java.io.IOException ex)
+ {
+ System.err.println("Error closing file: " + outputFile);
+ ex.printStackTrace();
+ System.exit(1);
+ }
+ fin2.close();
+ } catch (EOFException ex)
+ {
+ ex.printStackTrace();
+ }
+ }
+
+ public static void usage()
+ {
+ System.out.println("MuonColliderFlukaToStdhepConverter: \n an application to read in data files from Fluka and convert to stdhep format.\n");
+ System.out.println("Usage: \n\n >> java MuonColliderFlukaToStdhepConverter InputTextFile <output directory> \n");
+ System.out.println(" Where: \n InputTextFile is the input text file to process ");
+ System.out.println(" Writes to the current working directory unless output directory is specified");
+ System.out.println("\n e.g. >> java MuonColliderFlukaToStdhepConverter pairs.txt \n");
+ System.out.println(" will convert the particles listed in decays.txt to the same named decays.stdhep");
+
+ System.exit(0);
+ }
+
+ public static void setupFlukaParticleTable(int[] id2pdg, double[] mass, String[] idnames)
+ {
+ // id conversion
+ /*
+ PROTON 1 Proton 2212
+ APROTON 2 Antiproton -2212
+ ELECTRON 3 Electron 11
+ POSITRON 4 Positron -11
+ NEUTRIE 5 Electron Neutrino 12
+ ANEUTRIE 6 Electron Antineutrino -12
+ PHOTON 7 Photon 22
+ NEUTRON 8 Neutron 2112
+ ANEUTRON 9 Antineutron -2112
+ MUON+ 10 Positive Muon -13
+ MUON- 11 Negative Muon 13
+ KAONLONG 12 Kaon-zero long 130
+ PION+ 13 Positive Pion 211
+ PION- 14 Negative Pion -211
+ KAON+ 15 Positive Kaon 321
+ KAON- 16 Negative Kaon -321
+ LAMBDA 17 Lambda 3122
+ ALAMBDA 18 Antilambda -3122
+ KAONSHRT 19 Kaon zero short 310
+ SIGMA- 20 Negative Sigma 3112
+ SIGMA+ 21 Positive Sigma 3222
+ SIGMAZER 22 Sigma-zero 3212
+ PIZERO 23 Pion-zero 111
+ KAONZERO 24 Kaon-zero 311
+ AKAONZER 25 Antikaon-zero -311
+ Reserved 26 --- ---
+ NEUTRIM 27 Muon neutrino 14
+ ANEUTRIM 28 Muon antineutrino -14
+ Blank 29 --- ---
+ Reserved 30 --- ---
+ ASIGMA- 31 Antisigma-minus -3222
+ ASIGMAZE 32 Antisigma-zero -3212
+ ASIGMA+ 33 Antisigma-plus -3112
+ XSIZERO 34 Xi-zero 3322
+ AXSIZERO 35 Antixi-zero -3322
+ XSI- 36 Negative Xi 3312
+ AXSI+ 37 Positive Xi -3312
+ OMEGA- 38 Omega-minus 3334
+ AOMEGA+ 39 Antiomega -3334
+ Reserved 40 --- ---
+ TAU+ 41 Positive Tau -15
+ TAU- 42 Negative Tau 15
+ NEUTRIT 43 Tau neutrino 16
+ ANEUTRIT 44 Tau antineutrino -16
+ D+ 45 D-plus 411
+ D- 46 D-minus -411
+ D0 47 D-zero 421
+ D0BAR 48 AntiD-zero -421
+ DS+ 49 D_s-plus 431
+ DS- 50 D_s-minus -431
+ LAMBDAC+ 51 Lambda_c-plus 4122
+ XSIC+ 52 Xi_c-plus 4232
+ XSIC0 53 Xi_c-zero 4112
+ XSIPC+ 54 Xi'_c-plus 4322
+ XSIPC0 55 Xi'_c-zero 4312
+ OMEGAC0 56 Omega_c-zero 4332
+ ALAMBDC- 57 Antilambda_c-minus -4122
+ AXSIC- 58 AntiXi_c-minus -4232
+ AXSIC0 59 AntiXi_c-zero -4132
+ AXSIPC- 60 AntiXi'_c-minus -4322
+ AXSIPC0 61 AntiXi'_c-zero -4312
+ AOMEGAC0 62 AntiOmega_c-zero -4332
+ */
+ id2pdg[1] = 2212; // p
+ id2pdg[2] = -2212; // pbar
+ id2pdg[3] = 11; //e-
+ id2pdg[4] = -11; //e+
+ id2pdg[5] = 12; //nue
+ id2pdg[6] = -12; //nuebar
+ id2pdg[7] = 22; //gamma
+ id2pdg[8] = 2112; //n
+ id2pdg[9] = -2112; //nbar
+ id2pdg[10] = -13; //mu+
+ id2pdg[11] = 13; //mu-
+ id2pdg[12] = 130; //K0L
+ id2pdg[13] = 211; //pi+
+ id2pdg[14] = -211; //pi-
+ id2pdg[15] = 321; //K+
+ id2pdg[16] = -321; //K-
+ id2pdg[17] = 3122; //Lambda
+ id2pdg[18] = -3122; //Lambdabar
+ id2pdg[19] = 310; //K0S
+ id2pdg[20] = 3112; // Sigma-
+ id2pdg[21] = -3112; //Sigma+
+ id2pdg[22] = 3212; //Sigma0
+ id2pdg[23] = 111; //pi0
+ id2pdg[24] = 311; //K0
+ id2pdg[25] = -311; //K0bar
+
+//masses
+
+ ParticlePropertyProvider ppp = ParticlePropertyManager.getParticlePropertyProvider();
+ for (int i = 0; i < id2pdg.length; ++i)
+ {
+ ParticleType t = ppp.get(id2pdg[i]);
+// System.out.println("Fluka id= " + i + " pdg= " + id2pdg[i] + " " + t.getName());
+ if (!t.getName().equals("unknown"))
+ {
+ mass[i] = t.getMass();
+ idnames[i] = t.getName();
+// System.out.println(mass[i]);
+ }
+
+ }
+ }
+}
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