mcd-analysis/src/main/java/org/lcsim/mcd/evtgen
diff -N MuonColliderMarsToStdhepConverter.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MuonColliderMarsToStdhepConverter.java 6 Jun 2011 21:16:17 -0000 1.1
@@ -0,0 +1,363 @@
+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 static java.lang.Math.sqrt;
+
+/**
+ *
+ * @author Norman A. Graf
+ *
+ * @version $Id: MuonColliderMarsToStdhepConverter.java,v 1.1 2011/06/06 21:16:17 ngraf Exp $
+ */
+public class MuonColliderMarsToStdhepConverter
+{
+
+ public static void main(String args[]) throws IOException
+ {
+ System.out.println(args.length);
+ // 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[41];
+ double[] idmasses = new double[41];
+ setupParticleTable(id2pdg, idmasses);
+ 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
+ boolean first = true;
+ 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
+ int oldDecayNum = 999;
+ int decayNum = 0;
+ double eSum = 0.;
+ while ((thisLine = myInput.readLine()) != null)
+ {
+ String[] nums = thisLine.trim().split("\\s+");
+ //should be 18 entries per line
+//Every row consists of
+// NI,JJ,X,Y,Z,PX,PY,PZ,TOFF,W,ZDEC,XORIG,YORIG,ZORIG,WORIG,EORIG,IORIG,KORIG
+// FORMAT(I11,I7,3F11.3,11(1PE14.6),I7,I3)
+//
+//
+//0 NI - decay #,
+//1 JJ - particle ID, reaching detector
+//2 X,
+//3 Y,
+//4 Z - coordinates on interface surface (cm)
+//5 PX,
+//6 PY,
+//7 PZ - momentum components (GeV/c)
+//8 TOFF - time with respect to bunch crossing.
+//9 W - particle weight.
+//10 ZDEC - decay coordinate along beam line.
+//11 XORIG,
+//12 YORIG,
+//13 ZORIG - coordinate of particle origin (cm),
+//14 WORIG,
+//15 EORIG,
+//16 IORIG - weight, energy, particle ID of particle origin
+//17 KORIG - type of process in which particle were created
+
+ decayNum = Integer.parseInt(nums[0]);
+ if (decayNum != oldDecayNum)
+ {
+ if (!first)
+ {
+ //
+ // Create an StdhepEvent for old decay event and write it out...
+ //
+ StdhepEvent event = new StdhepEvent(_eventnum++, _nhep, _isthep, _idhep, _jmohep, _jdahep, _phep, _vhep);
+ w.writeRecord(event);
+ System.out.println("writing event for event " + _eventnum + " with decayNum " + oldDecayNum + " with " + _nhep + " particles " + eSum + " GeV");
+ oldDecayNum = decayNum;
+ _nevents++;
+ numEventsInFile++;
+ if (numEventsInFile >= maxnumEventsInFile)
+ {
+ System.out.println("wrote " + (numEventsInFile - 1) + " events to " + outputFile);
+ numEventsInFile = 0;
+ outPutFileNum++;
+ // write an end run record
+ w.writeRecord(new StdhepEndRun(nevtreq, nevtgen, nevtwrt, stdecom, stdxsec, stdseed1, stdseed2));
+ // now close this file
+ w.close();
+ // and open the next...
+ outputFile = outputDir + "/" + f.getName() + "_" + outPutFileNum + "_"+maxnumEventsInFile +".stdhep";
+ System.out.println(outputFile);
+ w = new StdhepWriter(outputFile, "Imported Stdhep Events v1.0", "From file " + f.getName(), 10);
+ w.setCompatibilityMode(false);
+ }
+ _nhep = 0;
+ eSum = 0.;
+ } else
+ {
+ oldDecayNum = decayNum;
+ first = false;
+ }
+ }
+// System.out.println(decayNum);
+// int j = 0;
+ int id = Integer.parseInt(nums[1]);
+ if (id > 31)
+ System.out.println("Bad particle ID!!!");
+ int pdgId = id2pdg[id];
+ double mass = idmasses[id];
+ //MARS is cm, stdhep is mm
+ pos[0] = Double.parseDouble(nums[2]) * 10.;
+ pos[1] = Double.parseDouble(nums[3]) * 10.;
+ pos[2] = Double.parseDouble(nums[4]) * 10.;
+ //MARS and stdhep in GeV
+ mom[0] = Double.parseDouble(nums[5]);
+ mom[1] = Double.parseDouble(nums[6]);
+ mom[2] = Double.parseDouble(nums[7]);
+ //MARS time in s, stdhep in mm/c
+ // 299,792,458 m 1 1000 mm
+ // s * -------------- * --- * --------
+ // s c m
+// double time = Double.parseDouble(nums[8]) * 299792458000.;
+ double time = 0.;
+ double energy = sqrt(mom[0] * mom[0] + mom[1] * mom[1] + mom[2] * mom[2] + mass * mass);
+ eSum += energy;
+ // 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 old decay event and write it out...
+ //
+ StdhepEvent event = new StdhepEvent(_eventnum++, _nhep, _isthep, _idhep, _jmohep, _jdahep, _phep, _vhep);
+ w.writeRecord(event);
+ System.out.println("writing event for event " + _eventnum + " with decayNum " + oldDecayNum + " with " + _nhep + " particles " + eSum + " GeV");
+ // 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("MuonColliderMarsToStdhepConverter: \n an application to read in data files from MARS and convert to stdhep format.\n");
+ System.out.println("Usage: \n\n >> java MuonColliderMarsToStdhepConverter 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 MuonColliderMarsToStdhepConverter 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 setupParticleTable(int[] id2pdg, double[] mass)
+ {
+ // id conversion
+// PARTICLE ID:
+// 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
+// p n pi+ pi- K+ K- mu+ mu- g e- e+ pbar pi0 d t He3 He4 num nuam nue nuae
+// 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
+// K0L K0S K0 AK0 LAM ALAM SIG+ SIG0 SIG- nbar KSI0 KSI- OM- sib- si0b sib+ AKSI0 ksib+ omb+
+ id2pdg[1] = 2212; // p
+ id2pdg[2] = 2112; // n
+ id2pdg[3] = 211; //pi+
+ id2pdg[4] = -211; //pi-
+ id2pdg[5] = 321; //K+
+ id2pdg[6] = -321; //K-
+ id2pdg[7] = -13; //mu+
+ id2pdg[8] = 13; //mu-
+ id2pdg[9] = 22; //g
+ id2pdg[10] = 11; //e-
+ id2pdg[11] = -11; //e+
+ id2pdg[12] = -2212; //pbar
+ id2pdg[13] = 111; //pi0
+ id2pdg[14] = 1000010020; //d
+ id2pdg[15] = 1000010030; //t
+ id2pdg[16] = 1000020030; //He3
+ id2pdg[17] = 1000020040; //He4
+ id2pdg[18] = 14; // num
+ id2pdg[19] = -14; //nuam
+ id2pdg[20] = 12; // nu1
+ id2pdg[21] = -12; //nuae
+ id2pdg[22] = 130; //K0L
+ id2pdg[23] = 310; //K0S
+ id2pdg[24] = 311; //K0
+ id2pdg[25] = -311; //AK0
+ id2pdg[26] = 3122; //LAM
+ id2pdg[27] = -3122; //ALAM
+ id2pdg[28] = 3222; //SIG+
+ id2pdg[29] = 3212; //SIG0
+ id2pdg[30] = 3112; //SIG-
+ id2pdg[31] = -2112; //nbar
+ id2pdg[32] = 3322;
+ id2pdg[33] = 3312;
+ id2pdg[34] = 3334;
+ id2pdg[35] = -3222;
+ id2pdg[36] = -3212;
+ id2pdg[37] = -3112;
+ id2pdg[38] = -3322;
+ id2pdg[39] = -3312;
+ id2pdg[40] = -3334;
+//masses
+
+ ParticlePropertyProvider ppp = ParticlePropertyManager.getParticlePropertyProvider();
+ for (int i = 0; i < id2pdg.length; ++i)
+ {
+ ParticleType t = ppp.get(id2pdg[i]);
+ System.out.println("i= " + i + " pdg= " + id2pdg[i] + " " + t.getName());
+ if (!t.getName().equals("unknown"))
+ {
+ mass[i] = t.getMass();
+// System.out.println(mass[i]);
+ }
+
+ }
+ }
+}
