Merge lp:~maddevelopers/mg5amcnlo/upgrade_pythia_compatibility into lp:mg5amcnlo/lts

[Olivier Mattelaer]

Hi,

This is the first of the big changes in MG5.

This branch has quite a huge number of modifications into it.

This change quite deeply the cmd_interface by splitting/refactoring it in two file.
One source of problem during merges is the function split_arg which is now inside the class
=> self.split_arg

For the Madevent Template all the binaries are now in ./bin/internal!
A couple of them are simply removed and replace by the madevent interface.

This introduces also the launch -i which might be nice to have also in your development.

Please don't wait too long before proposing your merges.
The sooner we all agree on a structure, the better it is for all our project.
This merge doesn't need to be the final one, neither have all functionalities.
The idea is to fix the structure and move together.

My advice for this merge is that you focus on your branch.
1) Make the merge with your branch and check that everything is still working.
 (that's the time to add tests if you don't have enough)
2) If smtg goes wrong/if you notice a change of structure that you don't like, please discuss it here.
3) If you would like a modification to prepare your branch, that's also the time.

Please don't put this on too low priority. I'd say that this should be done in a matter of two weeks. I'm there to help you for the merge and to help you on every step.

Cheers,

Olivier

[lp:~maddevelopers/mg5amcnlo/upgrade_pythia_compatibility]

307. By mattelaer-olivier on 2011-11-29

Correct the bug reported by Sho

308. By mattelaer-olivier on 2011-11-29

improve crossxhtml for decqy
run automaticaly calculate_width in case of 1>X process

309. By mattelaer-olivier on 2011-11-30

correct sho problem (recurence)

310. By Olivier Mattelaer on 2011-12-05

use the tag for the pythia/pgs/delphes in the links and reorganize the output dir

311. By Olivier Mattelaer on 2011-12-06

improve auto-completion
pass sum_html in python
new organisation for the HTML/Event directory
Use the tag name for the event file for pythia/pgs/delphes

312. By mattelaer-olivier on 2011-12-07

Modify gencrossxhtml in order to prepare the possibiity of multi-tag for a run.

313. By Olivier Mattelaer on 2011-12-07

add css files

314. By Olivier Mattelaer on 2011-12-07

correct small bug + prepare the multi-pythia/pgs/delphes

315. By Olivier Mattelaer on 2011-12-09

new module for the banner
multi-pythia working
reintroduce the link for error log in the html

316. By Olivier Mattelaer on 2011-12-10

multi pythia/pgs/delphes is now fully working.

317. By Olivier Mattelaer on 2011-12-10

make the specification of the tag works + auto-completion + help

318. By Olivier Mattelaer on 2011-12-11

Have a remove command supporting the --tag argument.
Have a correct computation of the pythia error.

319. By Olivier Mattelaer on 2011-12-12

make banner_run compatible with tag (and fix some bug associate to the previous version)
improve customize auto-completion (avoid redondancy + fix wrong column spliting)

320. By mattelaer-olivier on 2011-12-12

have a nice history even when using ; trick
if set command is used in mg5 command pass it to Madevent command
update/fix/improve the acceptance tests

321. By Olivier Mattelaer on 2011-12-12

fix launch command (for SA and pythia8)

322. By Olivier Mattelaer on 2011-12-13

merge with rambo_bug branch
fix Johan's bugs
  - problem with autolinking of ME5_debug
  - problem with the MGMEVersion file

323. By Olivier Mattelaer on 2011-12-13

Add the fermion flow check when loading a new model.

324. By Olivier Mattelaer on 2011-12-15

merge with modif on the other computer
fix case sensitive path in acceptance tests
remove 4fermion model (tch)

325. By Olivier Mattelaer on 2011-12-15

change date

326. By Olivier Mattelaer on 2011-12-15

forget to add aloha_fct file

327. By Olivier Mattelaer on 2011-12-15

add results.dat for P directory

328. By Olivier Mattelaer on 2011-12-16

modify collect_result in order
1) to have it running event in groupsubprocess
2) to have it called automatically with the launch command

329. By Olivier Mattelaer on 2011-12-17

change launch scheme for 1>X

330. By Olivier Mattelaer on 2011-12-17

allow arguments in program definitions bug #905053

331. By Olivier Mattelaer on 2011-12-21

allow groupsubprocess with auto value
create numbdiag for the web interface

332. By Olivier Mattelaer on 2011-12-21

merge with 1.3.32

333. By Olivier Mattelaer on 2012-01-04

Fixed problem with acceptance tests (problem with grouping)

334. By Olivier Mattelaer on 2012-01-04

Add the display options output to the MG5_debug/ME5_debug.

335. By Olivier Mattelaer on 2012-01-05

include Johan's comments

336. By Olivier Mattelaer on 2012-01-09

fix plot for pgs in multiple pgs run
fix Gamma5 problem(thanks to Peter Manning)

337. By Olivier Mattelaer on 2012-01-10

fix the widths problem.
change the syntax to remove the run_name in advanced case and to allow to not put it for pythia/pgs/delphes
change delphes+pgs question
...

338. By Olivier Mattelaer on 2012-01-12

Check the 3/3bar coherence at the import model stage
use that information to fix the Identity -> T object
improve display interactions with particle filter
fix problem in the ordering of the cross-section
improve display style of html page.

339. By Olivier Mattelaer on 2012-01-19

improve autocompletion (especially the - problem)
improve SGE cluster

340. By Olivier Mattelaer on 2012-01-19

put correct version number/Update note

341. By Olivier Mattelaer on 2012-01-20

allow both -modelname and --modelname for the import command

342. By Olivier Mattelaer on 2012-01-20

improve the creation of the gridpack (change store routine)

343. By Olivier Mattelaer on 2012-01-20

tentive try to have a gridpack version working
improve autocompletion, by activating the correct autocompletion after the ';'

344. By Olivier Mattelaer on 2012-01-21

improve auto-completion (adding spaces if this is the only solutions)

345. By Olivier Mattelaer on 2012-01-23

some progress on gridpack (not fully working yet)
improve auto-completion

346. By Olivier Mattelaer on 2012-01-23

gridpack working

347. By Olivier Mattelaer on 2012-01-23

merge with the web (don't ask me why)

348. By Olivier Mattelaer on 2012-01-29

small improvment for display + for false import

349. By Olivier Mattelaer on 2012-01-31

- ensure that RunWeb is removed in case of gridpack run.
- correct a bug in the tag name automatic assignment in case
  of long chain with multiple call to pythia/pgs.
- remove a bug in the removal of events.lhe.gz when already present.

350. By Olivier Mattelaer on 2012-01-31

adding crossxhtml to gridrun

351. By Olivier Mattelaer on 2012-01-31

add action both for online and standard mode.
(use an alert in standard mode)

352. By Olivier Mattelaer on 2012-01-31

simplify treatment of tag in aloha

353. By Olivier Mattelaer on 2012-02-01

merge with trunk

354. By Olivier Mattelaer on 2012-02-01

merge with Johan modifications

355. By Olivier Mattelaer on 2012-02-01

fix a stupid bug

356. By Olivier Mattelaer on 2012-02-02

modification for cp3wks05 and for condor cluster

357. By Olivier Mattelaer on 2012-02-02

working version both on cp3wks05 (pbs) and ingrid-ui1 (condor)

358. By Johan Alwall on 2012-02-03

Merged with fix_new_gridpack branch: Reworked the gridpacks so that 1) the granularity is set to 1, 2) channels can run down to a single iteration during the gridrun. This removes the problem with huge variance for granularity channels while actually further reducing the gridpack run time.

359. By mattelaer-olivier on 2012-02-03

improve configuration file treatment (autocompletion for set / save options/ ...)

360. By mattelaer-olivier on 2012-02-03

return a better error if the cluster is wrongly define

361. By mattelaer-olivier on 2012-02-03

fix to solve UIUC and ingrid trouble.

362. By mattelaer-olivier on 2012-02-03

fix the cluster (try #2)

363. By mattelaer-olivier on 2012-02-03

final fix for the UIUC problem.

[Olivier Mattelaer (olivier-mattelaer)]

Since this has receive no feedback and since our meeting about the new interface stated that this one will be based on 1.4.x.
I merge this in the 2.0 series

review: Approve

=== removed symlink 'AUTHORS'
=== target was 'madgraph/AUTHORS'
=== modified file 'Template/Cards/run_card.dat'
--- Template/Cards/run_card.dat	2011-06-01 21:40:33 +0000
+++ Template/Cards/run_card.dat	2012-02-04 00:00:25 +0000
@@ -20,7 +20,7 @@
 #*********************************************************************
 # Tag name for the run (one word)                                    *
 #*********************************************************************
-  'fermi'     = run_tag ! name of the run 
+  tag_1     = run_tag ! name of the run 
 #*********************************************************************
 # Run to generate the grid pack                                      *
 #*********************************************************************
@@ -29,7 +29,7 @@
 # Number of events and rnd seed                                      *
 # Warning: Do not generate more than 100K event in a single run      *
 #*********************************************************************
-  10000       = nevents ! Number of unweighted events requested 
+  10000 = nevents ! Number of unweighted events requested 
       0       = iseed   ! rnd seed (0=assigned automatically=default))
 #*********************************************************************
 # Collider type and energy                                           *
@@ -69,25 +69,25 @@
 # Automatic ptj and mjj cuts if xqcut > 0
 # (turn off for VBF and single top processes)
 #**********************************************************
-   T  = auto_ptj_mjj
+   T  = auto_ptj_mjj  ! Automatic setting of ptj and mjj
 #**********************************************************
 #                                                                    
 #**********************************
 # BW cutoff (M+/-bwcutoff*Gamma)
 #**********************************
-  15  = bwcutoff
+  15  = bwcutoff      ! (M+/-bwcutoff*Gamma)
 #**********************************************************
 # Apply pt/E/eta/dr/mij cuts on decay products or not
-# (note that etmiss/ptll/ptheavy/sorted cuts always apply)
+# (note that etmiss/ptll/ptheavy/ht/sorted cuts always apply)
 #**********************************************************
-   T  = cut_decays
+   T  = cut_decays    ! Cut decay products 
 #*************************************************************
 # Number of helicities to sum per event (0 = all helicities)
 # 0 gives more stable result, but longer run time (needed for
 # long decay chains e.g.).
 # Use >=2 if most helicities contribute, e.g. pure QCD.
 #*************************************************************
-   0  = nhel
+   0  = nhel          ! Number of helicities used per event
 #*******************                                                 
 # Standard Cuts
 #*******************                                                 
@@ -215,7 +215,7 @@
 # maximal pdg code for quark to be considered as a light jet         *
 # (otherwise b cuts are applied)                                     *
 #*********************************************************************
- 4 = maxjetflavor
+ 4 = maxjetflavor    ! Maximum jet pdg code
 #*********************************************************************
 # Jet measure cuts                                                   *
 #*********************************************************************

=== added file 'Template/HTML/crossx.html'
=== added file 'Template/HTML/mgstyle.css'
--- Template/HTML/mgstyle.css	1970-01-01 00:00:00 +0000
+++ Template/HTML/mgstyle.css	2012-02-04 00:00:25 +0000
@@ -0,0 +1,26 @@
+th a:link{
+    color:black;
+}
+
+table {
+             border-width: 1px;
+             border-spacing: 0px;
+             border-style: outset;
+             border-color: black;
+             border-collapse: collapse;
+             background-color: white;
+}
+table th {
+             border-width: 2px;
+             padding: 1px;
+             border-style: solid;
+             border-color: black;
+             background-color: white;
+}
+
+table td {
+             border-width: 2px;
+             padding: 1px;
+             border-style: solid;
+             border-color: black;
+}

=== added file 'Template/HTML/sortable.js'
--- Template/HTML/sortable.js	1970-01-01 00:00:00 +0000
+++ Template/HTML/sortable.js	2012-02-04 00:00:25 +0000
@@ -0,0 +1,335 @@
+/*
+Table sorting script  by Joost de Valk, check it out at http://www.joostdevalk.nl/code/sortable-table/.
+Based on a script from http://www.kryogenix.org/code/browser/sorttable/.
+Distributed under the MIT license: http://www.kryogenix.org/code/browser/licence.html .
+Modify by the MG team for compatibility issue
+
+Copyright (c) 1997-2007 Stuart Langridge, Joost de Valk.
+
+Version 1.5.7
+*/
+
+/* You can change these values */ 
+var europeandate = true;
+var alternate_row_colors = true;
+
+/* Don't change anything below this unless you know what you're doing */
+addEvent(window, "load", sortables_init);
+
+var SORT_COLUMN_INDEX;
+var thead = false;
+
+function sortables_init() {
+	// Find all tables with class sortable and make them sortable
+	if (!document.getElementsByTagName) return;
+	tbls = document.getElementsByTagName("table");
+	for (ti=0;ti<tbls.length;ti++) {
+		thisTbl = tbls[ti];
+		if (((' '+thisTbl.className+' ').indexOf("sortable") != -1) && (thisTbl.id)) {
+			ts_makeSortable(thisTbl);
+			// make it sortable according to the second column
+			if (thisTbl.tHead && thisTbl.tHead.rows.length > 0) {
+			    var firstRow = thisTbl.tHead.rows[thisTbl.tHead.rows.length-1];
+			} else {
+			    var firstRow = thisTbl.rows[0];
+			}
+
+			for (var ci=0;ci<firstRow.cells[1].childNodes.length;ci++) {
+			    if (firstRow.cells[1].childNodes[ci].tagName && firstRow.cells[1].childNodes[ci].tagName.toLowerCase() == 'a') var lnk = firstRow.cells[1].childNodes[ci];
+			}
+                        ts_resortTable(lnk, 1); //order by cross-section
+                        ts_resortTable(lnk, 1); //order by cross-section
+		}
+	}
+}
+
+function ts_makeSortable(t) {
+	if (t.rows && t.rows.length > 0) {
+		if (t.tHead && t.tHead.rows.length > 0) {
+			var firstRow = t.tHead.rows[t.tHead.rows.length-1];
+			thead = true;
+		} else {
+			var firstRow = t.rows[0];
+		}
+	}
+	if (!firstRow) return;
+	
+	// We have a first row: assume it's the header, and make its contents clickable links
+	for (var i=0;i<firstRow.cells.length;i++) {
+		var cell = firstRow.cells[i];
+		var txt = ts_getInnerText(cell);
+		if (cell.className != "unsortable" && cell.className.indexOf("unsortable") == -1) {
+			cell.innerHTML = '<a href="#" class="sortheader" onclick="ts_resortTable(this, '+i+');return false;">'+txt+'<span class="sortarrow"></span></a>';
+		}
+	}
+	if (alternate_row_colors) {
+		alternate(t);
+	}
+}
+
+function ts_getInnerText(el) {
+	if (typeof el == "string") return el;
+	if (typeof el == "undefined") { return el };
+	if (el.innerText) return el.innerText;	//Not needed but it is faster
+	var str = "";
+	
+	var cs = el.childNodes;
+	var l = cs.length;
+	for (var i = 0; i < l; i++) {
+		switch (cs[i].nodeType) {
+			case 1: //ELEMENT_NODE
+				str += ts_getInnerText(cs[i]);
+				break;
+			case 3:	//TEXT_NODE
+				str += cs[i].nodeValue;
+				break;
+		}
+	}
+	return str;
+}
+
+function ts_resortTable(lnk, clid) {
+	var span;
+	for (var ci=0;ci<lnk.childNodes.length;ci++) {
+		if (lnk.childNodes[ci].tagName && lnk.childNodes[ci].tagName.toLowerCase() == 'span') span = lnk.childNodes[ci];
+	}
+	var spantext = ts_getInnerText(span);
+	var td = lnk.parentNode;
+	var column = clid || td.cellIndex;
+	var t = getParent(td,'TABLE');
+	// Work out a type for the column
+	if (t.rows.length <= 1) return;
+	var itm = "";
+	var i = 1;
+	while (itm == "" && i < t.tBodies[0].rows.length) {
+		var itm = ts_getInnerText(t.tBodies[0].rows[i].cells[column]);
+		itm = trim(itm);
+		if (itm.substr(0,4) == "<!--" || itm.length == 0) {
+			itm = "";
+		}
+		i++;
+	}
+	if (itm == "") return; 
+	sortfn = ts_sort_caseinsensitive;
+	if (itm.match(/^\d\d[\/\.-][a-zA-z][a-zA-Z][a-zA-Z][\/\.-]\d\d\d\d$/)) sortfn = ts_sort_date;
+	if (itm.match(/^\d\d[\/\.-]\d\d[\/\.-]\d\d\d{2}?$/)) sortfn = ts_sort_date;
+	if (itm.match(/^-?[£$€Û¢´]\d/)) sortfn = ts_sort_numeric;
+	if (itm.match(/^-?(\d+[,\.]?)+(E[-+][\d]+)?%?$/)) sortfn = ts_sort_numeric;
+	if (itm.match(/^-?(\d+[,\.]?\d+e[-+][\d]+)?%?$/)) sortfn = ts_sort_numeric;
+	if (itm.match(/^-?(\d+[,\.]?)+([\d]+)?%?$/)) sortfn = ts_sort_numeric;
+	SORT_COLUMN_INDEX = column;
+	var firstRow = new Array();
+	var newRows = new Array();
+	for (k=0;k<t.tBodies.length;k++) {
+		for (i=0;i<t.tBodies[k].rows[0].length;i++) { 
+			firstRow[i] = t.tBodies[k].rows[0][i]; 
+		}
+	}
+	for (k=0;k<t.tBodies.length;k++) {
+		if (!thead) {
+			// Skip the first row
+			for (j=1;j<t.tBodies[k].rows.length;j++) { 
+				newRows[j-1] = t.tBodies[k].rows[j];
+			}
+		} else {
+			// Do NOT skip the first row
+			for (j=0;j<t.tBodies[k].rows.length;j++) { 
+				newRows[j] = t.tBodies[k].rows[j];
+			}
+		}
+	}
+	newRows.sort(sortfn);
+	if (span.getAttribute("sortdir") == 'down') {
+			ARROW = '&nbsp;&nbsp;&darr;';
+			newRows.reverse();
+			span.setAttribute('sortdir','up');
+	} else {
+			ARROW = '&nbsp;&nbsp;&uarr;';
+			span.setAttribute('sortdir','down');
+	} 
+    // We appendChild rows that already exist to the tbody, so it moves them rather than creating new ones
+    // don't do sortbottom rows
+    for (i=0; i<newRows.length; i++) { 
+		if (!newRows[i].className || (newRows[i].className && (newRows[i].className.indexOf('sortbottom') == -1))) {
+			t.tBodies[0].appendChild(newRows[i]);
+		}
+	}
+    // do sortbottom rows only
+    for (i=0; i<newRows.length; i++) {
+		if (newRows[i].className && (newRows[i].className.indexOf('sortbottom') != -1)) 
+			t.tBodies[0].appendChild(newRows[i]);
+	}
+	// Delete any other arrows there may be showing
+	var allspans = document.getElementsByTagName("span");
+	for (var ci=0;ci<allspans.length;ci++) {
+		if (allspans[ci].className == 'sortarrow') {
+			if (getParent(allspans[ci],"table") == getParent(lnk,"table")) { // in the same table as us?
+				allspans[ci].innerHTML = '';
+			}
+		}
+	}		
+	span.innerHTML = ARROW;
+	alternate(t);
+}
+
+function getParent(el, pTagName) {
+	if (el == null) {
+		return null;
+	} else if (el.nodeType == 1 && el.tagName.toLowerCase() == pTagName.toLowerCase()) {
+		return el;
+	} else {
+		return getParent(el.parentNode, pTagName);
+	}
+}
+
+function sort_date(date) {	
+	// y2k notes: two digit years less than 50 are treated as 20XX, greater than 50 are treated as 19XX
+	dt = "00000000";
+	if (date.length == 11) {
+		mtstr = date.substr(3,3);
+		mtstr = mtstr.toLowerCase();
+		switch(mtstr) {
+			case "jan": var mt = "01"; break;
+			case "feb": var mt = "02"; break;
+			case "mar": var mt = "03"; break;
+			case "apr": var mt = "04"; break;
+			case "may": var mt = "05"; break;
+			case "jun": var mt = "06"; break;
+			case "jul": var mt = "07"; break;
+			case "aug": var mt = "08"; break;
+			case "sep": var mt = "09"; break;
+			case "oct": var mt = "10"; break;
+			case "nov": var mt = "11"; break;
+			case "dec": var mt = "12"; break;
+			// default: var mt = "00";
+		}
+		dt = date.substr(7,4)+mt+date.substr(0,2);
+		return dt;
+	} else if (date.length == 10) {
+		if (europeandate == false) {
+			dt = date.substr(6,4)+date.substr(0,2)+date.substr(3,2);
+			return dt;
+		} else {
+			dt = date.substr(6,4)+date.substr(3,2)+date.substr(0,2);
+			return dt;
+		}
+	} else if (date.length == 8) {
+		yr = date.substr(6,2);
+		if (parseInt(yr) < 50) { 
+			yr = '20'+yr; 
+		} else { 
+			yr = '19'+yr; 
+		}
+		if (europeandate == true) {
+			dt = yr+date.substr(3,2)+date.substr(0,2);
+			return dt;
+		} else {
+			dt = yr+date.substr(0,2)+date.substr(3,2);
+			return dt;
+		}
+	}
+	return dt;
+}
+
+function ts_sort_date(a,b) {
+	dt1 = sort_date(ts_getInnerText(a.cells[SORT_COLUMN_INDEX]));
+	dt2 = sort_date(ts_getInnerText(b.cells[SORT_COLUMN_INDEX]));
+	
+	if (dt1==dt2) {
+		return 0;
+	}
+	if (dt1<dt2) { 
+		return -1;
+	}
+	return 1;
+}
+function ts_sort_numeric(a,b) {
+	var aa = ts_getInnerText(a.cells[SORT_COLUMN_INDEX]);
+        aa = clean_num(aa);
+	var bb = ts_getInnerText(b.cells[SORT_COLUMN_INDEX]);
+	bb = clean_num(bb);
+	return compare_numeric(aa,bb);
+}
+function compare_numeric(a,b) {
+	var a = parseFloat(a);
+	a = (isNaN(a) ? 0 : a);
+	var b = parseFloat(b);
+	b = (isNaN(b) ? 0 : b);
+	return a - b;
+}
+function ts_sort_caseinsensitive(a,b) {
+	aa = ts_getInnerText(a.cells[SORT_COLUMN_INDEX]).toLowerCase();
+	bb = ts_getInnerText(b.cells[SORT_COLUMN_INDEX]).toLowerCase();
+	if (aa==bb) {
+		return 0;
+	}
+	if (aa<bb) {
+		return -1;
+	}
+	return 1;
+}
+function ts_sort_default(a,b) {
+	aa = ts_getInnerText(a.cells[SORT_COLUMN_INDEX]);
+	bb = ts_getInnerText(b.cells[SORT_COLUMN_INDEX]);
+	if (aa==bb) {
+		return 0;
+	}
+	if (aa<bb) {
+		return -1;
+	}
+	return 1;
+}
+
+function addEvent(elm, evType, fn, useCapture)
+// addEvent and removeEvent
+// cross-browser event handling for IE5+,	NS6 and Mozilla
+// By Scott Andrew
+{
+	if (elm.addEventListener){
+		elm.addEventListener(evType, fn, useCapture);
+		return true;
+	} else if (elm.attachEvent){
+		var r = elm.attachEvent("on"+evType, fn);
+		return r;
+	} else {
+		alert("Handler could not be removed");
+	}
+}
+function clean_num(str) {
+	str = str.replace(new RegExp(/[^-?0-9.eE]/g),"");
+	return str;
+}
+function trim(s) {
+	return s.replace(/^\s+|\s+$/g, "");
+}
+function alternate(table) {
+	// Take object table and get all it's tbodies.
+	var tableBodies = table.getElementsByTagName("tbody");
+	// Loop through these tbodies
+	for (var i = 0; i < tableBodies.length; i++) {
+		// Take the tbody, and get all it's rows
+		var tableRows = tableBodies[i].getElementsByTagName("tr");
+		// Loop through these rows
+		// Start at 1 because we want to leave the heading row untouched
+		for (var j = 0; j < tableRows.length; j++) {
+			// Check if j is even, and apply classes for both possible results
+			if ( (j % 2) == 0  ) {
+				if ( !(tableRows[j].className.indexOf('odd') == -1) ) {
+					tableRows[j].className = tableRows[j].className.replace('odd', 'even');
+				} else {
+					if ( tableRows[j].className.indexOf('even') == -1 ) {
+						tableRows[j].className += " even";
+					}
+				}
+			} else {
+				if ( !(tableRows[j].className.indexOf('even') == -1) ) {
+					tableRows[j].className = tableRows[j].className.replace('even', 'odd');
+				} else {
+					if ( tableRows[j].className.indexOf('odd') == -1 ) {
+						tableRows[j].className += " odd";
+					}
+				}
+			} 
+		}
+	}
+}

=== modified file 'Template/Source/PDF/pdg2pdf.f'
--- Template/Source/PDF/pdg2pdf.f	2011-03-02 20:10:29 +0000
+++ Template/Source/PDF/pdg2pdf.f	2012-02-04 00:00:25 +0000
@@ -27,6 +27,12 @@
       if(iabs(ipart).eq.22) ipart=7
       iporg=ipart
 
+c     This will be called for any PDG code, but we only support up to 7
+      if(iabs(ipart).gt.7)then
+         pdg2pdf=0d0
+         return 
+      endif
+
       ireuse = 0
       do i=1,2
 c     Check if result can be reused since any of last two calls

=== modified file 'Template/Source/PDF/pdg2pdf_lhapdf.f'
--- Template/Source/PDF/pdg2pdf_lhapdf.f	2010-10-30 03:26:37 +0000
+++ Template/Source/PDF/pdg2pdf_lhapdf.f	2012-02-04 00:00:25 +0000
@@ -20,6 +20,13 @@
 
       ipart=ipdg
       if(ipart.eq.21) ipart=0
+      if(iabs(ipart).eq.22) ipart=7
+
+c     This will be called for any PDG code, but we only support up to 7
+      if(iabs(ipart).gt.7)then
+         pdg2pdf=0d0
+         return 
+      endif
 
       if(ih.eq.ihlast.and.x.eq.xlast.and.xmu.eq.xmulast)then
          pdg2pdf=pdflast(ipart);

=== renamed file 'Template/Events/banner_header.txt' => 'Template/Source/banner_header.txt'
=== modified file 'Template/Source/dsample.f'
--- Template/Source/dsample.f	2011-08-19 00:06:00 +0000
+++ Template/Source/dsample.f	2012-02-04 00:00:25 +0000
@@ -1,4 +1,4 @@
-      subroutine sample_full(ndim,ncall,itmax,dsig,ninvar,nconfigs)
+      subroutine sample_full(ndim,ncall,itmax,itmin,dsig,ninvar,nconfigs)
 c**************************************************************************
 c     Driver for sample which does complete integration
 c     This is done in double precision, and should be told the
@@ -6,7 +6,8 @@
 c     Arguments:
 c     ndim       Number of dimensions for integral(number or random #'s/point)
 c     ncall      Number of times to evaluate the function/iteration
-c     itmax      Number of iterations
+c     itmax      Max number of iterations
+c     itmin      Min number of iterations
 c     ninvar     Number of invarients to keep grids on (s,t,u, s',t' etc)
 c     nconfigs   Number of different pole configurations 
 c     dsig       Function to be integrated
@@ -16,7 +17,7 @@
 c
 c Arguments
 c
-      integer ndim,ncall,itmax,ninvar,nconfigs
+      integer ndim,ncall,itmax,itmin,ninvar,nconfigs
       external         dsig
       double precision dsig
 c
@@ -24,7 +25,7 @@
 c
       double precision x(maxinvar),wgt,p(4*maxdim/3+14)
       double precision tdem, chi2, dum
-      integer ievent,kevent,nwrite,iter,nun,luntmp
+      integer ievent,kevent,nwrite,iter,nun,luntmp,itsum
       integer jmax,i,j,ipole
       integer itmax_adjust
 c
@@ -51,8 +52,8 @@
       common /to_accuracy/accur
 
       double precision twgt, maxwgt,swgt(maxevents)
-      integer                             lun, nw
-      common/to_unwgt/twgt, maxwgt, swgt, lun, nw
+      integer                             lun, nw, itminx
+      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itminx
 
       integer nzoom
       double precision  tx(1:3,maxinvar)
@@ -100,6 +101,7 @@
       kevent = 0
       nzoom = 0
       xzoomfact = 1d0
+      itminx = itmin
       if (nsteps .lt. 1) nsteps=1
       nwrite = itmax*ncall/nsteps
 c      open(unit=66,file='.sample_warn',status='unknown')
@@ -140,7 +142,7 @@
                wgt=0d0
             endif
             if (nzoom .le. 0) then
-               call sample_put_point(wgt,x(1),iter,ipole) !Store result
+               call sample_put_point(wgt,x(1),iter,ipole,itmin) !Store result
             else
                nzoom = nzoom -1
                ievent=ievent-1
@@ -165,7 +167,9 @@
          i=i+1
       enddo
       cur_it = i
-      i = cur_it - 3
+c     Use the last 3 iterations or cur_it-1 if cur_it-1 >= itmin but < 3
+      itsum = min(max(itmin,cur_it-1),3)
+      i = cur_it - itsum
       if (i .gt. 0) then
       tmean = 0d0
       tsigma = 0d0
@@ -183,13 +187,14 @@
       nun = neventswritten
 
       chi2 = 0d0
-      do i = cur_it-3,cur_it-1
+      do i = cur_it-itsum,cur_it-1
          chi2 = chi2+(xmean(i)-tmean)**2/xsigma(i)**2
       enddo
       chi2 = chi2/2d0   !Since using only last 3, n-1=2
       write(*,'(a)') '-----------------------------------------------------'
       write(*,'(a)') '---------------------------'
-      write(*,'(a,e12.4)') ' Results Last 3 iters:  Integral = ',tmean
+      write(*,'(a,i3,a,e12.4)') ' Results Last ',itsum,
+     $     ' iters:  Integral = ',tmean
       write(*,'(25x,a,e12.4)') 'Std dev = ',tsigma
       write(*,'(17x,a,f12.4)') 'Chi**2 per DoF. =',chi2
       write(*,'(a)') '-----------------------------------------------------'
@@ -197,15 +202,17 @@
 
       if (nun .lt. 0) nun=-nun   !Case when wrote maximun number allowed
       if (chi2 .gt. 1) tsigma=tsigma*sqrt(chi2)
+c     JA 02/2011 Added twgt to results.dat to allow event generation in
+c     first iteration for gridpack runs
       if (icor .eq. 0) then
-         write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,tsigma,0.0,kevent,nw,
-     &     cur_it-1,nun, nun/max(tmean,1d-99)
+         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5)')tmean,tsigma,0.0,
+     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt
       else
-         write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,0.0,tsigma,kevent,nw,
-     &     cur_it-1,nun, nun/max(tmean,1d-99)
+         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5)')tmean,0.0,tsigma,
+     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt
       endif
 c      do i=1,cur_it-1
-      do i=cur_it-3,cur_it-1
+      do i=cur_it-itsum,cur_it-1
          write(66,'(i4,4e15.5)') i,xmean(i),xsigma(i),xeff(i),xwmax(i)
       enddo
       close(66)
@@ -297,7 +304,7 @@
                wgt=0d0
             endif
             if (nzoom .le. 0) then
-               call sample_put_point(wgt,x(1),iter,ipole) !Store result
+               call sample_put_point(wgt,x(1),iter,ipole,itmin) !Store result
             else
                nzoom = nzoom -1
                ievent=ievent-1
@@ -314,7 +321,9 @@
          i=i+1
       enddo
       cur_it = i
-      i = cur_it - 3
+c     Use the last 3 iterations or cur_it-1 if cur_it-1 >= itmin
+      itsum = min(max(itmin,cur_it-1),3)
+      i = cur_it - itsum
       if (i .gt. 0) then
       tmean = 0d0
       tsigma = 0d0
@@ -334,13 +343,13 @@
       nun = neventswritten
 
       chi2 = 0d0
-      do i = cur_it-3,cur_it-1
+      do i = cur_it-itsum,cur_it-1
          chi2 = chi2+(xmean(i)-tmean)**2/xsigma(i)**2
       enddo
       chi2 = chi2/2d0   !Since using only last 3, n-1=2
       write(*,'(a)') '-----------------------------------------------------'
       write(*,'(a)') '---------------------------'
-      write(*,'(a,e12.4)') ' Results Last 3 iters:  Integral = ',tmean
+      write(*,'(a,i3,a,e12.4)') ' Results Last ',itsum,' iters:  Integral = ',tmean
       write(*,'(25x,a,e12.4)') 'Std dev = ',tsigma
       write(*,'(17x,a,f12.4)') 'Chi**2 per DoF. =',chi2
       write(*,'(a)') '-----------------------------------------------------'
@@ -348,15 +357,17 @@
 
       if (nun .lt. 0) nun=-nun   !Case when wrote maximun number allowed
       if (chi2 .gt. 1) tsigma=tsigma*sqrt(chi2)
+c     JA 02/2011 Added twgt to results.dat to allow event generation in
+c     first iteration for gridpack runs
       if (icor .eq. 0) then
-         write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,tsigma,0.0,kevent,nw,
-     &     cur_it-1,nun, nun/max(tmean,1d-99)
+         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5)')tmean,tsigma,0.0,
+     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt
       else
-         write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,0.0,tsigma,kevent,nw,
-     &     cur_it-1,nun, nun/max(tmean,1d-99)
+         write(66,'(3e12.5,2i9,i5,i9,e10.3,e12.5)')tmean,0.0,tsigma,
+     &     kevent, nw, cur_it-1, nun, nun/max(tmean,1d-99), twgt
       endif
 c      do i=1,cur_it-1
-      do i=cur_it-3,cur_it-1
+      do i=cur_it-itsum,cur_it-1
          write(66,'(i4,4e15.5)') i,xmean(i),xsigma(i),xeff(i),xwmax(i)
       enddo
       close(66)
@@ -371,237 +382,6 @@
 
       end
 
-      subroutine sample_fullx(ndim,ncall,itmax,dsig,ninvar,nconfigs)
-c**************************************************************************
-c     Driver for sample which does complete integration
-c     This is done in double precision, and should be told the
-c     number of possible phasespace choices.
-c     Arguments:
-c     ndim       Number of dimensions for integral(number or random #'s/point)
-c     ncall      Number of times to evaluate the function/iteration
-c     itmax      Number of iterations
-c     ninvar     Number of invarients to keep grids on (s,t,u, s',t' etc)
-c     nconfigs   Number of different pole configurations 
-c     dsig       Function to be integrated
-c**************************************************************************
-      implicit none
-      include 'genps.inc'
-c
-c Arguments
-c
-      integer ndim,ncall,itmax,ninvar,nconfigs
-      external         dsig
-      double precision dsig
-c
-c Local
-c
-      double precision x(maxinvar),wgt,p(4*maxdim/3+14)
-      double precision tdem, chi2
-      integer ievent,kevent,nwrite,iter,nun
-      integer jmax,i,j,ipole
-c
-c     External
-c
-      integer  n_unwgted
-      external n_unwgted
-c
-c Global
-c
-      integer                                      nsteps
-      character*40          result_file,where_file
-      common /sample_status/result_file,where_file,nsteps
-      double precision fx
-      common /to_fx/   fx
-
-      integer           mincfig, maxcfig
-      common/to_configs/mincfig, maxcfig
-
-      double precision     xmean(99),xsigma(99),xwmax(99),xeff(99)
-      common/to_iterations/xmean,    xsigma,    xwmax,    xeff
-
-      double precision    accur
-      common /to_accuracy/accur
-
-      double precision twgt, maxwgt,swgt(maxevents)
-      integer                             lun, nw
-      common/to_unwgt/twgt, maxwgt, swgt, lun, nw
-
-      integer nzoom
-      double precision  tx(1:3,maxinvar)
-      common/to_xpoints/tx, nzoom
-
-      double precision xzoomfact
-      common/to_zoom/  xzoomfact
-
-      double precision tmean, tsigma
-      integer             dim, events, itm, kn, cur_it, invar, configs
-      common /sample_common/
-     .     tmean, tsigma, dim, events, itm, kn, cur_it, invar, configs
-
-      integer              icor
-      common/to_correlated/icor
-c
-c     External
-c
-      logical pass_point
-c
-c     Data
-c
-c      data result_file,where_file,nsteps/'SAMPLE','WHERE.AMI',100/
-c      data accur/-1d0/
-c      data mincfig /1/
-c      data maxcfig /1/
-c      data twgt/-1d0/              !Dont write out events
-c      data lun/27/                 !Unit number for events
-c      data maxwgt/0d0/
-c      data nw/0/                   !Number of events written
-      
-
-c-----
-c Begin Code
-c-----
-      ievent = 0
-      kevent = 0
-      nzoom = 0
-      xzoomfact = 1d0
-      if (nsteps .lt. 1) nsteps=1
-      nwrite = itmax*ncall/nsteps
-c      open(unit=66,file='.sample_warn',status='unknown')
-c      write(66,*) 'Warnings from sample run.',itmax,ncall
-c      close(66)
-      call sample_init(ndim,ncall,itmax,ninvar,nconfigs)
-      call graph_init
-      do i=1,itmax
-         xmean(i)=0d0
-         xsigma(i)=0d0
-      enddo
-c      mincfig=1
-c      maxcfig=nconfigs
-      wgt = 0d0
-c
-c     Main Integration Loop
-c
-      iter = 1
-      do while(iter .le. itmax)
-c
-c     Get integration point
-c
-         call sample_get_config(wgt,iter,ipole)
-         if (iter .le. itmax) then
-            ievent=ievent+1
-            call x_to_f_arg(ndim,ipole,mincfig,maxcfig,ninvar,wgt,x,p)
-            if (pass_point(p)) then
-               xzoomfact = 1d0
-               fx = dsig(p,wgt,0) !Evaluate function
-               if (xzoomfact .gt. 0d0) then
-                  wgt = wgt*fx*xzoomfact
-               else
-                  wgt = -xzoomfact
-               endif
-               if (wgt .gt. 0d0) call graph_point(p,wgt) !Update graphs
-            else
-               fx =0d0
-               wgt=0d0
-            endif
-            if (nzoom .le. 0) then
-               call sample_put_point(wgt,x(1),iter,ipole) !Store result
-            else
-               nzoom = nzoom -1
-               ievent=ievent-1
-            endif
-         endif
-         if (wgt .gt. 0d0) kevent=kevent+1    
-c
-c     Write out progress/histograms
-c
-         if (kevent .ge. nwrite) then
-            nwrite = nwrite+ncall*itmax/nsteps
-            nwrite = min(nwrite,ncall*itmax)
-c            open(unit=22,file=where_file,status='old',
-c     &           access='append',err=99)
-c            write(22,'(2i15)') ievent,kevent
-c            close(22)
-            call graph_store
-         endif
- 99   enddo
-c
-c     All done
-c
-c      open(unit=66,file='.sample_warn',status='old',access='append')      
-c      write(66,*) 'Finished sample',ievent,kevent,wgt
-c      if (wgt .ne. -1) then
-c         do j=1,5
-c            jmax = min(ndim,4*j)
-c            write(66,'(4e19.12)') (x(i),i=(j-1)*4+1,jmax)
-c         enddo
-c      endif
-c      close(66)
-
-      open(unit=66,file='results.dat',status='unknown')
-      i=1
-      do while(xmean(i) .ne. 0 .and. i .lt. cur_it)
-         i=i+1
-      enddo
-      cur_it = i
-      i = cur_it - 3
-      if (i .gt. 0) then
-      tmean = 0d0
-      tsigma = 0d0
-      tdem = 0d0
-      do while (xmean(i) .ne. 0 .and. i .lt. cur_it)
-         tmean = tmean+xmean(i)*xmean(i)**2/xsigma(i)**2
-         tdem = tdem+xmean(i)**2/xsigma(i)**2
-         tsigma = tsigma + xmean(i)**2/ xsigma(i)**2
-         i=i+1
-      enddo
-c      tmean = tmean/dble(i-1)
-c      tsigma= sqrt(tsigma)/dble(i-1)
-      tmean = tmean/tsigma
-c      tsigma= sqrt(tsigma)/dble(3)
-      tsigma= tmean/sqrt(tsigma)
-      nun = n_unwgted()
-
-      chi2 = 0d0
-      do i = cur_it-3,cur_it-1
-         chi2 = chi2+(xmean(i)-tmean)**2/xsigma(i)**2
-      enddo
-      chi2 = chi2/2d0   !Since using only last 3, n-1=2
-c      tsigma = tsigma*sqrt(chi2)
-c      write(*,*) "chi2 / dof=", chi2
-      write(*,'(a)') '-----------------------------------------------------'
-      write(*,'(a)') '---------------------------'
-      write(*,'(a,e12.4)') ' Results Last 3 iters:  Integral = ',tmean
-      write(*,'(25x,a,e12.4)') 'Std dev = ',tsigma
-      write(*,'(17x,a,f12.4)') 'Chi**2 per DoF. =',chi2
-      write(*,'(a)') '-----------------------------------------------------'
-      write(*,'(a)') '---------------------------'
-
-      if (nun .lt. 0) nun=-nun   !Case when wrote maximun number allowed
-      if (chi2 .gt. 1) tsigma=tsigma*sqrt(chi2)
-      if (icor .eq. 0) then
-         write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,tsigma,0.0,kevent,nw,
-     &     cur_it-1,nun, nun/max(tmean,1d-99)
-      else
-         write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,0.0,tsigma,kevent,nw,
-     &     cur_it-1,nun, nun/max(tmean,1d-99)
-      endif
-c      do i=1,cur_it-1
-      do i=cur_it-3,cur_it-1
-         write(66,'(i4,4e15.5)') i,xmean(i),xsigma(i),xeff(i),xwmax(i)
-      enddo
-      close(66)
-      else
-         open(unit=66,file='results.dat',status='unknown')
-         write(66,'(3e12.5,2i9,i5,i9,e10.3)')0,0,0.0,kevent,nw,
-     &     1,0, 0
-         write(66,'(i4,4e15.5)') 1,0,0,0,0
-         close(66)
-
-      endif
-      
-      end
-
-
       subroutine sample_writehtm()
 c***********************************************************************
 c     Writes out results of run in html format
@@ -1336,10 +1116,10 @@
       endif
       end
 
-      subroutine sample_result(mean, sigma)
+      subroutine sample_result(mean, sigma, itmin)
       implicit none
       double precision mean, sigma
-      integer i,cur_it
+      integer i,cur_it,itmin,itsum
       double precision tsigma,tmean,tsig,tdem
 
       double precision     xmean(99),xsigma(99),xwmax(99),xeff(99)
@@ -1351,7 +1131,9 @@
          i=i+1
       enddo
       cur_it = i
-      i = cur_it - 3
+c     Use the last 3 iterations or cur_it-1 if cur_it-1 >= itmin
+      itsum = min(max(itmin,cur_it-1),3)
+      i = cur_it - itsum
       tmean = 0d0
       tsigma = 0d0
       if (i .gt. 0) then
@@ -1393,7 +1175,7 @@
 c
 c     Local
 c
-      integer i, j, k, knt, non_zero, nun
+      integer i, j, k, knt, non_zero, nun,itsum
       double precision vol,xnmin,xnmax,tot,xdum,tmp1,chi2tmp
       double precision rc, dr, xo, xn, x(maxinvar), dum(ng-1)
       save vol,knt
@@ -1447,8 +1229,8 @@
       common /to_error/reliable
 
       double precision twgt, maxwgt,swgt(maxevents)
-      integer                             lun, nw
-      common/to_unwgt/twgt, maxwgt, swgt, lun, nw
+      integer                             lun, nw, itmin
+      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 
 
       real*8             wmax                 !This is redundant
@@ -1889,9 +1671,9 @@
             xdum=dsig(0,0,2)
 c
 c     Add test to see if we have achieved desired accuracy 
-c     JA 8/17/2011 allow minimum 3 iterations instead of 5
+c     Allow minimum itmin iterations
 c
-            if (tsigma .gt. 0d0 .and. cur_it .gt. 3 .and. accur .gt. 0d0) then
+            if (tsigma .gt. 0d0 .and. cur_it .gt. itmin .and. accur .gt. 0d0) then
 
                xmean = tmean/tsigma
                xchi2 = (chi2/xmean/xmean-tsigma)/dble(cur_it-2)               
@@ -1930,7 +1712,7 @@
 c
 c New check to see if we need to keep integrating this one or not.
 c
-            if (cur_it .gt. 3 .and. accur .lt. 0d0) then  !Check luminocity
+            if (cur_it .gt. itmin .and. accur .lt. 0d0) then  !Check luminocity
 c
 c             Lets get the actual number instead 
 c             tjs 5/22/2007
@@ -1941,16 +1723,18 @@
                nun = neventswritten
 c               tmp1 = tmean / tsigma
 c               chi2tmp = (chi2/tmp1/tmp1-tsigma)/dble(cur_it-2)
-c     Calculate chi2 for last three events (ja 03/11)
+c     Calculate chi2 for last few iterations (ja 03/11)
                tmeant = 0d0
                tsigmat = 0d0
-               do i=cur_it-3,cur_it-1
+c     Use the last 3 iterations or cur_it-1 if cur_it-1 >= itmin but < 3
+               itsum = min(max(itmin,cur_it-1),3)
+               do i=cur_it-itsum,cur_it-1
                   tmeant = tmeant+ymean(i)*ymean(i)**2/ysigma(i)**2
                   tsigmat = tsigmat + ymean(i)**2/ ysigma(i)**2
                enddo
                tmeant = tmeant/tsigmat
                chi2tmp = 0d0
-               do i = cur_it-3,cur_it-1
+               do i = cur_it-itsum,cur_it-1
                   chi2tmp = chi2tmp+(ymean(i)-tmeant)**2/ysigma(i)**2
                enddo
                chi2tmp = chi2tmp/2d0  !Since using only last 3, n-1=2

=== modified file 'Template/Source/gen_ximprove.f'
--- Template/Source/gen_ximprove.f	2011-08-19 00:06:00 +0000
+++ Template/Source/gen_ximprove.f	2012-02-04 00:00:25 +0000
@@ -19,8 +19,8 @@
 c
 c     global
 c
-      integer max_np
-      common/max_np/max_np
+      integer max_np,min_iter
+      common/max_np/max_np,min_iter
 
 c
 c     local
@@ -41,7 +41,7 @@
       logical parallel, gen_events
       character*20 param(maxpara),value(maxpara)
       integer npara, nreq, ngran, nhel_refine
-      integer ij, kl, iseed
+      integer ij, kl, iseed, ioffset
       logical Gridpack,gridrun
       logical split_channels
       common /to_split/split_channels
@@ -62,6 +62,7 @@
      &        ', or number events (>1), max processes per job',
      &        ', and whether to split channels (T/F)'
          read(*,*) err_goal, max_np, split_channels
+         min_iter=3
          parallel = .false.
          if (err_goal .lt. 1) then
             write(*,'(a,f8.2,a)') 'Running for accuracy of ',
@@ -79,11 +80,13 @@
       else
          gen_events=.true.
          split_channels=.false.
+c        Allow all the way down to a single iteration for gridruns
+         min_iter=1
          call get_integer(npara,param,value," gevents "  ,nreq  ,2000   )
-         err_goal = 1.2*nreq
+         err_goal = 1.2*nreq ! extra factor to ensure works
          call get_integer(npara,param,value," gseed "  ,iseed  ,4321   )
          call get_integer(npara,param,value," ngran "  ,ngran  , -1)
-         if (ngran.eq.-1) ngran = int(sqrt(real(nreq)))
+         if (ngran.eq.-1) ngran = 1
          write(*,*) "Running on Grid to generate ",nreq," additional events"
          write(*,*) "   with granularity equal to ",ngran
 c
@@ -91,9 +94,12 @@
 c     Modified to allow for more sequences
 c     iseed can be between 0 and 31328*30081
 c     before patern repeats
+c     JA 11/2/2011: Check for ioffset, as in ntuple (ranmar.f)
 c
+         call get_offset(ioffset)
+         iseed = iseed * 31300       
          ij=1802 + mod(iseed,30081)
-         kl=9373 + (iseed/31328)
+         kl=9373 + (iseed/31328) + ioffset
 c         write(*,'($a,i6,a3,i6)') 'Using random seed offsets',jconfig," : ",ioffset
 c         write(*,*) ' with seed', iseed
          do while (ij .gt. 31328)
@@ -193,9 +199,8 @@
      $        i,nevents,gname,nhel_refine)
       else
          open(unit=25,file='../results.dat',status='old',err=199)
+         read(25,*) xtot
          write(*,'(a,e12.3)') 'Reading total xsection ',xtot
-         read(25,*) xtot
-         write(*,'(e12.3)') xtot
  199     close(25)
          if (gridpack) then
             call write_gen_grid(err_goal,dble(ngran),i,nevents,gname,
@@ -227,8 +232,8 @@
 c
 c     global
 c
-      integer max_np
-      common/max_np/max_np
+      integer max_np,min_iter
+      common/max_np/max_np,min_iter
 c
 c     Arguments
 c
@@ -320,7 +325,6 @@
 c
 c     Now write the commands
 c      
-         write(26,20) 'echo $j'
          write(26,20) 'if [[ ! -e $j ]]; then'
          write(26,25) 'mkdir $j'
          write(26,20) 'fi'
@@ -328,7 +332,7 @@
          write(26,20) 'rm -f $k'
 c         write(26,20) 'rm -f moffset.dat'
 
-         write(26,'(5x,a,2i8,a)') 'echo "',npoints,max_iter,
+         write(26,'(5x,a,3i8,a)') 'echo "',npoints,max_iter,min_iter,
      $        '" >& input_sg.txt' 
          write(26,'(5x,a,f8.3,a)') 'echo "',max(elimit/ysec,0.001d0),
      $        '" >> input_sg.txt'
@@ -338,7 +342,7 @@
      &        '"  >> input_sg.txt' !Helicity 
          write(26,'(5x,3a)')'echo "',gn(io(i))(2:ip-1),
      $        '" >>input_sg.txt'
-         write(26,20) 'time ../madevent >> $k <input_sg.txt'
+         write(26,20) '../madevent >> $k <input_sg.txt'
          write(26,20) 'mv ftn26 ftn25'
 c         write(26,20) 'rm ftn26'
          write(26,20) 'cat $k >> log.txt'
@@ -391,16 +395,16 @@
       write(*,*) 'Opening file ',fname
       open (unit=26, file = fname, status='unknown')
       write(26,15) '#!/bin/bash'
-      write(26,15) '#PBS -q ' // PBS_QUE
-      write(26,15) '#PBS -o /dev/null'
-      write(26,15) '#PBS -e /dev/null'
+c      write(26,15) '#PBS -q ' // PBS_QUE
+c      write(26,15) '#PBS -o /dev/null'
+c      write(26,15) '#PBS -e /dev/null'
       write(26,15) 'if [[ "$PBS_O_WORKDIR" != "" ]]; then' 
       write(26,15) '    cd $PBS_O_WORKDIR'
       write(26,15) 'fi'
       write(26,15) 'k=run1_app.log'
       write(lun,15) 'script=' // fname
-      write(lun,15) 'rm -f wait.$script >& /dev/null'
-      write(lun,15) 'touch run.$script'
+c      write(lun,15) 'rm -f wait.$script >& /dev/null'
+c      write(lun,15) 'touch run.$script'
  15   format(a)
       end
 
@@ -431,7 +435,6 @@
 c
 c     Now write the commands
 c      
-c      write(lun,20) 'echo $i'
 c      write(lun,20) 'j=G$i'
 c      write(lun,20) 'if (! -e $j) then'
 c      write(lun,25) 'mkdir $j'
@@ -445,7 +448,7 @@
 c         write(lun,20) 'cp ../../public.sh .'
 c         write(lun,20) 'qsub -N $1$i public.sh >> ../../running_jobs'
 c      else
-c         write(lun,20) 'time ../madevent > $k <input_app.txt'
+c         write(lun,20) '../madevent > $k <input_app.txt'
 c         write(lun,20) 'rm -f ftn25 ftn99'
 c         write(lun,20) 'cp $k log.txt'
 c      endif
@@ -476,8 +479,8 @@
 c
 c     global
 c
-      integer max_np
-      common/max_np/max_np
+      integer max_np,min_iter
+      common/max_np/max_np,min_iter
 c      integer    max_np     !now set in run_config.inc
 c      parameter (max_np = 5)  !number of channels/job
 
@@ -632,7 +635,6 @@
 c
 c     Now write the commands
 c      
-         write(26,20) 'echo $j'
          write(26,20) 'if [[ ! -e $j ]]; then'
          write(26,25) 'mkdir $j'
          write(26,20) 'fi'
@@ -655,7 +657,7 @@
          write(26,20) 'if [[ ! -e ftn25 ]]; then'
 
 
-         write(26,'(9x,a,2i8,a)') 'echo "',npoints,max_iter,
+         write(26,'(9x,a,3i8,a)') 'echo "',npoints,max_iter,min_iter,
      $        '" >& input_sg.txt' 
 c
 c     tjs 8/7/2007-JA 8/17/11 Allow stop when have enough luminocity
@@ -669,7 +671,7 @@
      &        ' " >> input_sg.txt' !Helicity 0=exact
          write(26,'(9x,3a)')'echo "',gn(io(np))(2:ip-1),
      $        '" >>input_sg.txt'
-         write(26,25) 'time ../madevent >> $k <input_sg.txt'
+         write(26,25) '../madevent >> $k <input_sg.txt'
          write(26,25) 'cat $k >> log.txt'
          write(26,25) 'if [[ -e ftn26 ]]; then'
          write(26,25) '     cp ftn26 ftn25'
@@ -678,7 +680,7 @@
 
          write(26,25) 'rm -f $k'
 
-         write(26,'(9x,a,2i8,a)') 'echo "',npoints,max_iter,
+         write(26,'(9x,a,3i8,a)') 'echo "',npoints,max_iter,min_iter,
      $        '" >& input_sg.txt' 
 c
 c tjs 8/7/2007-JA 8/17/11    Change to request luminocity not accuracy
@@ -706,7 +708,7 @@
          write(26,25) 'if [[ -e ftn26 ]]; then'
          write(26,25) '     cp ftn26 ftn25'
          write(26,25) 'fi'
-         write(26,25) 'time ../madevent >> $k <input_sg.txt'
+         write(26,25) '../madevent >> $k <input_sg.txt'
          write(26,25) 'cat $k >> log.txt'
          write(26,20) 'fi'
          write(26,20) 'cd ../'
@@ -741,8 +743,8 @@
 c
 c   global
 c
-      integer max_np
-      common/max_np/max_np
+      integer max_np,min_iter
+      common/max_np/max_np,min_iter
 c
 c     Arguments
 c
@@ -788,7 +790,7 @@
                write(27,*) xtot*ngran/xsec(i)/goal_lum
             endif
             npoints = goal_lum * xsec(i) / xtot
-            if (npoints .lt. min_gevents_wu) npoints = min_gevents_wu
+            if (npoints .lt. ngran) npoints = ngran
             np = np+1
             if (np .gt. max_np) then
                if (fopened) then
@@ -799,15 +801,14 @@
                np = 1
             endif
             ip = index(gn(i),'/')
-            write(*,*) 'Channel ',gn(i)(2:ip-1),
-     $           yerr, jpoints(i),npoints
+            write(*,*) 'Channel ',gn(i)(2:ip-1), goal_lum * xsec(i) / xtot,
+     $           npoints
 
             ip = index(gn(i),'/')
             write(26,'(2a)') 'j=',gn(i)(1:ip-1)
 c
 c           Now write the commands
 c      
-            write(26,20) 'echo $j'
             write(26,20) 'if [[ ! -e $j ]]; then'
             write(26,25) 'mkdir $j'
             write(26,20) 'fi'
@@ -824,21 +825,21 @@
             write(26,20) 'if [[ ! -e ftn25 ]]; then'
 
 
-            write(26,'(9x,a,2i8,a)') 'echo "',npoints,max_iter,
-     $           '" >& input_sg.txt' 
+            write(26,'(9x,a,3i8,a)') 'echo "',max(npoints,min_events),
+     $           max_iter,min_iter,'" >& input_sg.txt' 
 c
 c     tjs 8/7/2007  Allow stop when have enough events
 c
             write(*,*) "Cross section",i,xsec(i),mfact(i)
-            write(26,'(9x,a,f11.3,a)') 'echo "',-npoints*1.2,
-     $        '" >> input_sg.txt'                       !Accuracy
+            write(26,'(9x,a,e13.5,a)') 'echo "',-npoints/xsec(i),
+     $        '" >> input_sg.txt'                       !Luminocity
             write(26,'(9x,a)') 'echo "2" >> input_sg.txt' !Grid Adjustment
             write(26,'(9x,a)') 'echo "1" >> input_sg.txt' !Suppression
             write(26,'(9x,a,i4,a)') 'echo "',nhel_refine,
      &           ' " >> input_sg.txt' !Helicity 0=exact
             write(26,'(9x,3a)')'echo "',gn(i)(2:ip-1),
      $           '" >>input_sg.txt'
-            write(26,25) 'time ../madevent >> $k <input_sg.txt'
+            write(26,25) '../madevent >> $k <input_sg.txt'
             write(26,25) 'cat $k >> log.txt'
             write(26,25) 'if [[ -e ftn26 ]]; then'
             write(26,25) '     cp ftn26 ftn25'
@@ -847,13 +848,13 @@
 
             write(26,25) 'rm -f $k'
             
-            write(26,'(9x,a,2i8,a)') 'echo "',npoints,max_iter,
-     $           '" >& input_sg.txt' 
+            write(26,'(9x,a,3i8,a)') 'echo "',max(npoints,min_events),
+     $           max_iter,min_iter,'" >& input_sg.txt' 
 c
 c tjs 8/7/2007    Change to request events not accuracy
 c
-            write(26,'(9x,a,f11.3,a)') 'echo "',-npoints*1.2,
-     $           '" >> input_sg.txt' !Accuracy
+            write(26,'(9x,a,e13.5,a)') 'echo "',-npoints / xsec(i),
+     $           '" >> input_sg.txt' ! Luminocity
             write(26,'(9x,a)') 'echo "0" >> input_sg.txt'
             write(26,'(9x,a)') 'echo "1" >> input_sg.txt'
 
@@ -866,7 +867,7 @@
             write(26,25) 'if [[ -e ftn26 ]]; then'
             write(26,25) '     cp ftn26 ftn25'
             write(26,25) 'fi'
-            write(26,25) 'time ../madevent >> $k <input_sg.txt'
+            write(26,25) '../madevent >> $k <input_sg.txt'
             write(26,25) 'cat $k >> log.txt'
             write(26,20) 'fi'
             write(26,20) 'cd ../'

=== modified file 'Template/Source/make_opts'
--- Template/Source/make_opts	2011-07-12 03:57:50 +0000
+++ Template/Source/make_opts	2012-02-04 00:00:25 +0000
@@ -8,7 +8,7 @@
 
 # Set FC unless it's defined by an environment variable
 ifeq ($(origin FC),default)
-  FC=g77
+  FC=gfortran
 endif
 
 # Options: dynamic, lhapdf

=== removed file 'Template/Source/scale_events.f'
--- Template/Source/scale_events.f	2010-10-30 03:26:37 +0000
+++ Template/Source/scale_events.f	1970-01-01 00:00:00 +0000
@@ -1,150 +0,0 @@
-      program scale_events
-c********************************************************************
-c     Takes events from events.dat and scales them according to 
-c     the correct cross section for each diagram
-c********************************************************************
-      implicit none
-c
-c     Constants
-c
-      character*(*) symfile
-      parameter (symfile='symfact.dat')
-      character*(*) scaled_file
-      parameter (scaled_file='scaled.dat')
-c
-c     local
-c
-      character*30 dirname
-      double precision xi
-      integer j,k,ip
-c-----
-c  Begin Code
-c-----
-      open(unit=16,file=scaled_file,status='unknown',err=999)
-      open(unit=35,file=symfile,status='old',err=59)
-      do while (.true.)
-         read(35,*,err=99,end=99) xi,j
-         if (j .gt. 0) then
-            if ( (xi-int(xi+.01)) .lt. 1d-5) then
-               k = int(xi+.01)
-               if (k .lt. 10) then
-                  write(dirname,'(a,i1,a)') 'G',k,'/'
-               else if (k .lt. 100) then
-                  write(dirname,'(a,i2,a)') 'G',k,'/'
-               else if (k .lt. 1000) then
-                  write(dirname,'(a,i3,a)') 'G',k,'/'
-               else if (k .lt. 10000) then
-                  write(dirname,'(a,i4,a)') 'G',k,'/'
-               endif
-            else               !Handle B.W.
-               if (xi .lt. 10) then
-                  write(dirname,'(a,f5.3,a,a)') 'G',xi,'/'
-               else if (xi .lt. 100) then
-                  write(dirname,'(a,f6.3,a,a)') 'G',xi,'/'
-               else if (xi .lt. 1000) then
-                  write(dirname,'(a,f7.3,a,a)') 'G',xi,'/'
-               else if (xi .lt. 10000) then
-                  write(dirname,'(a,3f8.3,a,a)') 'G',xi,'/'
-               endif
-            endif
-            ip = index(dirname,'/')
-            write(*,*) 'Scaling ',dirname(1:ip)
-            call scale_dir(dirname,j,ip)
-         endif
-      enddo
- 99   close(35)
-      close(16)
-      stop
-c
-c     Come here if there isn't a symfact file. Means we will work on
-c     this file alone
-c
- 59   dirname="./"
-      j = 1
-      ip = 2
-      write(*,*) 'Scaling ',dirname(1:ip)
-      call scale_dir(dirname,j,ip)
-      close(16)
-      stop   
- 999  write(*,*) 'Error opening file ',scaled_file
-      close(16)
-      end
-
-      subroutine scale_dir(dirname,mfact, ip)
-c********************************************************************
-c     Takes events from events.dat and scales them according to 
-c     the correct cross section. 
-c********************************************************************
-      implicit none
-c
-c     parameters
-c     
-      include 'nexternal.inc'
-      character*(*) event_file,           xsec_file
-      parameter (event_file='events.lhe', xsec_file='results.dat')
-      integer maxexternal
-      parameter (maxexternal=15)
-
-c
-c     Arguments
-c
-      integer mfact, ip
-      character*(30) dirname
-c
-c     Local
-c
-      double precision xsec, sum, wgt, mxwgt
-      double precision x1,x2,p(0:4,nexternal)
-      integer i,j,k, kevent,m, ic(7,maxexternal),n
-      double precision scale,aqed,aqcd
-      integer ievent
-      character*79 buff
-      logical done
-c-----
-c  Begin Code
-c-----     
-      open(unit=15,file=dirname(1:ip) // xsec_file,status='old',err=998)
-      read(15,*) xsec
-      xsec = xsec * mfact 
-      close(15)
-      sum=0d0
-      mxwgt=-1d0
-      kevent = 0
-      open(unit=15,file=dirname(1:ip)//event_file,status='old',err=999)
-      done = .false.
-      do while (.not. done)
-         call read_event(15,p,wgt,n,ic,ievent,scale,aqcd,aqed,done)
-         if (.not. done) then
-            sum=sum+wgt
-            mxwgt = max(wgt,mxwgt)
-            kevent = kevent+1
-         endif
-      enddo      
- 99   close(15)
-      write(*,*) 'Found ',kevent,' events'
-      write(*,*) 'total weight',sum
-      write(*,*) 'Integrated weight',xsec
-c      stop
-c
-c     Now write out scaled events
-c
-      call write_comments(16)
-      open(unit=15,file=dirname(1:ip) //event_file,status='old',err=999)
-      done=.false.
-      m = 0
-      do while (.not. done)
-         m=m+1
-         call read_event(15,p,wgt,n,ic,ievent,scale,aqcd,aqed,done)
-c         call read_event(15,p,wgt,n,ic,done)
-         if (.not. done) then
-            call write_event(16,p,wgt*xsec/sum,n,ic,m,scale,aqcd,aqed)
-         endif
-      enddo 
- 900  close(15)
-      return
- 55   format(i3,4e19.12)         
- 998  write(*,*) 'Error opening file ',dirname(1:ip) // xsec_file
-      return
- 999  write(*,*) 'Error opening file ',dirname(1:ip) //event_file
-      return
-      end

=== removed file 'Template/Source/select_events.f'
--- Template/Source/select_events.f	2010-10-30 03:26:37 +0000
+++ Template/Source/select_events.f	1970-01-01 00:00:00 +0000
@@ -1,300 +0,0 @@
-      program pick_events
-c********************************************************************
-c     Takes events from scaled.dat and plots distributions 
-c********************************************************************
-      implicit none
-      integer nevents
-      double precision mxwgt,mxwgt1, wgt_goal, eff,xsec
-      integer gevents
-c-----
-c   Begin Code
-c-----
-      call read_all_events(nevents,mxwgt,mxwgt1,eff,xsec)
-      write(*,*) 'Number of events',nevents,mxwgt,mxwgt1
-      write(*,'($a,i7,a)') 'Enter number of events desired ( <',
-     &     int(nevents*eff),') '
-      read(*,*) gevents
-      wgt_goal = mxwgt1*Max((nevents*eff)/gevents,1d0)
-c      wgt_goal =  mxwgt1      
-      call write_events(wgt_goal,xsec,gevents)
-      end
-
-      subroutine write_events(wgt_goal,xsec,gevents)
-c********************************************************************
-c     
-c********************************************************************
-      implicit none
-c
-c     parameters
-c     
-      character*(*) scaled_file
-      parameter (scaled_file='events.lhe')
-      character*(*) unwgt_file
-      parameter (unwgt_file='unweighted_events.lhe')
-      integer    max_write
-      parameter (max_write=100000)
-      integer    maxexternal
-      parameter (maxexternal=15)
-c
-c     Arguments
-c
-      double precision wgt_goal,xsec
-      integer gevents
-c
-c     Local
-c
-      double precision sum, wgt, s_over
-      double precision p(0:4,maxexternal),r
-      integer i,j,k,m, kevent, ic(7,maxexternal)
-      integer iseed, n_over, n
-      integer i4,r8,record_length,irc,ng,iol
-      real xran1
-      character*79 buff
-      logical lwrite
-      logical*1 l_store(max_write)
-      logical done
-      double precision scale, aqcd, aqed
-      integer ievent
-
-      external xran1
-c-----
-c  Begin Code
-c-----     
-      sum=0d0
-      s_over = 0d0
-      n_over = 0
-      iseed = 0
-c      mxwgt=-1d0
-      kevent = 0
-      rewind(37)        !This is comment block
-      open(unit=15,file=scaled_file, status='old',err=999)
-      I4 = 4
-      R8 = 8
-      record_length = 3*I4+maxexternal*I4*7+maxexternal*4*R8+3*R8
-      open(unit=16,access='direct',status='scratch',err=999,
-     &     recl=record_length)
-      done=.false.
-      do while (.not. done)
-c         write(*,*) 'Reading Event'
-         call read_event(15,P,wgt,n,ic,ievent,scale,aqcd,aqed,done)
-c         write(*,*) 'read event',n,wgt,wgt_goal
-         if (.not. done) then
-            r = xran1(iseed)
-            lwrite = (wgt / r .gt. wgt_goal)
-c            write(*,*) 'lwrite',lwrite,r
-            if (lwrite) then
-               sum=sum+wgt_goal
-               kevent = kevent+1
-               l_store(kevent)=.false.
-c               write(*,*) 'writing event',kevent
-               write(16,rec=kevent) wgt_goal,n,
-     &              ((ic(i,j),j=1,maxexternal),i=1,7),
-     &              ((p(i,j),i=0,3),j=1,maxexternal),scale,aqcd,aqed 
-               if (wgt .gt. wgt_goal) then
-                  s_over = s_over-wgt_goal+wgt
-                  n_over = n_over + 1
-               endif
-            endif
-         endif
-      enddo
- 99   close(15)
-c      rewind(16)
-c      write(*,*) 'Done',kevent
-      open(unit=26,file=unwgt_file, status='unknown',err=999)
-      call write_comments(26)
-      write(26,'(a2,75a1)') '##',('*', i=1,68)
-      write(26,'(a)') '##'
-      write(26,'(a)') '##-----------------------'
-      write(26,'(a)') '## Unweighting Statistics'
-      write(26,'(a)') '##-----------------------'
-      write(26,'(a)') '##'
-      write(26,'(a,i7)') '##   Number of Events Written   :' , kevent
-      write(26,'(a,i7)') '##   Number of Events Truncated :' , n_over
-      write(26,'(a,f7.2,a)') '##   Truncated Cross Section    :' ,
-     $     s_over*100./(xsec),"%"
-
-      write(26,'(a2,75a1)') '##',('*', i=1,68)
-
-      k=0
-      do while(k .lt. kevent)
-         r = xran1(iseed)
-         ng=0
-         i=0
-c         write(*,*) r,kevent,k
-         do while (r .gt. real(ng)/(kevent-k) .and. i .lt. kevent)
-c            write(*,*) r,kevent,k
-            i=i+1
-            if (.not. l_store(i)) ng=ng+1
-         enddo
-         k=k+1
-         l_store(i)=.true.
-c         write(*,*) 'Writing event ',k,i
-         read(16,rec=i) wgt_goal,n,
-     &        ((ic(i,j),j=1,maxexternal),i=1,7),
-     &        ((p(i,j),i=0,3),j=1,maxexternal),scale,aqcd,aqed 
-
-         wgt_goal = wgt_goal*xsec/sum
-         call write_event(26,p,wgt_goal,n,ic,k,scale,aqcd,aqed)
-      enddo
- 92   close(16,status='delete')
-      close(26)
- 55   format(i3,4e19.11)         
-      write(*,*) 'Wrote ',kevent,' events'
-      write(*,*) 'Integrated weight',xsec
-      write(*,*) 'Truncated events ',n_over,
-     &     s_over, s_over/xsec*100,'%'      
-      return
- 999  write(*,*) 'Error opening file ',scaled_file
-      end
-
-      subroutine read_all_events(kevent, mxwgt,mxwgt1, eff, xsec)
-c********************************************************************
-c********************************************************************
-      implicit none
-c
-c     parameters
-c     
-      character*(*) scaled_file
-      parameter (scaled_file='events.lhe')
-      integer maxexternal
-      parameter (maxexternal=15)
-      include 'run_config.inc'
-      integer    max_read
-      parameter (max_read = 2000000)
-c
-c     Arguments
-c
-      integer kevent
-      double precision mxwgt,mxwgt1,eff,xsec
-c
-c     Local
-c
-      double precision sum, wgt
-      double precision p(0:3,maxexternal)
-      real xwgt(max_read),xtot
-      integer i,j,k,m, ic(7,maxexternal),n
-      double precision scale,aqcd,aqed
-      integer ievent
-      logical done
-      character*79 buff
-c-----
-c  Begin Code
-c-----     
-      sum=0d0
-      mxwgt=-1d0
-      kevent = 0
-      open(unit=15,file=scaled_file, status='old',err=999)
-      done=.false.
-      do while (.not. done)
-         call read_event(15,P,wgt,n,ic,ievent,scale,aqcd,aqed,done)
-c         call read_event(15,p,wgt,n,ic,done)
-         if (.not. done) then
-            sum=sum+wgt
-            kevent = kevent+1
-            xwgt(kevent) = wgt
-            mxwgt = max(wgt,mxwgt)
-         endif
-         if (kevent .ge. max_read) then
-            write(*,*) 'Error too many events to read in select_events'
-     $           , kevent
-            write(*,*) 'Reset max_read in Source/select_events.f'
-            stop
-         endif
-      enddo
- 99   close(15)
- 55   format(i3,4e19.11)         
-      write(*,*) 'Sorting',kevent
-      call sort(kevent,xwgt)
-      xtot = 0d0
-      i = kevent
-      do while (xtot-xwgt(i)*(kevent-i) .lt. sum*trunc_max
-     $     .and. i .gt. 2)      !Find out minimum target
-         xtot = xtot + xwgt(i)
-         i=i-1
-      enddo
-      eff = sum/kevent/xwgt(i)
-      write(*,*) 'Found ',kevent,' events'
-      write(*,*) 'Integrated weight',sum
-      write(*,*) 'Maximum wgt',mxwgt, xwgt(i)
-      write(*,*) 'Average wgt', sum/kevent
-      write(*,*) 'Unweight Efficiency', eff
-      mxwgt1=xwgt(i)
-      xsec = sum
-      return
- 999  write(*,*) 'Error opening file ',scaled_file
-
-      end
-
-
-
-      subroutine sort(n,ra)
-      real ra(n)
-      l=n/2+1
-      ir=n
-10    continue
-        if(l.gt.1)then
-          l=l-1
-          rra=ra(l)
-        else
-          rra=ra(ir)
-          ra(ir)=ra(1)
-          ir=ir-1
-          if(ir.eq.1)then
-            ra(1)=rra
-            return
-          endif
-        endif
-        i=l
-        j=l+l
-20      if(j.le.ir)then
-          if(j.lt.ir)then
-            if(ra(j).lt.ra(j+1))j=j+1
-          endif
-          if(rra.lt.ra(j))then
-            ra(i)=ra(j)
-            i=j
-            j=j+j
-          else
-            j=ir+1
-          endif
-        go to 20
-        endif
-        ra(i)=rra
-      go to 10
-      end
-
-      function xran1(idum)
-      dimension r(97)
-      parameter (m1=259200,ia1=7141,ic1=54773,rm1=3.8580247e-6)
-      parameter (m2=134456,ia2=8121,ic2=28411,rm2=7.4373773e-6)
-      parameter (m3=243000,ia3=4561,ic3=51349)
-      data iff /0/
-      save r, ix1,ix2,ix3
-      if (idum.lt.0.or.iff.eq.0) then
-        iff=1
-        ix1=mod(ic1-idum,m1)
-        ix1=mod(ia1*ix1+ic1,m1)
-        ix2=mod(ix1,m2)
-        ix1=mod(ia1*ix1+ic1,m1)
-        ix3=mod(ix1,m3)
-        do 11 j=1,97
-          ix1=mod(ia1*ix1+ic1,m1)
-          ix2=mod(ia2*ix2+ic2,m2)
-          r(j)=(float(ix1)+float(ix2)*rm2)*rm1
-11      continue
-        idum=1
-      endif
-      ix1=mod(ia1*ix1+ic1,m1)
-      ix2=mod(ia2*ix2+ic2,m2)
-      ix3=mod(ia3*ix3+ic3,m3)
-      j=1+(97*ix3)/m3
-      if(j.gt.97.or.j.lt.1)then
-         write(*,*) 'j is bad in ran1.f',j, 97d0*ix3/m3
-         pause
-      endif
-      xran1=r(j)
-      r(j)=(float(ix1)+float(ix2)*rm2)*rm1
-      return
-      end
-
-

=== modified file 'Template/Source/setrun.f'
--- Template/Source/setrun.f	2011-08-24 23:32:31 +0000
+++ Template/Source/setrun.f	2012-02-04 00:00:25 +0000
@@ -20,6 +20,11 @@
       double precision D
       common/to_dj/D
 c
+c     PARAM_CARD
+c
+      character*30 param_card_name
+      common/to_param_card_name/param_card_name
+c
 c     local
 c     
       integer npara
@@ -372,13 +377,13 @@
 
       if(lp1.ne.0.or.lp2.ne.0) then
           write(*,*) 'A PDF is used, so alpha_s(MZ) is going to be modified'
-          call setpara('param_card.dat',.true.)
+          call setpara(param_card_name)
           asmz=G**2/(16d0*atan(1d0))
           write(*,*) 'Old value of alpha_s from param_card: ',asmz
           call pdfwrap
           write(*,*) 'New value of alpha_s from PDF ',pdlabel,':',asmz
       else
-          call setpara('param_card.dat',.true.)
+          call setpara(param_card_name)
           asmz=G**2/(16d0*atan(1d0))
           nloop=2
           pdlabel='none'

=== modified file 'Template/SubProcesses/addmothers.f'
--- Template/SubProcesses/addmothers.f	2011-09-30 04:08:40 +0000
+++ Template/SubProcesses/addmothers.f	2012-02-04 00:00:25 +0000
@@ -19,18 +19,23 @@
       integer isym(nexternal,99), jsym
       integer i,j,k,ida(2),ns,nres,ires,icl,ito2,idenpart,nc,ic
       integer mo_color,da_color(2),itmp
-      integer ito(-nexternal+3:nexternal),iseed,maxcolor
-      integer icolalt(2,-nexternal+3:2*nexternal-3)
+      integer ito(-nexternal+3:nexternal),iseed,maxcolor,maxorg
+      integer icolalt(2,-nexternal+2:2*nexternal-3)
       double precision qicl(-nexternal+3:2*nexternal-3), factpm
       double precision xtarget
       data iseed/0/
 
+c     Variables for combination of color indices (including multipart. vert)
+      integer maxcolmp
+      parameter(maxcolmp=20)
+      integer ncolmp,icolmp(2,maxcolmp)
+
       double precision ZERO
       parameter (ZERO=0d0)
       double precision pmass(-nexternal:0,lmaxconfigs)
       double precision pwidth(-nexternal:0,lmaxconfigs)
       integer pow(-nexternal:0,lmaxconfigs)
-      logical first_time
+      logical first_time,tchannel
       save pmass,pwidth,pow
       data first_time /.true./
 
@@ -67,9 +72,9 @@
 c      common/to_colstats/ncols,ncolflow,ncolalt,icorg
 
       double precision pt
-      integer get_color
+      integer get_color,elim_indices,set_colmp,fix_tchannel_color
       real ran1
-      external pt,ran1,get_color
+      external pt,ran1,get_color,elim_indices,set_colmp,fix_tchannel_color
 
       if (first_time) then
          include 'props.inc'
@@ -137,6 +142,9 @@
 
       endif ! nc.gt.0
 
+c     Store original maxcolor to know if we have epsilon vertices
+        maxorg=maxcolor
+
 c     
 c     Get mother information from chosen graph
 c     
@@ -144,22 +152,70 @@
 c     First check number of resonant s-channel propagators
         ns=0
         nres=0
-
+        tchannel=.false.
 c     Loop over propagators to find mother-daughter information
-        do i=-1,-nexternal+3,-1
+        do i=-1,-nexternal+2,-1
 c       Daughters
-          ida(1)=iforest(1,i,iconfig)
-          ida(2)=iforest(2,i,iconfig)
-          do j=1,2
-            if(ida(j).gt.0) ida(j)=isym(ida(j),jsym)
-          enddo
+          if(i.gt.-nexternal+2)then
+             ida(1)=iforest(1,i,iconfig)
+             ida(2)=iforest(2,i,iconfig)
+             do j=1,2
+                if(ida(j).gt.0) ida(j)=isym(ida(j),jsym)
+             enddo
+          endif
 c       Decide s- or t-channel
-          if(iabs(sprop(numproc,i,iconfig)).gt.0) then ! s-channel propagator
+          if(i.gt.-nexternal+2.and.
+     $         iabs(sprop(numproc,i,iconfig)).gt.0) then ! s-channel propagator
             jpart(1,i)=sprop(numproc,i,iconfig)
             ns=ns+1
+          else if(nres.gt.0.and.maxcolor.gt.maxorg) then
+c         For t-channel propagators, just check that the colors are ok
+             if(i.eq.-nexternal+2) then
+c            This is the final t-channel, combining with leg 2
+                mo_color=0
+                if(.not.tchannel)then
+c                  There are no previous t-channels, so this is a combination of
+c                  2, 1 and the last s-channel
+                   ida(1)=1
+                   ida(2)=i+1
+                else
+c                  The daughter info is in iforest
+                   ida(1)=iforest(1,i,iconfig)
+                   ida(2)=iforest(2,i,iconfig)
+                endif
+c            Reverse colors of t-channels to get right color ordering
+                ncolmp=0
+                ncolmp=set_colmp(ncolmp,icolmp,2,jpart,
+     $               iforest(1,-max_branch,iconfig),icolalt,
+     $               icolalt(2,2),icolalt(1,2))
+             else
+                jpart(1,i)=tprid(i,iconfig)
+                mo_color=get_color(jpart(1,i))
+                ncolmp=0
+             endif
+             if(mo_color.gt.1.and.
+     $            mo_color.ne.3.and.mo_color.ne.8)then
+                da_color(1)=get_color(jpart(1,ida(1)))
+                da_color(2)=get_color(jpart(1,ida(2)))
+                call write_error(da_color(1), da_color(2), mo_color)
+             endif
+c            Set icolmp for daughters
+             ncolmp=set_colmp(ncolmp,icolmp,ida(2),jpart,
+     $            iforest(1,-max_branch,iconfig),icolalt,
+     $            icolalt(1,ida(2)),icolalt(2,ida(2)))
+c            Reverse colors of t-channels to get right color ordering
+             ncolmp=set_colmp(ncolmp,icolmp,ida(1),jpart,
+     $            iforest(1,-max_branch,iconfig),icolalt,
+     $            icolalt(2,ida(1)),icolalt(1,ida(1)))
+c            Fix t-channel color
+c             print *,'t-channel: ',i,ida(1),ida(2),mo_color
+c             print *,'colors: ',((icolmp(j,k),j=1,2),k=1,ncolmp)
+             maxcolor=fix_tchannel_color(mo_color,maxcolor,
+     $                                   ncolmp,icolmp,i,icolalt)
+             tchannel=.true.
+             cycle
           else
-c         Don't care about t-channel propagators
-            goto 100
+            goto 100             
           endif
 c       Set status codes for propagator
 c          if((igscl(0).ne.0.and.
@@ -186,140 +242,43 @@
 c          endif
 c       Set color info for all s-channels
           mo_color = get_color(jpart(1,i))
-          da_color(1) = get_color(jpart(1,ida(1)))
-          da_color(2) = get_color(jpart(1,ida(2)))
-          if(da_color(2).lt.da_color(1))then
-c            Order daughters according to color
-             itmp=ida(1)
-             ida(1)=ida(2)
-             ida(2)=itmp
-             itmp=da_color(1)
-             da_color(1)=da_color(2)
-             da_color(2)=itmp
-          endif
-c          print *,'graph: ',iconfig
-c          print *,'Resonance: ',i,' daughters ',ida(1),ida(2),
-c     $         ' ids ',jpart(1,i),jpart(1,ida(1)),jpart(1,ida(2)),
-c     $         ' colors ',mo_color,da_color(1),da_color(2)
-          
+c     If inside multipart. vertex (indicated by color 2) cycle
+c         Set tentative mothers
+          jpart(2,i) = 1
+          jpart(3,i) = 2
+c         Set mother info for daughters
+          do j=1,2
+            jpart(2,ida(j)) = i
+            jpart(3,ida(j)) = i
+          enddo
+          if(mo_color.eq.2) cycle
+c     Reset list of color indices
+          ncolmp=0
+c     Add new color indices to list of color indices
+          do j=1,2
+             ncolmp=set_colmp(ncolmp,icolmp,ida(j),jpart,
+     $            iforest(1,-max_branch,iconfig),icolalt,
+     $            icolalt(1,ida(j)),icolalt(2,ida(j)))
+          enddo
+c          print *,'s-channel: ',i,mo_color,ida(1),ida(2)
+c          print *,'colors: ',((icolmp(j,k),j=1,2),k=1,ncolmp)
           if(mo_color.eq.1) then ! color singlet
-             icolalt(1,i) = 0
-             icolalt(2,i) = 0
+             maxcolor=elim_indices(0,0,ncolmp,icolmp,i,icolalt,maxcolor)
           elseif(mo_color.eq.-3) then ! color anti-triplet
-                icolalt(1,i) = 0
-             if(da_color(1).eq.-3.and.da_color(2).eq.1)then
-                icolalt(2,i) = icolalt(2,ida(1))
-             elseif(da_color(1).eq.-3.and.da_color(2).eq.8)then
-                icolalt(2,i) = icolalt(2,ida(2))
-             elseif(da_color(1).eq.-6.and.da_color(2).eq.3)then
-                if(-icolalt(1,ida(1)).eq.icolalt(1,ida(2)))then
-                   icolalt(2,i) = icolalt(2,ida(1))
-                else
-                   icolalt(2,i) = -icolalt(1,ida(1))
-                endif
-             elseif(da_color(1).eq.3.and.da_color(2).eq.3)then
-                maxcolor=maxcolor+1
-                icolalt(2,i) = maxcolor
-             else
-                call write_error(da_color(1), da_color(2), mo_color)
-             endif
+             maxcolor=elim_indices(0,1,ncolmp,icolmp,i,icolalt,maxcolor)
           elseif(mo_color.eq.3) then ! color triplet
-                icolalt(2,i) = 0
-             if(da_color(1).eq.1.and.da_color(2).eq.3)then
-                icolalt(1,i) = icolalt(1,ida(2))
-             elseif(da_color(1).eq.3.and.da_color(2).eq.8)then
-                icolalt(1,i) = icolalt(1,ida(2))
-             elseif(da_color(1).eq.-3.and.da_color(2).eq.6)then
-                if(icolalt(2,ida(1)).eq.icolalt(1,ida(2)))then
-                   icolalt(1,i) = -icolalt(2,ida(2))
-                else
-                   icolalt(1,i) = icolalt(1,ida(2))
-                endif
-             elseif(da_color(1).eq.-3.and.da_color(2).eq.-3)then
-                maxcolor=maxcolor+1
-                icolalt(1,i) = maxcolor
-             else
-                call write_error(da_color(1), da_color(2), mo_color)
-             endif
+             maxcolor=elim_indices(1,0,ncolmp,icolmp,i,icolalt,maxcolor)
           elseif(mo_color.eq.-6) then ! color anti-sextet
-             if(da_color(1).eq.-6.and.da_color(2).eq.1)then
-                icolalt(1,i) = icolalt(1,ida(1))
-                icolalt(2,i) = icolalt(2,ida(1))
-             elseif(da_color(1).eq.-6.and.da_color(2).eq.8)then
-                if(icolalt(2,ida(1)).eq.icolalt(1,ida(2)))then
-                   icolalt(1,i) = icolalt(1,ida(1))
-                   icolalt(2,i) = icolalt(2,ida(2))
-                else
-                   icolalt(1,i) = -icolalt(2,ida(2))
-                   icolalt(2,i) = icolalt(2,ida(1))
-                endif
-             elseif(da_color(1).eq.-3.and.da_color(2).eq.-3)then
-                icolalt(1,i) = -icolalt(2,ida(1))
-                icolalt(2,i) = icolalt(2,ida(2))
-             else
-                call write_error(da_color(1), da_color(2), mo_color)
-             endif
+             maxcolor=elim_indices(0,2,ncolmp,icolmp,i,icolalt,maxcolor)
           elseif(mo_color.eq.6) then ! color sextet
-             if(da_color(1).eq.1.and.da_color(2).eq.6)then
-                icolalt(1,i) = icolalt(1,ida(2))
-                icolalt(2,i) = icolalt(2,ida(2))
-             elseif(da_color(1).eq.6.and.da_color(2).eq.8)then
-                if(icolalt(1,ida(1)).eq.icolalt(2,ida(2)))then
-                   icolalt(1,i) = icolalt(1,ida(2))
-                   icolalt(2,i) = icolalt(2,ida(1))
-                else
-                   icolalt(1,i) = icolalt(1,ida(1))
-                   icolalt(2,i) = -icolalt(1,ida(2))
-                endif
-             elseif(da_color(1).eq.3.and.da_color(2).eq.3)then
-                icolalt(1,i) = icolalt(1,ida(1))
-                icolalt(2,i) = -icolalt(1,ida(2))
-             else
-                call write_error(da_color(1), da_color(2), mo_color)
-             endif
+             maxcolor=elim_indices(2,0,ncolmp,icolmp,i,icolalt,maxcolor)
           elseif(mo_color.eq.8) then ! color octet
-             if(da_color(1).eq.-3.and.da_color(2).eq.3)then
-                icolalt(1,i) = icolalt(1,ida(2))
-                icolalt(2,i) = icolalt(2,ida(1))             
-             elseif(da_color(1).eq.1.and.da_color(2).eq.8)then
-                icolalt(1,i) = icolalt(1,ida(2))
-                icolalt(2,i) = icolalt(2,ida(2))             
-             elseif(da_color(1).eq.8.and.da_color(2).eq.8)then
-                if(icolalt(1,ida(1)).eq.icolalt(2,ida(2)))then
-                   icolalt(1,i) = icolalt(1,ida(2))
-                   icolalt(2,i) = icolalt(2,ida(1))
-                else
-                   icolalt(1,i) = icolalt(1,ida(1))
-                   icolalt(2,i) = icolalt(2,ida(2))
-                endif
-             elseif(da_color(1).eq.-6.and.da_color(2).eq.6)then
-                if(-icolalt(1,ida(1)).eq.icolalt(1,ida(2)))then
-                   icolalt(1,i) = -icolalt(1,ida(2))
-                   icolalt(2,i) = icolalt(2,ida(1))
-                elseif(icolalt(1,ida(1)).eq.icolalt(2,ida(2)))then
-                   icolalt(1,i) = icolalt(1,ida(2))
-                   icolalt(2,i) = icolalt(2,ida(1))
-                elseif(icolalt(2,ida(1)).eq.icolalt(1,ida(2)))then
-                   icolalt(1,i) = -icolalt(2,ida(2))
-                   icolalt(2,i) = -icolalt(1,ida(1))
-                else
-                   icolalt(1,i) = icolalt(1,ida(2))
-                   icolalt(2,i) = -icolalt(1,ida(1))
-                endif
-             else
-                call write_error(da_color(1), da_color(2), mo_color)
-             endif
-          else
+             maxcolor=elim_indices(1,1,ncolmp,icolmp,i,icolalt,maxcolor)
+          else ! 2 indicates multipart. vertex
+             da_color(1) = get_color(jpart(1,ida(1)))
+             da_color(2) = get_color(jpart(1,ida(2)))
              call write_error(da_color(1), da_color(2), mo_color)
           endif
-c         Set tentative mothers
-          jpart(2,i) = 1
-          jpart(3,i) = 2
-c         Set mother info for daughters
-          do j=1,2
-            jpart(2,ida(j)) = i
-            jpart(3,ida(j)) = i
-          enddo
 c       Just zero helicity info for intermediate states
           jpart(7,i) = 0
         enddo                   ! do i
@@ -371,11 +330,25 @@
       return
       end
 
+c     *************************************
       subroutine write_error(ida1,ida2,imo)
+c     *************************************
       implicit none
       integer ida1,ida2,imo
 
       open(unit=26,file='../../../error',status='unknown',err=999)
+      if (ida1.eq.1000)then
+         write(26,*) 'Error: too many particles in multipart. vertex,',
+     $        ' please increase maxcolmp in addmothers.f'
+         write(*,*) 'Error: too many particles in multipart. vertex,',
+     $        ' please increase maxcolmp in addmothers.f'
+         stop
+      endif
+      if (ida1.eq.1001)then
+         write(26,*) 'Error: failed to reduce to color indices: ',ida2,imo
+         write(*,*) 'Error: failed to reduce to color indices: ',ida2,imo
+         stop
+      endif
       write(26,*) 'Error: Color combination ',ida1,ida2,
      $     '->',imo,' not implemented in addmothers.f'
       write(*,*) 'Error: Color combination ',ida1,ida2,
@@ -384,3 +357,467 @@
 
  999  write(*,*) 'error'
       end
+
+c     *********************************************************************
+      function set_colmp(ncolmp,icolmp,npart,jpart,forest,icol,icol1,icol2)
+c     *********************************************************************
+      implicit none
+      integer maxcolmp
+      parameter(maxcolmp=20)
+      include 'nexternal.inc'
+      include 'genps.inc'
+c     Arguments
+      integer set_colmp
+      integer ncolmp,icolmp(2,*),npart,icol1,icol2
+      integer icol(2,-nexternal+2:2*nexternal-3)
+      integer jpart(7,-nexternal+3:2*nexternal-3)
+      integer forest(2,-max_branch:-1)
+c     Local
+      integer da_color(2),itmp,ida(2),icolor,ipart,i,j
+      integer get_color,set1colmp
+
+      set_colmp=ncolmp
+      icolor=get_color(jpart(1,npart))
+      if(icolor.eq.1) return
+      if(icolor.eq.2) then
+c     Multiparticle vertex - need to go through daughters and collect all colors
+         ipart=npart
+        do while(icolor.eq.2)
+          ida(1)=forest(1,ipart)
+          ida(2)=forest(2,ipart)
+          da_color(1)=get_color(jpart(1,ida(1)))
+          da_color(2)=get_color(jpart(1,ida(2)))
+c          print *,'iforest: ',ipart,ida(1),ida(2),da_color(1),da_color(2)
+          if(da_color(1).ne.2.and.da_color(2).lt.da_color(1).or.
+     $         da_color(2).eq.2)then
+c            Order daughters according to color, but always color 2 first
+             itmp=ida(1)
+             ida(1)=ida(2)
+             ida(2)=itmp
+             itmp=da_color(1)
+             da_color(1)=da_color(2)
+             da_color(2)=itmp
+          endif
+          do i=1,2
+             if(da_color(i).ne.1.and.da_color(i).ne.2)then
+                ncolmp=set1colmp(ncolmp,icolmp,icol(1,ida(i)),
+     $               icol(2,ida(i)))
+             endif
+          enddo
+          icolor=da_color(1)
+          ipart=ida(1)
+        enddo
+      else
+         ncolmp=set1colmp(ncolmp,icolmp,icol1,icol2)
+      endif
+      set_colmp=ncolmp
+      return
+      end
+
+c     ******************************************************
+      function set1colmp(ncolmp,icolmp,icol1,icol2)
+c     ******************************************************
+      implicit none
+c     Arguments
+      integer maxcolmp
+      parameter(maxcolmp=20)
+      integer set1colmp
+      integer ncolmp,icolmp(2,*),icol1,icol2
+
+      ncolmp=ncolmp+1
+      icolmp(1,ncolmp)=icol1
+      icolmp(2,ncolmp)=icol2
+c     Avoid color sextet-type negative indices
+      if(icolmp(1,ncolmp).lt.0)then
+         ncolmp=ncolmp+1
+         icolmp(2,ncolmp)=-icolmp(1,ncolmp-1)
+         icolmp(1,ncolmp-1)=0
+         icolmp(1,ncolmp)=0
+      elseif(icolmp(2,ncolmp).lt.0)then
+         ncolmp=ncolmp+1
+         icolmp(1,ncolmp)=-icolmp(2,ncolmp-1)
+         icolmp(2,ncolmp-1)=0
+         icolmp(2,ncolmp)=0
+      endif
+      if(ncolmp.gt.maxcolmp)
+     $     call write_error(1000,ncolmp,maxcolmp)
+      set1colmp=ncolmp
+      return
+      end
+
+c********************************************************************
+      function fix_tchannel_color(mo_color,maxcolor,ncolmp,icolmp,ires,
+     $                            icol)
+c********************************************************************
+c     Successively eliminate identical pairwise color indices from the
+c     icolmp list, until only (max) one triplet and one antitriplet remains
+c
+
+      implicit none
+      include 'nexternal.inc'
+      integer fix_tchannel_color
+      integer mo_color,maxcolor,ncolmp,icolmp(2,*)
+      integer ires,icol(2,-nexternal+2:2*nexternal-3)
+      integer i,j,i3,i3bar,max3,max3bar,min3,min3bar,maxcol,mincol
+
+c     Successively eliminate color indices in pairs until only the wanted
+c     indices remain
+      do i=1,ncolmp
+         do j=1,ncolmp
+            if(icolmp(1,i).ne.0.and.icolmp(1,i).eq.icolmp(2,j)) then
+               icolmp(1,i)=0
+               icolmp(2,j)=0
+            endif
+         enddo
+      enddo
+      
+      i3=0
+      i3bar=0
+      icol(1,ires)=0
+      icol(2,ires)=0
+      min3=1000
+      max3=0
+      min3bar=1000
+      max3bar=0
+      do i=1,ncolmp
+         if(icolmp(1,i).gt.0)then
+            i3=i3+1
+c           color for t-channels needs to be reversed
+            if(i3.eq.1) icol(2,ires)=icolmp(1,i)
+            if(i3.eq.2) icol(1,ires)=-icolmp(1,i)
+            if(icolmp(1,i).gt.max3) max3=icolmp(1,i)
+            if(icolmp(1,i).lt.min3) min3=icolmp(1,i)
+         endif
+         if(icolmp(2,i).gt.0)then
+            i3bar=i3bar+1
+c           color for t-channels needs to be reversed
+            if(i3bar.eq.1) icol(1,ires)=icolmp(2,i)
+            if(i3bar.eq.2) icol(2,ires)=-icolmp(2,i)
+            if(icolmp(2,i).gt.max3bar) max3bar=icolmp(2,i)
+            if(icolmp(2,i).lt.min3bar) min3bar=icolmp(2,i)
+         endif
+      enddo
+
+      if(mo_color.eq.0)then
+         icol(1,ires)=0
+         icol(2,ires)=0
+      endif
+
+      fix_tchannel_color=maxcolor
+      if(mo_color.le.1.and.i3.eq.0.and.i3bar.eq.0) return
+      if(mo_color.eq.3.and.(i3.eq.1.and.i3bar.eq.0
+     $     .or.i3bar.eq.1.and.i3.eq.0)) return
+      if(mo_color.eq.8.and.i3.eq.1.and.i3bar.eq.1) return
+
+c     Make sure that max and min don't come from the same octet
+      call clean_max_min(icolmp,ncolmp,max3,min3,max3bar,min3bar,i3,i3bar)
+
+      if(mo_color.le.1.and.i3-i3bar.eq.2.or.
+     $   mo_color.le.1.and.i3bar-i3.eq.2.or.
+     $   mo_color.le.1.and.i3.eq.1.and.i3bar.eq.1) then
+c     Replace the maximum index with the minimum one everywhere
+         maxcol=max(max3,max3bar)
+         if(maxcol.eq.max3) then
+            mincol=min3bar
+         else
+            mincol=min3
+         endif
+         do i=ires+1,-1
+            do j=1,2
+               if(icol(j,i).eq.maxcol)
+     $              icol(j,i)=mincol
+            enddo
+         enddo
+c         print *,'Replaced ',maxcol,' by ',mincol
+      elseif(mo_color.le.1.and.i3.eq.2.and.i3bar.eq.2) then
+c     Replace the maximum indices with the minimum ones everywhere
+         do i=ires+1,-1
+            do j=1,2
+               if(icol(j,i).eq.max3bar)
+     $              icol(j,i)=min3
+               if(icol(j,i).eq.max3)
+     $              icol(j,i)=min3bar
+            enddo
+         enddo
+c         print *,'Replaced ',max3bar,' by ',min3,' and ',max3,' by ',min3bar
+      elseif(mo_color.le.1.and.mod(i3,3).eq.0.and.mod(i3bar,3).eq.0)then
+c     This is epsilon index - do nothing
+         continue
+      else if(mo_color.eq.3.and.mod(i3-i3bar,3).eq.2) then
+c     This is an epsilon index
+         maxcolor=maxcolor+1
+         icol(1,ires)=maxcolor
+         icol(2,ires)=0
+c         print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)
+      else if(mo_color.eq.3.and.mod(i3bar-i3,3).eq.2) then
+c     This is an epsilon index
+         maxcolor=maxcolor+1
+         icol(1,ires)=0
+         icol(2,ires)=maxcolor
+c         print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)
+      else if(mo_color.eq.3.and.(i3-i3bar.eq.1.or.i3bar-i3.eq.1).or.
+     $        mo_color.eq.8.and.i3.eq.2.and.i3bar.eq.2) then
+c     Replace the maximum index with the minimum one everywhere
+c     (we don't know if we should replace i3 with i3bar or vice versa)
+         maxcol=max(max3,max3bar)
+         if(maxcol.eq.max3) then
+            mincol=min3bar
+         else
+            mincol=min3
+         endif
+         do i=ires+1,-1
+            do j=1,2
+               if(icol(j,i).eq.maxcol)
+     $              icol(j,i)=mincol
+            enddo
+         enddo
+c     Fix the color for ires (remember 3<->3bar for t-channels)
+         icol(1,ires)=0
+         icol(2,ires)=0
+c         print *,'Replaced ',maxcol,' by ',mincol
+         if(max3.eq.maxcol)then
+            if(i3-i3bar.ge.0) icol(2,ires)=min3
+            if(i3bar-i3.ge.0) icol(1,ires)=max3bar
+         else
+            if(i3-i3bar.ge.0) icol(2,ires)=max3
+            if(i3bar-i3.ge.0) icol(1,ires)=min3bar
+         endif
+c         print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)
+      else
+c     Don't know how to deal with this
+         call write_error(i3,i3bar,mo_color)
+      endif
+
+      fix_tchannel_color=maxcolor
+      
+      return
+      end
+
+c*******************************************************************
+      function elim_indices(n3,n3bar,ncolmp,icolmp,ires,icol,maxcolor)
+c*******************************************************************
+c     Successively eliminate identical pairwise color indices from the
+c     icolmp list, until only the wanted indices remain
+c     n3 gives the number of triplet indices, n3bar number of antitriplets
+c     n3=1 for triplet, n3bar=1 for antitriplet, 
+c     (n3,n3bar)=(1,1) for octet,
+c     n3=2 for sextet, n3bar=2 for antisextet 
+c     If there are epsilon^{ijk} or epsilonbar color couplings, we
+c     need to introduce new index based on maxcolor.
+c
+
+      implicit none
+      include 'nexternal.inc'
+      integer elim_indices
+      integer n3,n3bar,ncolmp,icolmp(2,*),maxcolor
+      integer ires,icol(2,-nexternal+2:2*nexternal-3)
+      integer i,j,i3,i3bar
+
+c     Successively eliminate color indices in pairs until only the wanted
+c     indices remain
+      do i=1,ncolmp
+         do j=1,ncolmp
+            if(icolmp(1,i).ne.0.and.icolmp(1,i).eq.icolmp(2,j)) then
+               icolmp(1,i)=0
+               icolmp(2,j)=0
+            endif
+         enddo
+      enddo
+      
+      i3=0
+      i3bar=0
+      icol(1,ires)=0
+      icol(2,ires)=0
+      do i=1,ncolmp
+         if(icolmp(1,i).gt.0)then
+            i3=i3+1
+            if(i3.eq.1) icol(1,ires)=icolmp(1,i)
+            if(i3.eq.2) icol(2,ires)=-icolmp(1,i)
+         endif
+         if(icolmp(2,i).gt.0)then
+            i3bar=i3bar+1
+            if(i3bar.eq.1) icol(2,ires)=icolmp(2,i)
+            if(i3bar.eq.2) icol(1,ires)=-icolmp(2,i)
+         endif
+      enddo
+
+      if(n3.eq.0) icol(1,ires)=0
+      if(n3bar.eq.0) icol(2,ires)=0
+
+      if(i3.ne.n3.or.i3bar.ne.n3bar) then
+         if(n3.gt.0.and.n3bar.eq.0.and.mod(i3bar+n3,3).eq.0)then
+c        This is an epsilon index interaction
+            maxcolor=maxcolor+1
+            icol(1,ires)=maxcolor
+            if(n3.eq.2)then
+               maxcolor=maxcolor+1
+               icol(2,ires)=-maxcolor
+            endif
+         elseif(n3bar.gt.0.and.n3.eq.0.and.mod(i3+n3bar,3).eq.0)then
+c        This is an epsilonbar index interaction
+            maxcolor=maxcolor+1
+            icol(2,ires)=maxcolor
+            if(n3.eq.2)then
+               maxcolor=maxcolor+1
+               icol(1,ires)=-maxcolor
+            endif
+         elseif(n3.gt.0.and.n3bar.eq.0.and.i3-i3bar.eq.n3.or.
+     $          n3bar.gt.0.and.n3.eq.0.and.i3bar-i3.eq.n3bar.or.
+     $          n3.eq.1.and.n3bar.eq.1.and.i3-i3bar.eq.0.or.
+     $          n3.eq.0.and.n3bar.eq.0.and.i3-i3bar.eq.0)then
+c        We have a previous epsilon which gives the wrong pop-up index
+            call fix_s_color_indices(n3,n3bar,i3,i3bar,ncolmp,icolmp,ires,icol)
+         else
+c           Don't know how to deal with this
+            call write_error(1001,n3,n3bar)
+         endif
+      endif
+
+      elim_indices=maxcolor
+      
+      return
+      end
+
+c*******************************************************************
+      subroutine fix_s_color_indices(n3,n3bar,i3,i3bar,ncolmp,icolmp,
+     $                             ires,icol)
+c*******************************************************************
+c
+c     Fix color flow if some particle has got the wrong pop-up color
+c     due to epsilon-ijk vertices
+c
+
+      implicit none
+      include 'nexternal.inc'
+      integer n3,n3bar,ncolmp,icolmp(2,*),maxcolor
+      integer ires,icol(2,-nexternal+2:2*nexternal-3)
+      integer i,j,i3,i3bar
+      integer max_n3,max_n3bar,min_n3,min_n3bar,maxcol,mincol
+
+      max_n3=0
+      max_n3bar=0
+      min_n3=1000
+      min_n3bar=1000
+      do i=1,ncolmp
+         if(icolmp(1,i).gt.max_n3)
+     $        max_n3=icolmp(1,i)
+         if(icolmp(2,i).gt.max_n3bar)
+     $        max_n3bar=icolmp(2,i)
+         if(icolmp(1,i).gt.0.and.icolmp(1,i).lt.min_n3)
+     $        min_n3=icolmp(1,i)
+         if(icolmp(2,i).gt.0.and.icolmp(2,i).lt.min_n3bar)
+     $        min_n3bar=icolmp(2,i)
+      enddo
+
+      icol(1,ires)=0
+      icol(2,ires)=0
+
+c     Make sure that max and min don't come from the same octet
+      call clean_max_min(icolmp,ncolmp,max_n3,min_n3,
+     $                   max_n3bar,min_n3bar,i3,i3bar)
+
+      if(n3.eq.1.and.n3bar.eq.0.and.i3-i3bar.eq.n3.or.
+     $   n3bar.eq.1.and.n3.eq.0.and.i3bar-i3.eq.n3bar.or.
+     $   n3bar.eq.1.and.n3.eq.1.and.i3bar-i3.eq.0.or.
+     $   n3bar.eq.0.and.n3.eq.0.and.i3bar-i3.eq.0)then
+c     Replace the highest 3bar-index with the lowest 3-index,
+c     and vice versa
+         maxcol=max(max_n3,max_n3bar)
+         if(maxcol.eq.max_n3) then
+            mincol=min_n3bar
+         else
+            mincol=min_n3
+         endif
+         do i=ires,-1
+            do j=1,2
+               if(icol(j,i).eq.maxcol)
+     $              icol(j,i)=mincol
+            enddo
+         enddo
+c         print *,'Replaced ',maxcol,' with ',mincol
+         if(max_n3.eq.maxcol)then
+            if(n3.eq.1) icol(1,ires)=min_n3
+            if(n3bar.eq.1) icol(2,ires)=max_n3bar
+         else
+            if(n3.eq.1) icol(1,ires)=max_n3
+            if(n3bar.eq.1) icol(2,ires)=min_n3bar
+         endif
+c         print *,'Set mother color for ',ires,' to ',(icol(j,ires),j=1,2)
+      else
+c     Don't know how to deal with this
+         call write_error(1001,n3,n3bar)
+      endif
+      return
+      end
+
+c*******************************************************************************
+      subroutine clean_max_min(icolmp,ncolmp,max3,min3,max3bar,min3bar,i3,i3bar)
+c*******************************************************************************
+      implicit none
+      integer ncolmp,icolmp(2,*)
+      integer i,j,max3,max3bar,min3,min3bar,i3,i3bar
+      
+c     Make sure that max and min don't come from the same octet
+      do i=1,ncolmp
+         if(icolmp(1,i).eq.max3.and.icolmp(2,i).eq.min3bar)then
+            min3bar=1000
+            do j=1,ncolmp
+               if(j.eq.i) cycle
+               if(icolmp(2,j).lt.min3bar) min3bar=icolmp(2,j)
+            enddo
+         endif
+         if(icolmp(2,i).eq.max3bar.and.icolmp(1,i).eq.min3)then
+            min3=1000
+            do j=1,ncolmp
+               if(j.eq.i) cycle
+               if(icolmp(1,j).lt.min3) min3=icolmp(1,j)
+            enddo
+         endif
+      enddo
+
+c     ...and that max and min are different
+      if(i3.gt.1.and.max3.eq.min3.or.i3bar.gt.1.and.max3bar.eq.min3bar)then
+         if(max3.gt.max3bar)then
+c        Need to change min3, while still ensuring that max3bar is ok
+            max3bar=0
+            min3=1000
+            do i=1,ncolmp
+               if(icolmp(2,i).gt.max3bar.and.icolmp(2,i).ne.min3bar)
+     $              max3bar=icolmp(2,i)
+               if(icolmp(1,i).lt.min3.and.icolmp(1,i).ne.max3)
+     $              min3=icolmp(1,i)
+            enddo
+c     Make sure that max3bar and min3 don't come from the same octet
+            do i=1,ncolmp
+               if(icolmp(2,i).eq.max3bar.and.icolmp(1,i).eq.min3)then
+                  min3=1000
+                  do j=1,ncolmp
+                     if(j.eq.i.or.icolmp(1,j).eq.max3) cycle
+                     if(icolmp(1,j).lt.min3) min3=icolmp(1,j)
+                  enddo
+               endif
+            enddo
+         else
+c        Need to change min3bar, while still ensuring that max3 is ok
+            max3=0
+            min3bar=1000
+            do i=1,ncolmp
+               if(icolmp(1,i).gt.max3.and.icolmp(1,i).ne.min3)
+     $              max3=icolmp(1,i)
+               if(icolmp(2,i).lt.min3bar.and.icolmp(2,i).ne.max3bar)
+     $              min3bar=icolmp(2,i)
+            enddo
+c     Make sure that max3 and min3bar don't come from the same octet
+            do i=1,ncolmp
+               if(icolmp(1,i).eq.max3.and.icolmp(2,i).eq.min3bar)then
+                  min3bar=1000
+                  do j=1,ncolmp
+                     if(j.eq.i.or.icolmp(2,j).eq.max3bar) cycle
+                     if(icolmp(2,j).lt.min3bar) min3bar=icolmp(2,j)
+                  enddo
+               endif
+            enddo
+         endif
+      endif
+
+      end

=== removed file 'Template/SubProcesses/check_dip.f'
--- Template/SubProcesses/check_dip.f	2010-10-30 03:26:37 +0000
+++ Template/SubProcesses/check_dip.f	1970-01-01 00:00:00 +0000
@@ -1,537 +0,0 @@
-      PROGRAM DRIVER
-C**************************************************************************
-C     THIS IS THE DRIVER FOR CHECKING THE STANDALONE MATRIX ELEMENT
-C     INCLUDING SUBTRACTION TERMS IN THE SOFT/COLINEAR LIMITS.
-C     IT USES A SIMPLE PHASE SPACE GENERATOR
-C     Fabio Maltoni - 3rd Febraury 2007
-C     Updated by Nicolas Greiner and Rikkert Frederix for the
-C     dipole terms - July 2008
-C**************************************************************************
-      IMPLICIT NONE
-C     
-C     CONSTANTS  
-C     
-      REAL*8 ZERO
-      PARAMETER (ZERO=0D0)
-C     
-C     INCLUDE FILES
-C     
-C---  the include file with the values of the parameters and masses	
-      INCLUDE "coupl.inc"
-C---  integer nexternal ! number particles (incoming+outgoing) in the me 
-      INCLUDE "nexternal.inc" 
-C---  particle masses
-      REAL*8 PMASS(NEXTERNAL)	
-C---  integer    n_max_cg
-      INCLUDE "ngraphs.inc"     !how many diagrams (could be useful to know...)
-      INCLUDE "dipole.inc"
-C     
-C     LOCAL
-C     
-      INTEGER I,J,K,l,m,l1,m1
-      REAL*8 P(0:3,NEXTERNAL)      ! four momenta. Energy is the zeroth component.
-      REAL*8 SQRTS                 ! sqrt(s)= center of mass energy 
-      REAL*8 ME,SUBTRACT,SUBTRACT2 ! Matrix Element and subtraction terms
-      real*8 s(nexternal,nexternal)! Final state invariants
-      real*8 start,number          ! Starting value for invariants and the number of phase space points
-      logical cut                  ! Cut to prevent going into double limits
-      real*8 sum,x(2)
-C     
-C     EXTERNAL
-C     
-      REAL*8 DOT
-      EXTERNAL DOT
-      
-C-----
-C     BEGIN CODE
-C-----
-C     
-C---  INITIALIZATION CALLS
-C     
-c---  Call to initialize the values of the couplings, masses and widths 
-c     used in the evaluation of the matrix element. The primary parameters of the
-c     models are read from Cards/param_card.dat. The secondary parameters are calculated
-c     in Source/MODEL/couplings.f. The values are stored in common blocks that are listed
-c     in coupl.inc .
-
-      call setpara('param_card.dat',.true.)  !first call to setup the paramaters
-      mu  = 91.188d0  !renormalization scale
-      muf = 91.188d0  !factorization scale
-      include "pmass.inc"                    !set up masses
-
-      SQRTS=1000d0              !CMS energy in GEV
-      x(1)=1d0                  !Bjorken x's
-      x(2)=1d0
-      number=1000000            !Number of phase space points
-
-
-c********************************************************************************************
-C First do the collinear limits:
-c********************************************************************************************
-      write (*,*) 'Check of all the collinear limits:'
-c Loop over all possible limits of final sate 
-      do l=1,nexternal-1
-         do m=l+1,nexternal
-
-            start=100d0
-            write (*,*) ' '
-            write (*,'(2A)')'-----------------------------------------',
-     .          '------------------------------------------------------'
-            write (*,*) ' Limit: p(',l,').p(',m,') goes to zero'
-            write (*,'(2A)')'-----------------------------------------',
-     .          '------------------------------------------------------'
-            write (*,11) 'p(',l,').p(',m,')/s(1,2) ,',' sqrt(s(',l
-     .           ,',',m,')) ,','|M|^2 ,'
-     .           ,'|Sub.term|^2 ,','|M|^2/|Sub.term|^2'
-            write (*,'(2A)')'-----------------------------------------',
-     .          '------------------------------------------------------'
-
-
-c     Loop over phase space points
-            do i=1,number
-
-c---  Now use a simple multipurpose PS generator (RAMBO) just to get a 
-c     RANDOM set of four momenta of given masses pmass(i) to be used to evaluate 
-c     the madgraph matrix-element and the subtraction terms.
-c     	
-               CALL GET_MOMENTA(SQRTS,PMASS,P)	
-
-c     Calculate all the final state invariants
-               do l1=1,nexternal-1
-                  do m1=l1+1,nexternal
-                     s(l1,m1)=dot(p(0,l1),p(0,m1))
-                  enddo
-               enddo
-               
-               cut =.false.
-c     Only continue if we are closer to the limit then the starting value.
-c     If we pass the if statement, update 'start' with the new value.
-               if(s(l,m)/s(1,2) .lt. start) then
-                  
-c     Make sure that we are only looking in a single collinear limit. Double
-c     logarithms are not cancelled by the dipoles. Do this by cutting
-c     away the phase-space point if any of the other final state invariants
-c     is smaller than 10% of the CMS energy.
-                  do j=1,nexternal-1
-                     do k=j+1,nexternal
-                        if((j.ne.l).or.(k.ne.m)) then
-                           if(s(j,k)/s(1,2).lt.0.01d0) then
-                              cut=.true.
-                           endif
-                        endif
-                     enddo
-                  enddo
-
-c     If we found a point closer to the limit, update 'start' and
-c     calculate the matrix element and subtraction terms.
-                  if(.not.cut) then
-                     start=s(l,m)/s(1,2)
-                     CALL SMATRIX(P,ME)
-                     Call DIPOLSUM(P,X,1d0,SUBTRACT)
-c     Uncomment the following 2 lines to include the possible non-divergent dipoles,
-c     i.e. dipoles for which the unresolved is massive (= in particles.dat the mass
-c     is set to something different from 'ZERO').
-c     Call DIPOLSUMFINITE(P,X,1d0,SUBTRACT2)
-c     SUBTRACT=SUBTRACT+SUBTRACT2
-
-c     Write out the results
-                     Write (*,10)s(l,m)/s(1,2),' ,',dsqrt(abs(
-     &                    (p(0,l)+p(0,m))**2-(p(1,l)+p(1,m))**2
-     &                    -(p(2,l)+p(2,m))**2-(p(3,l)+p(3,m))**2))
-     &                    ,' ,',ME,' ,',SUBTRACT,' ,',ME/SUBTRACT
-                  endif
-               endif
-            enddo
-         enddo
-      enddo
-
-
-
-
-
-
-
-c********************************************************************************
-C Now check also the soft limits:
-c********************************************************************************
-      do i=1,5
-         write (*,*) ' '
-      enddo
-      write (*,*) 'Check off all the soft limits:'
-c loop over all final state particles:
-      do l=3, nexternal
-
-            start=100d0
-            write (*,*) ' '
-            write (*,'(2A)')'-----------------------------------------',
-     .          '------------------------------------------------------'
-            write (*,*) ' Limit: p(',l,') goes soft'
-            write (*,'(2A)')'-----------------------------------------',
-     .          '------------------------------------------------------'
-            write (*,13) 'p(0,',l,')^2/s(1,2) ,','|M|^2 ,'
-     .           ,'|Sub.term|^2 ,','|M|^2/|Sub.term|^2'
-            write (*,'(2A)')'-----------------------------------------',
-     .          '------------------------------------------------------'
-
-c     Loop over phase space points
-            do i=1,number
-
-c     Get the random set of 4-momenta
-               CALL GET_MOMENTA(SQRTS,PMASS,P)	
-
-c     Calculate all the final state invariants
-               do l1=1,nexternal-1
-                  do m1=l1+1,nexternal
-                     s(l1,m1)=dot(p(0,l1),p(0,m1))
-                  enddo
-               enddo
-               
-               cut =.false.
-               sum = 0d0
-c     Only continue if we are closer to the limit then the starting value.
-c     If we pass the if statement, update 'start' with the new value.
-               if(p(0,l)**2/s(1,2) .lt. start) then
-                  
-c     Make sure that we are only looking at the soft limit. Do this by cutting
-c     away the phase-space point if any of the other final state invariants
-c     is smaller than 10% of the CMS energy.
-                  do j=1,nexternal-1
-                     do k=j+1,nexternal
-                        if((j.ne.l).and.(k.ne.l)) then
-                           if(s(j,k)/s(1,2).lt.0.01d0) then
-                              cut=.true.
-                           endif
-                        endif
-                        if((j.eq.l).or.(k.eq.l))then
-                           sum=sum+s(j,k)
-                        endif
-                     enddo
-                  enddo
-c     cut away points that are collinear by forcing that all the invariants
-c     with the soft particle contribute at least 2*nexternal'th to the sum. This
-c     makes sure that all these invariants go to zero at the same time.              
-                  do j=1,nexternal-1
-                     do k=j+1,nexternal
-                        if(((j.eq.l).or.(k.eq.l)).and.
-     .                       (s(j,k).lt.sum/real(2*nexternal)))then
-                           cut=.true.
-                        endif
-                     enddo
-                  enddo
-
-c     If we found a point closer to the limit, update 'start' and
-c     calculate the matrix element and subtraction terms.
-                  if(.not.cut) then
-                     start=p(0,l)**2/s(1,2)
-                     CALL SMATRIX(P,ME)
-                     Call DIPOLSUM(P,X,1d0,SUBTRACT)
-c     Uncomment the following 2 lines to include the possible non-divergent dipoles,
-c     i.e. dipoles for which the unresolved is massive (= in particles.dat the mass
-c     is set to something different from 'ZERO').
-c     Call DIPOLSUMFINITE(P,X,1d0,SUBTRACT2)
-c     SUBTRACT=SUBTRACT+SUBTRACT2
-
-c     Write out the results
-                     Write (*,14) p(0,l)/s(1,2),
-     .                    ' ,',ME,' ,',SUBTRACT,' ,',ME/SUBTRACT
-                  endif
-               endif
-            enddo
-
-      enddo
-
-
- 10   format(1X,e18.6,A2,e16.6,A2,e16.6,A2,e16.6,A2,e16.6)
- 11   format(1X,A2,I2,A4,I2,A,1X,A,I2,A1,I2,A,A18,A18,A20)
- 14   format(1X,e18.6,A2,e16.6,A2,e16.6,A2,e16.6)
- 13   format(3X,A4,I2,A12,A18,A18,A20)
-
-      end
-	
-	  
-	  
-	  
-      double precision function dot(p1,p2)
-C****************************************************************************
-C     4-Vector Dot product
-C****************************************************************************
-      implicit none
-      double precision p1(0:3),p2(0:3)
-      dot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
-      end
-
-
-      SUBROUTINE GET_MOMENTA(ENERGY,PMASS,P)
-C---- auxiliary function to change convention between madgraph and rambo
-c---- four momenta. 	  
-      IMPLICIT NONE
-      INCLUDE "nexternal.inc"
-C     ARGUMENTS
-      REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT
-C     LOCAL
-      INTEGER I
-      
-      P(0,1)=energy/2
-      P(1,1)=0d0
-      P(2,1)=0d0
-      P(3,1)=energy/2
-      
-      P(0,2)=energy/2
-      P(1,2)=0d0
-      P(2,2)=0d0
-      P(3,2)=-energy/2
-
-      call rambo(nexternal-2,energy,pmass(3),prambo,WGT)
-      DO I=3, NEXTERNAL
-         P(0,I)=PRAMBO(4,I-2)	
-         P(1,I)=PRAMBO(1,I-2)
-         P(2,I)=PRAMBO(2,I-2)
-         P(3,I)=PRAMBO(3,I-2)	
-      ENDDO
-
-      RETURN
-      END
-
-
-      SUBROUTINE RAMBO(N,ET,XM,P,WT)
-***********************************************************************
-*                       RAMBO                                         *
-*    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)                  *
-*                                                                     *
-*    A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR                *
-*    AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING                 *
-*    THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS                      *
-*    -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC (20-08-90)    *
-*                                                                     *
-*    N  = NUMBER OF PARTICLES                                         *
-*    ET = TOTAL CENTRE-OF-MASS ENERGY                                 *
-*    XM = PARTICLE MASSES ( DIM=NEXTERNAL-2 )                         *
-*    P  = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-2) )                    *
-*    WT = WEIGHT OF THE EVENT                                         *
-***********************************************************************
-      IMPLICIT REAL*8(A-H,O-Z)
-	  INCLUDE "nexternal.inc"
-      DIMENSION XM(NEXTERNAL-2),P(4,NEXTERNAL-2)
-      DIMENSION Q(4,NEXTERNAL-2),Z(NEXTERNAL-2),R(4),
-     .   B(3),P2(NEXTERNAL-2),XM2(NEXTERNAL-2),
-     .   E(NEXTERNAL-2),V(NEXTERNAL-2),IWARN(5)
-      SAVE ACC,ITMAX,IBEGIN,IWARN
-      DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
-*
-* INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
-      IF(IBEGIN.NE.0) GOTO 103
-      IBEGIN=1
-      TWOPI=8.*DATAN(1.D0)
-      PO2LOG=LOG(TWOPI/4.)
-      Z(2)=PO2LOG
-      DO 101 K=3,10
-  101 Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
-      DO 102 K=3,10
-  102 Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
-*
-* CHECK ON THE NUMBER OF PARTICLES
-  103 IF(N.GT.1.AND.N.LT.101) GOTO 104
-      PRINT 1001,N
-      STOP
-*
-* CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
-  104 XMT=0.
-      NM=0
-      DO 105 I=1,N
-      IF(XM(I).NE.0.D0) NM=NM+1
-  105 XMT=XMT+ABS(XM(I))
-      IF(XMT.LE.ET) GOTO 201
-      PRINT 1002,XMT,ET
-      STOP
-*
-* THE PARAMETER VALUES ARE NOW ACCEPTED
-*
-* GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
-  201 DO 202 I=1,N
-         r1=rn(1)
-      C=2.*r1-1.
-      S=SQRT(1.-C*C)
-      F=TWOPI*RN(2)
-      r1=rn(3)
-      r2=rn(4)
-      Q(4,I)=-LOG(r1*r2)
-      Q(3,I)=Q(4,I)*C
-      Q(2,I)=Q(4,I)*S*COS(F)
-  202 Q(1,I)=Q(4,I)*S*SIN(F)
-*
-* CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
-      DO 203 I=1,4
-  203 R(I)=0.
-      DO 204 I=1,N
-      DO 204 K=1,4
-  204 R(K)=R(K)+Q(K,I)
-      RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
-      DO 205 K=1,3
-  205 B(K)=-R(K)/RMAS
-      G=R(4)/RMAS
-      A=1./(1.+G)
-      X=ET/RMAS
-*
-* TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
-      DO 207 I=1,N
-      BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
-      DO 206 K=1,3
-  206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
-  207 P(4,I)=X*(G*Q(4,I)+BQ)
-*
-* CALCULATE WEIGHT AND POSSIBLE WARNINGS
-      WT=PO2LOG
-      IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
-      IF(WT.GE.-180.D0) GOTO 208
-      IF(IWARN(1).LE.5) PRINT 1004,WT
-      IWARN(1)=IWARN(1)+1
-  208 IF(WT.LE. 174.D0) GOTO 209
-      IF(IWARN(2).LE.5) PRINT 1005,WT
-      IWARN(2)=IWARN(2)+1
-*
-* RETURN FOR WEIGHTED MASSLESS MOMENTA
-  209 IF(NM.NE.0) GOTO 210
-* RETURN LOG OF WEIGHT
-      WT=WT
-      RETURN
-*
-* MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
-  210 XMAX=SQRT(1.-(XMT/ET)**2)
-      DO 301 I=1,N
-      XM2(I)=XM(I)**2
-  301 P2(I)=P(4,I)**2
-      ITER=0
-      X=XMAX
-      ACCU=ET*ACC
-  302 F0=-ET
-      G0=0.
-      X2=X*X
-      DO 303 I=1,N
-      E(I)=SQRT(XM2(I)+X2*P2(I))
-      F0=F0+E(I)
-  303 G0=G0+P2(I)/E(I)
-      IF(ABS(F0).LE.ACCU) GOTO 305
-      ITER=ITER+1
-      IF(ITER.LE.ITMAX) GOTO 304
-      PRINT 1006,ITMAX
-      GOTO 305
-  304 X=X-F0/(X*G0)
-      GOTO 302
-  305 DO 307 I=1,N
-      V(I)=X*P(4,I)
-      DO 306 K=1,3
-  306 P(K,I)=X*P(K,I)
-  307 P(4,I)=E(I)
-*
-* CALCULATE THE MASS-EFFECT WEIGHT FACTOR
-      WT2=1.
-      WT3=0.
-      DO 308 I=1,N
-      WT2=WT2*V(I)/E(I)
-  308 WT3=WT3+V(I)**2/E(I)
-      WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
-*
-* RETURN FOR  WEIGHTED MASSIVE MOMENTA
-      WT=WT+WTM
-      IF(WT.GE.-180.D0) GOTO 309
-      IF(IWARN(3).LE.5) PRINT 1004,WT
-      IWARN(3)=IWARN(3)+1
-  309 IF(WT.LE. 174.D0) GOTO 310
-      IF(IWARN(4).LE.5) PRINT 1005,WT
-      IWARN(4)=IWARN(4)+1
-* RETURN LOG OF WEIGHT
-  310 WT=WT
-      RETURN
-*
- 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
- 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',
-     . ' SMALLER THAN TOTAL ENERGY =',D15.6)
- 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
- 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
- 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE',
-     . ' DESIRED ACCURACY =',D15.6)
-      END
-
-      FUNCTION RN(IDUMMY)
-      REAL*8 RN,RAN
-      SAVE INIT
-      DATA INIT /1/
-      IF (INIT.EQ.1) THEN
-        INIT=0
-        CALL RMARIN(1802,9373)
-      END IF
-*
-  10  CALL RANMAR(RAN)
-      IF (RAN.LT.1D-16) GOTO 10
-      RN=RAN
-*
-      END
-
-
-
-      SUBROUTINE RANMAR(RVEC)
-*     -----------------
-* Universal random number generator proposed by Marsaglia and Zaman
-* in report FSU-SCRI-87-50
-* In this version RVEC is a double precision variable.
-      IMPLICIT REAL*8(A-H,O-Z)
-      COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
-      COMMON/ RASET2 / IRANMR,JRANMR
-      SAVE /RASET1/,/RASET2/
-      UNI = RANU(IRANMR) - RANU(JRANMR)
-      IF(UNI .LT. 0D0) UNI = UNI + 1D0
-      RANU(IRANMR) = UNI
-      IRANMR = IRANMR - 1
-      JRANMR = JRANMR - 1
-      IF(IRANMR .EQ. 0) IRANMR = 97
-      IF(JRANMR .EQ. 0) JRANMR = 97
-      RANC = RANC - RANCD
-      IF(RANC .LT. 0D0) RANC = RANC + RANCM
-      UNI = UNI - RANC
-      IF(UNI .LT. 0D0) UNI = UNI + 1D0
-      RVEC = UNI
-      END
- 
-      SUBROUTINE RMARIN(IJ,KL)
-*     -----------------
-* Initializing routine for RANMAR, must be called before generating
-* any pseudorandom numbers with RANMAR. The input values should be in
-* the ranges 0<=ij<=31328 ; 0<=kl<=30081
-      IMPLICIT REAL*8(A-H,O-Z)
-      COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
-      COMMON/ RASET2 / IRANMR,JRANMR
-      SAVE /RASET1/,/RASET2/
-* This shows correspondence between the simplified input seeds IJ, KL
-* and the original Marsaglia-Zaman seeds I,J,K,L.
-* To get the standard values in the Marsaglia-Zaman paper (i=12,j=34
-* k=56,l=78) put ij=1802, kl=9373
-      I = MOD( IJ/177 , 177 ) + 2
-      J = MOD( IJ     , 177 ) + 2
-      K = MOD( KL/169 , 178 ) + 1
-      L = MOD( KL     , 169 )
-      DO 300 II = 1 , 97
-        S =  0D0
-        T = .5D0
-        DO 200 JJ = 1 , 24
-          M = MOD( MOD(I*J,179)*K , 179 )
-          I = J
-          J = K
-          K = M
-          L = MOD( 53*L+1 , 169 )
-          IF(MOD(L*M,64) .GE. 32) S = S + T
-          T = .5D0*T
-  200   CONTINUE
-        RANU(II) = S
-  300 CONTINUE
-      RANC  =   362436D0 / 16777216D0
-      RANCD =  7654321D0 / 16777216D0
-      RANCM = 16777213D0 / 16777216D0
-      IRANMR = 97
-      JRANMR = 33
-      END
-
-
-
-
-
-

=== removed file 'Template/SubProcesses/check_intdip.f'
--- Template/SubProcesses/check_intdip.f	2010-10-30 03:26:37 +0000
+++ Template/SubProcesses/check_intdip.f	1970-01-01 00:00:00 +0000
@@ -1,401 +0,0 @@
-      PROGRAM checkintdip
-C**************************************************************************
-C     This is a program for testing the integrated dipoles. For a couple of
-C     phase space points, chosen by RAMBO, it prints the the singular and
-C     finite terms for the given process.
-C     R.Frederix and N.Greiner - February 2010
-C**************************************************************************
-      IMPLICIT NONE
-C     
-C     CONSTANTS  
-C     
-      REAL*8 ZERO
-      PARAMETER (ZERO=0D0)
-C     
-C     INCLUDE FILES
-C     
-C---  the include file with the values of the parameters and masses	
-      INCLUDE "coupl.inc"
-C---  integer nexternal ! number particles (incoming+outgoing) in the me 
-      INCLUDE "nexternal.inc" 
-C---  particle masses
-      REAL*8 PMASS(NEXTERNAL)	
-C---  integer    n_max_cg
-      INCLUDE "ngraphs.inc"     !how many diagrams (could be useful to know...)
-      INCLUDE "dipole.inc" 
-
-C     
-C     LOCAL
-C     
-      INTEGER I,J,K,l,m,l1,m1, number
-      REAL*8 P(0:3,NEXTERNAL)      ! four momenta. Energy is the zeroth component.
-      REAL*8 SQRTS                 ! sqrt(s)= center of mass energy 
-      REAL*8 ME,SUBTRACT,SUBTRACT2 ! Matrix Element and subtraction terms
-      real*8 start                 ! Starting value for invariants and the number of phase space points
-      logical cut                  ! Cut to prevent going into double limits
-      real*8 sum,x(2)
-      real*8 Z, eps,epssq, finite
-C     
-C     EXTERNAL
-C     
-      REAL*8 DOT
-      EXTERNAL DOT
-      
-C-----
-C     BEGIN CODE
-C-----
-C     
-C---  INITIALIZATION CALLS
-C     
-c---  Call to initialize the values of the couplings, masses and widths 
-c     used in the evaluation of the matrix element. The primary parameters of the
-c     models are read from Cards/param_card.dat. The secondary parameters are calculated
-c     in Source/MODEL/couplings.f. The values are stored in common blocks that are listed
-c     in coupl.inc .
-
-      call setpara('param_card.dat',.true.)  !first call to setup the paramaters
-      include "pmass.inc"                    !set up masses
-      mu  = 91.188d0  !renormalization scale
-      muf = 91.188d0  !factorization scale
-
-      SQRTS=1000d0              !CMS energy in GEV
-      x(1)=1d0                  !Bjorken x's
-      x(2)=1d0
-      number=3            !Number of phase space points
-
-c For this check routine to work properly,
-c the last parton should be a massless one
-      if (pmass(nexternal).ne.0d0) then
-         write (*,*) 'Error in process order:'
-         write (*,*) 'Please, change the order of the particles in the'
-         write (*,*)
-     &        'process such that the last one is a massless particle'
-         stop
-      endif
-
-c     Loop over phase space points
-      do i=1,number
-
-c---  Now use a simple multipurpose PS generator (RAMBO) just to get a 
-c     RANDOM set of four momenta of given masses pmass(i) to be used to evaluate 
-c     the integrated dipoles.
-c     	
-         CALL GET_INT_MOMENTA(SQRTS,PMASS,P,Z)
-
-         CALL INTDIPOLES(P,X,Z,1d0,EPS,EPSSQ,FINITE)
-c Note that for the intdipolesfinite the phase space point from
-c Rambo is, in general not correct
-c         CALL INTDIPOLESFINITE(P,X,Z,1d0,EPS,EPSSQ,FINITE)
-         write (*,*) ' '
-         write (*,*) 'PHASE SPACE POINT',i,':'
-         do j=1, nexternal-1
-            write (*,*)'   ',j, P(0,j), P(1,j), P(2,j), P(3,j)
-         enddo
-         write (*,*) '    Z =', Z
-         write (*,*) ' '
-         write (*,*) '    1/eps^2: ',epssq
-         write (*,*) '    1/eps:   ',eps
-         write (*,*) '    Finite:  ',finite
-         write (*,*) ' '
-         write (*,*) ' '
-      enddo
-
-      end
-
-
-
-      SUBROUTINE GET_INT_MOMENTA(ENERGY,PMASS,P,Z)
-C     Provide four-momenta for integrated dipoles.
-
-      IMPLICIT NONE
-      INCLUDE "nexternal.inc"
-C     ARGUMENTS
-      REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT,Z
-C     LOCAL
-      INTEGER I
-      
-      P(0,1)=energy/2
-      P(1,1)=0d0
-      P(2,1)=0d0
-      P(3,1)=energy/2
-      
-      P(0,2)=energy/2
-      P(1,2)=0d0
-      P(2,2)=0d0
-      P(3,2)=-energy/2
-
-      call rambo_int(nexternal-3,energy,pmass(3),prambo,Z)
-      DO I=3, NEXTERNAL-1
-         P(0,I)=PRAMBO(4,I-2)	
-         P(1,I)=PRAMBO(1,I-2)
-         P(2,I)=PRAMBO(2,I-2)
-         P(3,I)=PRAMBO(3,I-2)	
-      ENDDO
-
-      P(0,nexternal)=0d0
-      P(1,nexternal)=0d0
-      P(2,nexternal)=0d0
-      P(3,nexternal)=0d0
-
-      RETURN
-      END
-
-
-
-
-      double precision function dot(p1,p2)
-C****************************************************************************
-C     4-Vector Dot product
-C****************************************************************************
-      implicit none
-      double precision p1(0:3),p2(0:3)
-      dot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
-      end
-
-
-      SUBROUTINE RAMBO_INT(N,ET,XM,P,Zvar)
-C     Modified RAMBO suited for integrated dipoles
-C     R.Frederix and N.Greiner - February 2010
-***********************************************************************
-*                       RAMBO                                         *
-*    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)                  *
-*                                                                     *
-*    A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR                *
-*    AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING                 *
-*    THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS                      *
-*    -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC (20-08-90)    *
-*                                                                     *
-*    N  = NUMBER OF PARTICLES                                         *
-*    ET = TOTAL CENTRE-OF-MASS ENERGY                                 *
-*    XM = PARTICLE MASSES ( DIM=NEXTERNAL-2 )                         *
-*    P  = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-2) )                    *
-*    WT = WEIGHT OF THE EVENT                                         *
-***********************************************************************
-      IMPLICIT REAL*8(A-H,O-Z)
-	  INCLUDE "nexternal.inc"
-      DIMENSION XM(NEXTERNAL-2),P(4,NEXTERNAL-2)
-      DIMENSION Q(4,NEXTERNAL-2),Z(NEXTERNAL-2),R(4),
-     .   B(3),P2(NEXTERNAL-2),XM2(NEXTERNAL-2),
-     .   E(NEXTERNAL-2),V(NEXTERNAL-2),IWARN(5)
-      SAVE ACC,ITMAX,IBEGIN,IWARN
-      DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
-*
-* INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
-      IF(IBEGIN.NE.0) GOTO 103
-      IBEGIN=1
-      TWOPI=8.*DATAN(1.D0)
-      PO2LOG=LOG(TWOPI/4.)
-      Z(2)=PO2LOG
-      DO 101 K=3,10
-  101 Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2))
-      DO 102 K=3,10
-  102 Z(K)=(Z(K)-LOG(DFLOAT(K-1)))
-*
-* CHECK ON THE NUMBER OF PARTICLES
-  103 IF(N.GT.1.AND.N.LT.101) GOTO 104
-      PRINT 1001,N
-      STOP
-*
-* CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
-  104 XMT=0.
-      NM=0
-      DO 105 I=1,N
-      IF(XM(I).NE.0.D0) NM=NM+1
-  105 XMT=XMT+ABS(XM(I))
-      IF(XMT.LE.ET) GOTO 201
-      PRINT 1002,XMT,ET
-      STOP
-*
-* THE PARAMETER VALUES ARE NOW ACCEPTED
-*
-* GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
-  201 DO 202 I=1,N
-         r1=rn(1)
-      C=2.*r1-1.
-      S=SQRT(1.-C*C)
-      F=TWOPI*RN(2)
-      r1=rn(3)
-      r2=rn(4)
-      Q(4,I)=-LOG(r1*r2)
-      Q(3,I)=Q(4,I)*C
-      Q(2,I)=Q(4,I)*S*COS(F)
-  202 Q(1,I)=Q(4,I)*S*SIN(F)
-
-      Zvar=rn(1)
-*
-* CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
-      DO 203 I=1,4
-  203 R(I)=0.
-      DO 204 I=1,N
-      DO 204 K=1,4
-  204 R(K)=R(K)+Q(K,I)
-      RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
-      DO 205 K=1,3
-  205 B(K)=-R(K)/RMAS
-      G=R(4)/RMAS
-      A=1./(1.+G)
-      X=ET/RMAS
-*
-* TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
-      DO 207 I=1,N
-      BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
-      DO 206 K=1,3
-  206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
-  207 P(4,I)=X*(G*Q(4,I)+BQ)
-*
-* CALCULATE WEIGHT AND POSSIBLE WARNINGS
-      WT=PO2LOG
-      IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
-      IF(WT.GE.-180.D0) GOTO 208
-      IF(IWARN(1).LE.5) PRINT 1004,WT
-      IWARN(1)=IWARN(1)+1
-  208 IF(WT.LE. 174.D0) GOTO 209
-      IF(IWARN(2).LE.5) PRINT 1005,WT
-      IWARN(2)=IWARN(2)+1
-*
-* RETURN FOR WEIGHTED MASSLESS MOMENTA
-  209 IF(NM.NE.0) GOTO 210
-* RETURN LOG OF WEIGHT
-      WT=WT
-      RETURN
-*
-* MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
-  210 XMAX=SQRT(1.-(XMT/ET)**2)
-      DO 301 I=1,N
-      XM2(I)=XM(I)**2
-  301 P2(I)=P(4,I)**2
-      ITER=0
-      X=XMAX
-      ACCU=ET*ACC
-  302 F0=-ET
-      G0=0.
-      X2=X*X
-      DO 303 I=1,N
-      E(I)=SQRT(XM2(I)+X2*P2(I))
-      F0=F0+E(I)
-  303 G0=G0+P2(I)/E(I)
-      IF(ABS(F0).LE.ACCU) GOTO 305
-      ITER=ITER+1
-      IF(ITER.LE.ITMAX) GOTO 304
-      PRINT 1006,ITMAX
-      GOTO 305
-  304 X=X-F0/(X*G0)
-      GOTO 302
-  305 DO 307 I=1,N
-      V(I)=X*P(4,I)
-      DO 306 K=1,3
-  306 P(K,I)=X*P(K,I)
-  307 P(4,I)=E(I)
-*
-* CALCULATE THE MASS-EFFECT WEIGHT FACTOR
-      WT2=1.
-      WT3=0.
-      DO 308 I=1,N
-      WT2=WT2*V(I)/E(I)
-  308 WT3=WT3+V(I)**2/E(I)
-      WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
-*
-* RETURN FOR  WEIGHTED MASSIVE MOMENTA
-      WT=WT+WTM
-      IF(WT.GE.-180.D0) GOTO 309
-      IF(IWARN(3).LE.5) PRINT 1004,WT
-      IWARN(3)=IWARN(3)+1
-  309 IF(WT.LE. 174.D0) GOTO 310
-      IF(IWARN(4).LE.5) PRINT 1005,WT
-      IWARN(4)=IWARN(4)+1
-* RETURN LOG OF WEIGHT
-  310 WT=WT
-      RETURN
-*
- 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
- 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',
-     . ' SMALLER THAN TOTAL ENERGY =',D15.6)
- 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
- 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
- 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE',
-     . ' DESIRED ACCURACY =',D15.6)
-      END
-
-      FUNCTION RN(IDUMMY)
-      REAL*8 RN,RAN
-      SAVE INIT
-      DATA INIT /1/
-      IF (INIT.EQ.1) THEN
-        INIT=0
-        CALL RMARIN(1802,9373)
-      END IF
-*
-  10  CALL RANMAR(RAN)
-      IF (RAN.LT.1D-16) GOTO 10
-      RN=RAN
-*
-      END
-
-
-
-      SUBROUTINE RANMAR(RVEC)
-*     -----------------
-* Universal random number generator proposed by Marsaglia and Zaman
-* in report FSU-SCRI-87-50
-* In this version RVEC is a double precision variable.
-      IMPLICIT REAL*8(A-H,O-Z)
-      COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
-      COMMON/ RASET2 / IRANMR,JRANMR
-      SAVE /RASET1/,/RASET2/
-      UNI = RANU(IRANMR) - RANU(JRANMR)
-      IF(UNI .LT. 0D0) UNI = UNI + 1D0
-      RANU(IRANMR) = UNI
-      IRANMR = IRANMR - 1
-      JRANMR = JRANMR - 1
-      IF(IRANMR .EQ. 0) IRANMR = 97
-      IF(JRANMR .EQ. 0) JRANMR = 97
-      RANC = RANC - RANCD
-      IF(RANC .LT. 0D0) RANC = RANC + RANCM
-      UNI = UNI - RANC
-      IF(UNI .LT. 0D0) UNI = UNI + 1D0
-      RVEC = UNI
-      END
- 
-      SUBROUTINE RMARIN(IJ,KL)
-*     -----------------
-* Initializing routine for RANMAR, must be called before generating
-* any pseudorandom numbers with RANMAR. The input values should be in
-* the ranges 0<=ij<=31328 ; 0<=kl<=30081
-      IMPLICIT REAL*8(A-H,O-Z)
-      COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM
-      COMMON/ RASET2 / IRANMR,JRANMR
-      SAVE /RASET1/,/RASET2/
-* This shows correspondence between the simplified input seeds IJ, KL
-* and the original Marsaglia-Zaman seeds I,J,K,L.
-* To get the standard values in the Marsaglia-Zaman paper (i=12,j=34
-* k=56,l=78) put ij=1802, kl=9373
-      I = MOD( IJ/177 , 177 ) + 2
-      J = MOD( IJ     , 177 ) + 2
-      K = MOD( KL/169 , 178 ) + 1
-      L = MOD( KL     , 169 )
-      DO 300 II = 1 , 97
-        S =  0D0
-        T = .5D0
-        DO 200 JJ = 1 , 24
-          M = MOD( MOD(I*J,179)*K , 179 )
-          I = J
-          J = K
-          K = M
-          L = MOD( 53*L+1 , 169 )
-          IF(MOD(L*M,64) .GE. 32) S = S + T
-          T = .5D0*T
-  200   CONTINUE
-        RANU(II) = S
-  300 CONTINUE
-      RANC  =   362436D0 / 16777216D0
-      RANCD =  7654321D0 / 16777216D0
-      RANCM = 16777213D0 / 16777216D0
-      IRANMR = 97
-      JRANMR = 33
-      END
-
-
-
-
-
-

=== modified file 'Template/SubProcesses/cuts.f'
--- Template/SubProcesses/cuts.f	2011-09-16 05:33:32 +0000
+++ Template/SubProcesses/cuts.f	2012-02-04 00:00:25 +0000
@@ -53,7 +53,7 @@
       LOGICAL DEBUG
       integer i,j,njets,nheavyjets,hardj1,hardj2
       REAL*8 XVAR,ptmax1,ptmax2,htj,tmp,inclht
-      real*8 ptemp(0:3)
+      real*8 ptemp(0:3), ptemp2(0:3)
       character*20 formstr
 C
 C     PARAMETERS
@@ -207,7 +207,7 @@
 
          if(fixed_ren_scale) then
             G = SQRT(4d0*PI*ALPHAS(scale))
-            call setpara('param_card.dat',.false.)
+            call update_as_param()
          endif
 
 c     Put momenta in the common block to zero to start
@@ -262,12 +262,13 @@
          endif
       enddo
 c
-c    missing ET min & max cut 
-c    nb: missing et simply defined as the sum over the neutrino's 4 momenta
+c    missing ET min & max cut + Invariant mass of leptons and neutrino 
+c    nb: missing Et defined as the vector sum over the neutrino's pt
 c
 c-- reset ptemp(0:3)
       do j=0,3
-         ptemp(j)=0
+         ptemp(j)=0 ! for the neutrino
+         ptemp2(j)=0 ! for the leptons
       enddo
 c-  sum over the momenta
       do i=nincoming+1,nexternal
@@ -276,11 +277,19 @@
             do j=0,3
                ptemp(j)=ptemp(j)+p(j,i)
             enddo
+         elseif(is_a_l(i)) then            
+         if(debug) write (*,*) i,' -> lepton '
+            do j=0,3
+               ptemp2(j)=ptemp2(j)+p(j,i)
+            enddo
          endif
+
       enddo
 c-  check the et
       if(debug.and.ptemp(0).eq.0d0) write (*,*) 'No et miss in event'
       if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Et miss =',pt(ptemp(0)),'   ',misset,':',missetmax
+      if(debug.and.ptemp2(0).eq.0d0) write (*,*) 'No leptons in event'
+      if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Energy of leptons =',pt(ptemp2(0))
       if(ptemp(0).gt.0d0) then
          notgood=(pt(ptemp(0)) .lt. misset).or.
      &        (pt(ptemp(0)) .gt. missetmax)
@@ -290,6 +299,13 @@
             return
          endif
       endif
+      if (mmnl.gt.0d0.or.mmnlmax.lt.1d5)then
+         if(dsqrt(SumDot(ptemp,ptemp2,1d0)).lt.mmnl.or.dsqrt(SumDot(ptemp, ptemp2,1d0)).gt.mmnlmax) then
+            if(debug) write (*,*) 'lepton invariant mass -> fails'
+            passcuts=.false.
+            return
+         endif
+      endif
 c
 c     pt cut on heavy particles
 c     gives min(pt) for (at least) one heavy particle
@@ -680,7 +696,7 @@
                enddo
             enddo
          endif
-         call setpara('param_card.dat',.false.)
+         call update_as_param()
       endif
 
       IF (FIRSTTIME2) THEN

=== removed file 'Template/SubProcesses/driver.f'
--- Template/SubProcesses/driver.f	2011-04-26 18:45:41 +0000
+++ Template/SubProcesses/driver.f	1970-01-01 00:00:00 +0000
@@ -1,289 +0,0 @@
-      Program DRIVER
-c**************************************************************************
-c     This is the driver for the whole calulation
-c**************************************************************************
-      implicit none
-C
-C     CONSTANTS
-C
-      double precision zero
-      parameter       (ZERO = 0d0)
-      include 'genps.inc'
-      include 'maxconfigs.inc'
-      include 'nexternal.inc'
-      INTEGER    ITMAX,   NCALL
-C
-C     LOCAL
-C
-      integer i,ninvar,nconfigs,j,l,l1,l2,ndim
-      double precision dsig,tot,mean,sigma
-      integer npoints,lunsud
-      double precision x,y,jac,s1,s2,xmin
-      external dsig
-      character*130 buf
-      integer NextUnopen
-      external NextUnopen
-c
-c     Global
-c
-      integer                                      nsteps
-      character*40          result_file,where_file
-      common /sample_status/result_file,where_file,nsteps
-      integer           Minvar(maxdim,lmaxconfigs)
-      common /to_invar/ Minvar
-      integer ngroup
-      common/to_group/ngroup
-      data ngroup/0/
-cc
-      include 'run.inc'
-      
-      integer           mincfig, maxcfig
-      common/to_configs/mincfig, maxcfig
-
-
-      double precision twgt, maxwgt,swgt(maxevents)
-      integer                             lun, nw
-      common/to_unwgt/twgt, maxwgt, swgt, lun, nw
-
-c--masses
-      double precision pmass(nexternal)
-      common/to_mass/  pmass
-      double precision qmass(2)
-      common/to_qmass/  qmass
-
-c     $B$ new_def $E$  this is a tag for MadWeigth, Don't edit this line
-
-c      double precision xsec,xerr
-c      integer ncols,ncolflow(maxamps),ncolalt(maxamps),ic
-c      common/to_colstats/ncols,ncolflow,ncolalt,ic
-
-      include 'coupl.inc'
-
-C-----
-C  BEGIN CODE
-C-----  
-c
-c     Read process number
-c
-      open (unit=lun+1,file='../dname.mg',status='unknown',err=11)
-      read (lun+1,'(a130)',err=11,end=11) buf
-      l1=index(buf,'P')
-      l2=index(buf,'_')
-      if(l1.ne.0.and.l2.ne.0.and.l1.lt.l2-1)
-     $     read(buf(l1+1:l2-1),*,err=11) ngroup
- 11   print *,'Process in group number ',ngroup
-
-      lun = 27
-      twgt = -2d0            !determine wgt after first iteration
-      open(unit=lun,status='scratch')
-      nsteps=2
-      call setrun                !Sets up run parameters
-c     $B$ setpara $B$ ! this is a tag for MadWeight. Don't edit this line
-      call setpara('param_card.dat',.true.)   !Sets up couplings and masses
-c     $E$ setpara $E$ ! this is a tag for MadWeight. Don't edit this line
-      include 'pmass.inc'        !Sets up particle masses
-      include 'qmass.inc'        !Sets up particle masses inside onium state
-      call setcuts               !Sets up cuts 
-      call printout              !Prints out a summary of paramaters
-      call run_printout          !Prints out a summary of the run settings
-      nconfigs = 1
-
-c   If CKKW-type matching, read IS Sudakov grid
-      if(ickkw.eq.2 .and. (lpp(1).ne.0.or.lpp(2).ne.0))then
-        lunsud=NextUnopen()
-        open(unit=lunsud,file=issgridfile,status='old',ERR=20)
-        goto 40
- 20     issgridfile='lib/'//issgridfile
-        do i=1,5
-          open(unit=lunsud,file=issgridfile,status='old',ERR=30)          
-          exit
- 30       issgridfile='../'//issgridfile
-          if(i.eq.5)then
-            print *,'ERROR: No Sudakov grid file found in lib with ickkw=2'
-            stop
-          endif
-        enddo
-        print *,'Reading Sudakov grid file ',issgridfile
- 40     call readgrid(lunsud)
-        print *,'Done reading IS Sudakovs'
-      endif
-        
-      if(ickkw.eq.2)then
-        hmult=.false.
-        if(ngroup.ge.nhmult) hmult=.true.
-        if(hmult)then
-          print *,'Running CKKW as highest mult sample'
-        else
-          print *,'Running CKKW as lower mult sample'
-        endif
-      endif
-
-c     
-c     Get user input
-c
-      write(*,*) "getting user params"
-      call get_user_params(ncall,itmax,mincfig)
-      maxcfig=mincfig
-      minvar(1,1) = 0              !This tells it to map things invarients
-      write(*,*) 'Attempting mappinvarients',nconfigs,nexternal
-      call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
-      write(*,*) "Completed mapping",nexternal
-      ndim = 3*(nexternal-2)-4
-      if (abs(lpp(1)) .ge. 1) ndim=ndim+1
-      if (abs(lpp(2)) .ge. 1) ndim=ndim+1
-      ninvar = ndim
-      do j=mincfig,maxcfig
-         if (abs(lpp(1)) .ge. 1 .and. abs(lpp(1)) .ge. 1) then
-            minvar(ndim-1,j)=ninvar-1
-            minvar(ndim,j) = ninvar
-         elseif (abs(lpp(1)) .ge. 1 .or. abs(lpp(1)) .ge. 1) then
-            minvar(ndim,j) = ninvar
-         endif
-      enddo
-      write(*,*) "about to integrate ", ndim,ncall,itmax,ninvar,nconfigs
-      call sample_full(ndim,ncall,itmax,dsig,ninvar,nconfigs)
-c
-c     Now write out events to permanent file
-c
-      if (twgt .gt. 0d0) maxwgt=maxwgt/twgt
-      write(lun,'(a,f20.5)') 'Summary', maxwgt
-      
-c      call store_events
-
-c      write(*,'(a34,20I7)'),'Color flows originally chosen:   ',
-c     &     (ncolflow(i),i=1,ncols)
-c      write(*,'(a34,20I7)'),'Color flows according to diagram:',
-c     &     (ncolalt(i),i=1,ncols)
-c
-c      call sample_result(xsec,xerr)
-c      write(*,*) 'Final xsec: ',xsec
-
-      rewind(lun)
-
-      close(lun)
-      end
-
-c     $B$ get_user_params $B$ ! tag for MadWeight
-c     change this routine to read the input in a file
-c
-      subroutine get_user_params(ncall,itmax,iconfig)
-c**********************************************************************
-c     Routine to get user specified parameters for run
-c**********************************************************************
-      implicit none
-c
-c     Constants
-c
-      include 'nexternal.inc'
-      include 'maxparticles.inc'
-c
-c     Arguments
-c
-      integer ncall,itmax,iconfig
-c
-c     Local
-c
-      integer i, j, jconfig, ncode
-      double precision dconfig
-c
-c     Global
-c
-      integer           isum_hel
-      logical                   multi_channel
-      common/to_matrix/isum_hel, multi_channel
-      double precision    accur
-      common /to_accuracy/accur
-      integer           use_cut
-      common /to_weight/use_cut
-
-      integer        lbw(0:nexternal)  !Use of B.W.
-      common /to_BW/ lbw
-
-c-----
-c  Begin Code
-c-----
-      write(*,'(a)') 'Enter number of events and iterations: '
-      read(*,*) ncall,itmax
-      write(*,*) 'Number of events and iterations ',ncall,itmax
-      write(*,'(a)') 'Enter desired fractional accuracy: '
-      read(*,*) accur
-      write(*,*) 'Desired fractional accuracy: ',accur
-
-      write(*,'(a)') 'Enter 0 for fixed, 2 for adjustable grid: '
-      read(*,*) use_cut
-      if (use_cut .lt. 0 .or. use_cut .gt. 2) then
-         write(*,*) 'Bad choice, using 2',use_cut
-         use_cut = 2
-      endif
-
-      write(*,10) 'Suppress amplitude (0 no, 1 yes)? '
-      read(*,*) i
-      if (i .eq. 1) then
-         multi_channel = .true.
-         write(*,*) 'Using suppressed amplitude.'
-      else
-         multi_channel = .false.
-         write(*,*) 'Using full amplitude.'
-      endif
-
-      write(*,10) 'Exact helicity sum (0 yes, n = number/event)? '
-      read(*,*) i
-      if (i .eq. 0) then
-         isum_hel = 0
-         write(*,*) 'Explicitly summing over helicities'
-      else
-         isum_hel= i
-         write(*,*) 'Summing over',i,' helicities/event'
-      endif
-
-      write(*,10) 'Enter Configuration Number: '
-      read(*,*) dconfig
-c     ncode is number of digits needed for the BW code
-      ncode=int(dlog10(3d0)*(max_particles-3))+1
-      iconfig = int(dconfig*(1+10**(-ncode)))
-      write(*,12) 'Running Configuration Number: ',iconfig
-c
-c     Here I want to set up with B.W. we map and which we don't
-c
-      dconfig = dconfig-iconfig
-      if (dconfig .eq. 0) then
-         write(*,*) 'Not subdividing B.W.'
-         lbw(0)=0
-      else
-         lbw(0)=1
-         jconfig=dconfig*(10**ncode + 0.1)
-         write(*,*) 'Using dconfig=',jconfig
-         call DeCode(jconfig,lbw(1),3,nexternal)
-         write(*,*) 'BW Setting ', (lbw(j),j=1,nexternal-2)
-c         do i=nexternal-3,0,-1
-c            if (jconfig .ge. 2**i) then
-c               lbw(i+1)=1
-c               jconfig=jconfig-2**i
-c            else
-c               lbw(i+1)=0
-c            endif 
-c            write(*,*) i+1, lbw(i+1)
-c         enddo
-      endif
- 10   format( a)
- 12   format( a,i4)
-      end
-c     $E$ get_user_params $E$ ! tag for MadWeight
-c     change this routine to read the input in a file
-c
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-

=== modified file 'Template/SubProcesses/reweight.f'
--- Template/SubProcesses/reweight.f	2011-10-05 03:10:17 +0000
+++ Template/SubProcesses/reweight.f	2012-02-04 00:00:25 +0000
@@ -723,6 +723,9 @@
       integer              ISPROC 
       DOUBLE PRECISION PD(0:MAXPROC)
       COMMON /SubProc/ PD, ISPROC
+      INTEGER SUBDIAG(MAXSPROC),IB(2)
+      COMMON/TO_SUB_DIAG/SUBDIAG,IB
+      data IB/1,2/
 c     q2bck holds the central q2fact scales
       real*8 q2bck(2)
       common /to_q2bck/q2bck
@@ -758,12 +761,6 @@
       external alphas, isjetvx, getissud, pdg2pdf, xran1,  sudwgt
 
       rewgt=1.0d0
-
-      if((ickkw.gt.0.or..not.fixed_fac_scale.or..not.fixed_ren_scale)
-     $     .and..not.clustered)then
-        write(*,*)'Error: No clustering done when calling rewgt!'
-        stop
-      endif
       clustered=.false.
 
       if(ickkw.le.0) return
@@ -831,7 +828,7 @@
       endif
 c     Set x values for the two sides, for IS Sudakovs
       do i=1,2
-        xnow(i)=xbk(i)
+        xnow(i)=xbk(ib(i))
       enddo
       if(btest(mlevel,3))then
         write(*,*) 'Set x values to ',xnow(1),xnow(2)
@@ -1005,11 +1002,12 @@
      $                          ' to: ',sqrt(pt2pdf(imocl(n)))
                      else if(pt2pdf(idacl(n,i)).lt.q2now
      $                       .and.isjet(ipdgcl(idacl(n,i),igraphs(1),iproc))) then
-                        pdfj1=pdg2pdf(abs(ibeam(j)),ipdgcl(idacl(n,i),
-     $                       igraphs(1),iproc)*sign(1,ibeam(j)),xnow(j),sqrt(q2now))
-                        pdfj2=pdg2pdf(abs(ibeam(j)),ipdgcl(idacl(n,i),
-     $                       igraphs(1),iproc)*sign(1,ibeam(j)),xnow(j),
-     $                       sqrt(pt2pdf(idacl(n,i))))
+                        pdfj1=pdg2pdf(abs(lpp(IB(j))),ipdgcl(idacl(n,i),
+     $                       igraphs(1),iproc)*sign(1,lpp(IB(j))),
+     $                       xnow(j),sqrt(q2now))
+                        pdfj2=pdg2pdf(abs(lpp(IB(j))),ipdgcl(idacl(n,i),
+     $                       igraphs(1),iproc)*sign(1,lpp(IB(j))),
+     $                       xnow(j),sqrt(pt2pdf(idacl(n,i))))
                         if(pdfj2.lt.1d-10)then
 c                          Scale too low for heavy quark
                            rewgt=0d0
@@ -1026,12 +1024,12 @@
                            write(*,*)'           PDF: ',pdfj1,' / ',pdfj2
                            write(*,*)'        -> rewgt: ',rewgt
 c                           write(*,*)'  (compare for glue: ',
-c     $                          pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i)))),' / ',
-c     $                          pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2ijcl(n)))
+c     $                          pdg2pdf(lpp(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i)))),' / ',
+c     $                          pdg2pdf(lpp(j),21,xbk(j),sqrt(pt2ijcl(n)))
 c                           write(*,*)'       = ',pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i))))/
-c     $                          pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2ijcl(n)))
+c     $                          pdg2pdf(lpp(j),21,xbk(j),sqrt(pt2ijcl(n)))
 c                           write(*,*)'       -> ',pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2pdf(idacl(n,i))))/
-c     $                          pdg2pdf(ibeam(j),21,xbk(j),sqrt(pt2ijcl(n)))*rewgt,' )'
+c     $                          pdg2pdf(lpp(j),21,xbk(j),sqrt(pt2ijcl(n)))*rewgt,' )'
                         endif
 c                       Set scale for mother as this scale
                         pt2pdf(imocl(n))=q2now                           

=== modified file 'Template/SubProcesses/unwgt.f'
--- Template/SubProcesses/unwgt.f	2011-10-03 08:39:00 +0000
+++ Template/SubProcesses/unwgt.f	2012-02-04 00:00:25 +0000
@@ -27,8 +27,8 @@
 C     GLOBAL
 C
       double precision twgt, maxwgt,swgt(maxevents)
-      integer                             lun, nw
-      common/to_unwgt/twgt, maxwgt, swgt, lun, nw
+      integer                             lun, nw, itmin
+      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
       integer nzoom
       double precision  tx(1:3,maxinvar)
       common/to_xpoints/tx, nzoom
@@ -126,8 +126,8 @@
 C     GLOBAL
 C
       double precision twgt, maxwgt,swgt(maxevents)
-      integer                             lun, nw
-      common/to_unwgt/twgt, maxwgt, swgt, lun, nw
+      integer                             lun, nw, itmin
+      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 
 c-----
 c  Begin Code
@@ -163,8 +163,8 @@
 C     GLOBAL
 C
       double precision twgt, maxwgt,swgt(maxevents)
-      integer                             lun, nw
-      common/to_unwgt/twgt, maxwgt, swgt, lun, nw
+      integer                             lun, nw, itmin
+      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 
       double precision    matrix
       common/to_maxmatrix/matrix
@@ -245,8 +245,8 @@
 C     GLOBAL
 C
       double precision twgt, maxwgt,swgt(maxevents)
-      integer                             lun, nw
-      common/to_unwgt/twgt, maxwgt, swgt, lun, nw
+      integer                             lun, nw, itmin
+      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 
       integer                   neventswritten
       common /to_eventswritten/ neventswritten
@@ -269,7 +269,7 @@
 c     First scale all of the events to the total cross section
 c
       if (nw .le. 0) return
-      call sample_result(xsec,xerr)
+      call sample_result(xsec,xerr,itmin)
       if (xsec .le. 0) return   !Fix by TS 12/3/2010
       xtot=0
       call dsort(nw, swgt)
@@ -403,8 +403,8 @@
 C     GLOBAL
 C
       double precision twgt, maxwgt,swgt(maxevents)
-      integer                             lun, nw
-      common/to_unwgt/twgt, maxwgt, swgt, lun, nw
+      integer                             lun, nw, itmin
+      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 
       integer              IPROC 
       DOUBLE PRECISION PD(0:MAXPROC)
@@ -596,8 +596,8 @@
 C     GLOBAL
 C
       double precision twgt, maxwgt,swgt(maxevents)
-      integer                             lun, nw
-      common/to_unwgt/twgt, maxwgt, swgt, lun, nw
+      integer                             lun, nw, itmin
+      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 c-----
 c  Begin Code
 c-----

=== removed file 'Template/bin/TheChopper-pl'
--- Template/bin/TheChopper-pl	2010-10-30 03:26:37 +0000
+++ Template/bin/TheChopper-pl	1970-01-01 00:00:00 +0000
@@ -1,118 +0,0 @@
-#!/usr/bin/perl -w
-#------------------------------------------------------------
-# Simple script to chop out sub-processes  which have a small
-# impact on the cross section, in the grid pack.
-# This is meant to be run after the gridpack warming up run
-# and before packing it up.
-#
-# The script should be run from the main directory of madevent
-# and sits in the bin directory.
-#
-# This is still version 0, where no classes can be defined.
-# However, the result is a script, chop_procs sitting in 
-# SubProcesses, that can be edited by hand. 
-#
-# Author: Fabio Maltoni
-# Date: 10/4/08
-# Version: 0.1
-#
-# HISTORY
-# 
-# 0.2 Fix bug in writing subproc.mg
-#
-# 0.1 Added the possibility to list SubProc IDs
-#
-# 0.0 Original
-#
-#------------------------------------------------------------
-
-#------------------------------------------------------------
-# Auxiliary Ordering function
-#------------------------------------------------------------
-
-    sub OrderInc {
-     $all{$a} <=> $all{$b};
-    };
-
-#------------------------------------------------------------
-# Input the chopping value. Example : 0.01 = 1%
-#------------------------------------------------------------
-
-    print "input chopping off value (e.g., 0.01=1% ):\n";
-    $cutoff_lim=<STDIN>;
-
-#------------------------------------------------------------
-# Input the SubProcesses ids on which the script should run 
-#------------------------------------------------------------
-    
-    print "Input SubProcesses ids to process (space separated):\n";
-    @id_list=split('\s+',<STDIN>);
-    
-
-#------------------------------------------------------------
-# Parse the form, look for known variables
-#------------------------------------------------------------
-
-   
-    
-    chdir("SubProcesses");
-    open(SCRIPT,">chop_procs")  or die "Could not write chop_procs \n";
-     
-    foreach $id (@id_list) {
-    
-    print "\n\nNow processing SubProcesses with ID = $id\n";
-    
-    my @all_Pdirs = <P$id*>;
-    my $xsec = 0;
-    my $perc = 0;
-    %all = ();
-    %below = ();
-    
-    foreach $arg (@all_Pdirs) {
-    chdir("$arg");
-    open(RES,"<default_results.dat")  or die "Could not read default_results.dat\n";
-    my $buff;
-    my $number_read = read(RES,$buff,12);
-    print "$arg has xsec $buff\n";
-    $all{$arg} .= $buff;
-    $xsec = $xsec + $buff;
-    close(RES);
-    chdir("..");
-    };     	
-
-    print " total xsec for ID = $id : $xsec\n";
-    $cutoff=$cutoff_lim*$xsec;
-    print " cutoff for ID = $id : $cutoff \n";
-
-
-#   now I have to order them and find those procs with fall
-#   below xsec*cutoff.    
-
-    $i=0;
-    $integ[0]=0;
-    print "Processes in order FOR ID = $id:\n";
-    foreach $key (sort OrderInc (keys(%all))) {
-    $integ[$i+1] = $integ[$i] + $all{$key}; 	
-    if($integ[$i+1] < $cutoff) { 
-	$below{$key} .= $all{$key};
-	print "$integ[$i+1] : $all{$key} $key => below cutoff\n";}
-    else {
-	print "$integ[$i+1] : $all{$key} $key => above cutoff\n";}
-    $i =$i + 1;
-};
-#
-#  Now I move the corresponding procs out of the way
-#
-    print "The following SubProcess will be tagged to be removed:\n";
-    foreach $key (keys(%below)){
-	print "$key ";
-	print SCRIPT "mv $key XXX$key \# $below{$key}\n";
-#        system("mv $key XXX$key "); 
-    };
-#	print SCRIPT "ls -d P$id* >subproc.mg\n";
-#        system("ls -d P* >subproc.mg");  
- }
- print SCRIPT "ls -d P* >subproc.mg";  
- close(SCRIPT);
- system("chmod +x chop_procs");
-   

=== removed file 'Template/bin/addmasses.py'
--- Template/bin/addmasses.py	2010-10-30 03:26:37 +0000
+++ Template/bin/addmasses.py	1970-01-01 00:00:00 +0000
@@ -1,320 +0,0 @@
-#!/usr/bin/python
-#
-# doubleswitch.py
-# Author: Stephen Mrenna
-# Add masses to charged leptons, fix kinematics
-# Insert W/Z when missing
-
-import math,sys,re
-from xml.dom import minidom
-
-#tolerance for energy momentum conservation
-toler = 1e-4
-#lepton masses from PDGLive
-massDict = {}
-massDict[11]=0.000510998910
-massDict[-11]=massDict[11]
-massDict[13]=105.6583668e-3
-massDict[-13]=massDict[13]
-massDict[15]=1.77684
-massDict[-15]=massDict[15]
-
-#print out messages if true
-verbose=0
-
-#default is NOT to insert missing mothers
-motherFlag=1
-
-#useful class to describe 4 momentum
-class Momentum:
-    def __init__(self,px,py,pz,E,m):
-        self.px=px
-        self.py=py
-        self.pz=pz
-        self.E=E
-        self.m=m
-    def __add__(self,other):
-        t=Momentum(self.px+other.px,self.py+other.py,self.pz+other.pz,self.E+other.E,0)
-        t.m=t.calcMass()
-        return t
-    def __sub__(self,other):
-        t=Momentum(self.px-other.px,self.py-other.py,self.pz-other.pz,self.E-other.E,0)
-        t.m=t.calcMass()
-        return t
-    def calcMass(self):
-        tempMass2=self.E**2-self.px**2-self.py**2-self.pz**2
-        if tempMass2 > 0:
-            t=math.sqrt(tempMass2)
-            if t>toler:
-                return t
-            else:
-                return 0
-        else:
-            return 0
-    def boost(self,ref,rdir):
-        pmag=ref.E
-        DBX=ref.px*rdir/pmag
-        DBY=ref.py*rdir/pmag
-        DBZ=ref.pz*rdir/pmag
-        DB=math.sqrt(DBX**2+DBY**2+DBZ**2)
-        DGA=1.0/math.sqrt(1.0-DB**2)        
-        DBP=DBX*self.px+DBY*self.py+DBZ*self.pz
-        DGABP=DGA*(DGA*DBP/(1.0+DGA)+self.E)
-        self.px = self.px+DGABP*DBX
-        self.py = self.py+DGABP*DBY
-        self.pz = self.pz+DGABP*DBZ
-        self.E  = DGA*(self.E+DBP)
-    def reScale(self,pi,po):
-        self.px = self.px/pi*po
-        self.py = self.py/pi*po
-        self.pz = self.pz/pi*po
-    def printMe(self):
-        li = [self.px,self.py, self.pz, self.E, self.m]
-        print "| %18.10E %18.10E %18.10E %18.10E %18.10E |" % tuple(li)    
-
-#useful class to describe a particle
-class Particle:
-    def __init__(self,i,l):
-        self.no = i
-        self.id = l[0]
-        self.status = l[1]
-        self.mo1 = l[2]
-        self.mo2 = l[3]
-        self.co1 = l[4]
-        self.co2 = l[5]
-        self.mom = Momentum(l[6],l[7],l[8],l[9],l[10])
-        self.life = l[11]
-        self.polar = l[12]
-    def printMe(self):
-        li = [self.no, self.id, self.status,self.mo1, self.mo2, self.co1, self.co2, self.mom.px,self.mom.py, self.mom.pz, self.mom.E, self.mom.m, self.life, self.polar]
-        print "%2i | %9i | %4i | %4i %4i | %4i %4i | %18.10E %18.10E %18.10E %18.10E %18.10E | %1.0f. %2.0f" % tuple(li)
-    def writeMe(self):
-        li = [self.id, self.status,self.mo1, self.mo2, self.co1, self.co2, self.mom.px,self.mom.py, self.mom.pz, self.mom.E, self.mom.m, self.life, self.polar]
-        return "%9i %4i %4i %4i %4i %4i %18.10E %18.10E %18.10E %18.10E %18.10E  %1.0f. %2.0f\n" % tuple(li)        
-
-#useful function for converting a string to variables
-def parseStringToVars(input):
-    if input.find(".")>-1 :
-        return float(input)
-    else:
-        return int(input)
-
-#main part of analysis
-if len(sys.argv)!=3:
-    print "Usage: addmasses.py <infile> <outfile>    "
-    print " Last modified: Fri Nov 21 10:49:14 CST 2008 "
-    sys.exit(0)
-else:
-    print "Running addmasses.py to add masses and correct kinematics of fixed particles"
-    
-#first print out leading information
-try:
-    f=open(sys.argv[1],'r')
-except IOError:
-    print "need a file for reading"
-    sys.exit(0)
-
-try:
-    g=open(sys.argv[2],'w')
-except IOError:
-    print "need a file for writing"
-    sys.exit(0)
-
-while 1:
-    try:
-        line=f.readline()
-    except IOError:
-        print "Problem reading from file ",sys.argv[1]
-        sys.exit(0)
-    if line.find("<event>")==-1:
-        g.write(line)
-    else:
-        break
-
-f.close()
-
-#let xml find the event tags
-try:
-    xmldoc = minidom.parse(sys.argv[1])
-except IOError:
-    print " could not open file for xml parsing ",sys.argv[1]
-    sys.exit(0)
-    
-    
-reflist = xmldoc.getElementsByTagName('event')    
-
-for ref in reflist:
-    lines = ref.toxml()
-    slines = lines.split("\n")
-    next = 0
-    nup = 0
-    nlines = len(slines)
-    counter = 0
-    event = []
-    event_description=""
-    event_poundSign=""
-
-    while counter<nlines:
-        s=slines[counter]
-        if s.find("<event>")>-1:
-            next=1
-        elif s.find("</event>")>-1:
-            pass
-        elif s.find("#")>-1:
-            event_poundSign=s
-        elif next==1:
-            event_description=s
-            next=0
-        elif not s:
-            continue
-        else:
-            t=[]
-            for l in s.split(): t.append(parseStringToVars(l))
-            nup = nup+1
-            event.append(Particle(nup,t))
-        counter=counter+1
-
-
-#default is to skip this
-    if motherFlag:
-
-        noMotherList=[]
-        for p in event:
-            if p.status == -1: continue
-            idabs = abs(p.id)
-            idmo = p.mo1
-            pmo = event[idmo-1]
-            if idabs>=11 and idabs<=16:
-                if p.mo1==1 or (pmo.co1 !=0 or pmo.co2 !=0):
-                    noMotherList.append(p.no)
-            elif idabs<=5 and abs(pmo.id)==6:
-                if not ( p.co1==pmo.co1 and p.co2==pmo.co2):
-                    noMotherList.append(p.no)                
-
-        nAdded=0
-        if len(noMotherList)==0:
-            pass
-        elif len(noMotherList)%2 != 0:
-            print "single orphan; do not know how to process"
-        else:
-            ki=0
-            while ki<len(noMotherList)-1:
-                l1=event[noMotherList[ki]-1]
-                ki=ki+1
-                l2=event[noMotherList[ki]-1]
-                lq=l1.id + l2.id
-                if lq==-1 or lq==1:
-                    idMom=24*lq
-                elif lq==0:
-                    idMom=23
-                else:
-                    break
-                lMom=l1.mom+l2.mom            
-                lm=[idMom,2,l1.mo1,l1.mo2,0,0,lMom.px,lMom.py,lMom.pz,lMom.E,lMom.calcMass(),0,0]
-                nup=nup+1
-                nAdded=nAdded+1
-                event.append(Particle(nup,lm))
-                l1.mo1=l1.mo2=l2.mo1=l2.mo2=nup
-                ki=ki+1
-
-# update number of partons if mothers added
-            if nAdded>0:
-                s1=event_description.split()[0]
-                mySub = re.compile(s1)
-                event_description = mySub.sub(str(nup),event_description)
-
-
-    if nAdded>0:
-        totalLength = len(event)+1
-        newPosition = []
-        iPos=0
-        while iPos<totalLength:           
-            for p in event:
-                if p.mo1==iPos:
-                    newPosition.append(p.no)
-            iPos=iPos+1
-
-        iUp=1
-        event0=[]
-        for ip in newPosition:
-            event0.append(event[ip-1])
-
-        iUp=0
-        for p in event0:
-            iUp=iUp+1
-            if p.no != iUp:
-                ip=iUp                
-                while ip<totalLength-1:
-                    pp=event0[ip]
-                    if pp.mo1==p.no and pp.mo1>pp.no:
-                        pp.mo1=pp.mo2=iUp
-                    ip=ip+1
-            p.no=iUp
-
-        iUp=0
-        for p in event0:
-            event[iUp]=event0[iUp]
-            iUp=iUp+1
-
-
-
-
-# identify mothers
-    particleDict={}
-    for p in event:
-        idabs = abs(p.id)
-        if idabs>=11 and idabs<=16:
-            if p.mo1==1:
-                pass
-            else:
-                if p.mo1 in particleDict:
-                    l=[particleDict[p.mo1]]
-                    l.append(p.no)
-                else:
-                    l=p.no
-                particleDict[p.mo1]=l
-
-# repair kinematics
-    for k in particleDict:
-        t=particleDict[k]
-        p1=event[t[0]-1]
-        p1.mom.boost(event[k-1].mom,-1)
-        p2=event[t[1]-1]
-        p2.mom.boost(event[k-1].mom,-1)
-        rsh=event[k-1].mom.m
-        if p1.id in massDict: p1.mom.m=massDict[p1.id]
-        if p2.id in massDict: p2.mom.m=massDict[p2.id]
-        p1.mom.E = (rsh*rsh + p1.mom.m**2 - p2.mom.m**2)/(2.0*rsh)
-        pmagOld=math.sqrt(p1.mom.px**2+p1.mom.py**2+p1.mom.pz**2)
-        pmagNew=math.sqrt(p1.mom.E**2-p1.mom.m**2)
-        p1.mom.reScale(pmagOld,pmagNew)
-        p2.mom = Momentum(0,0,0,rsh,rsh) - p1.mom
-        p1.mom.boost(event[k-1].mom,1)
-        p2.mom.boost(event[k-1].mom,1)
-
-    pSum = Momentum(0,0,0,0,0)
-    for p in event:
-        if p.status== 2 :
-            pass
-        elif p.status==-1:
-            pSum = pSum - p.mom
-        elif p.status==1:
-            pSum = pSum + p.mom
-
-    if abs(pSum.px)>toler or abs(pSum.py)>toler or abs(pSum.pz)>toler or abs(pSum.E)>toler:
-        print "Event does not pass tolerance ",toler
-        pSum.printMe()
-
-    if 1:
-        g.write("<event>\n")
-        g.write(event_description+"\n")
-        for p in event:
-            g.write(p.writeMe())
-        g.write(event_poundSign+"\n")
-        g.write("</event>\n")
-#at the end
-g.write("</LesHouchesEvents>\n")
-g.close()
-
-    
-

=== removed file 'Template/bin/calculate_crossx'
--- Template/bin/calculate_crossx	2010-10-30 03:26:37 +0000
+++ Template/bin/calculate_crossx	1970-01-01 00:00:00 +0000
@@ -1,140 +0,0 @@
-#!/bin/bash
-#
-#  This runs survey,refine,unweight_events, one after the other
-#  
-#   First we need to get into the main directory
-#
-#
-#  Usage: generate_events compression events parallel [name]
-#
-
-if [[  ! -d ./bin  ]]; then
-    cd ../
-    if [[ ! -d ./bin ]]; then
-        echo "Error: it must be executed from the main, or bin directory"
-        exit
-    fi
-fi
-#
-#  Now let shell know where to find important executables
-#
-main=`pwd`
-dirbin=$main/bin
-pydir=$main/../pythia-pgs/src
-pgsdir=$pydir
-ERAdir=$main/../ExRootAnalysis
-webbin=$dirbin
-td=$dirbin/td
-web=0
-
-echo $$ >> myprocid
-
-#if ( "$1" == "" ) then
-#    echo 'Number of unweighted events. This is ingnored and read from run_card.dat ' 
-#    set a  = $< 
-#else
-#    set a = $1
-#endif
-if [[  "$1" == ""  ]]; then
-    echo 'Enter 1 for parallel 0 for serial run' 
-     read p 
-else
-     p=$1
-fi
-n=MadEvent
-if [[  $p -gt 0 ]]; then
-   if [[  "$2" == ""  ]]; then
-    echo 'Enter name for jobs on pbs queue'
-       read n
-   else
-       n=$2
-   fi
-   if [[  "$3" == ""  ]]; then
-       echo 'Enter run name'
-       read t
-   else
-       t=$3
-   fi
-else
-   if [[  "$2" == ""  ]]; then
-      echo 'Enter run name'
-      read t
-   else
-      t=$2
-   fi
-fi
-#set t = TeV2
-if [[  ${#argv} -gt 3  ]]; then
-    web=1
-    webbin="$MADGRAPH_BASE/MG_ME/WebBin"
-    pydir="$webbin/pythia-pgs"
-    pgsdir=$pydir
-    ERAdir="$MADGRAPH_BASE/MG_ME/ExRootAnalysis"
-    td=$webbin/td
-    touch Online
-fi
-date
-a=`awk '/^[^#].*=.*nevents/{print $1}' Cards/run_card.dat`
-echo Generating $a events
-#
-#  Check if run already exists. If so, store run w/ new name
-#  and remove old run before starting.
-#
-
-if [[  -e status  ]]; then
-  rm status
-fi
-if [[  -e error  ]]; then
-  rm error
-fi
-touch RunWeb
-echo "Cleaning directories" > status
-$dirbin/gen_crossxhtml-pl $t
-$dirbin/clean
-touch survey
-echo "Starting jobs" > status
-$dirbin/gen_crossxhtml-pl $t
-$dirbin/survey $p $n $t
-if [[  -e error  ]]; then
-  cat error
-  date
-  cp -f error status
-  rm refine
-  rm RunWeb
-  $dirbin/gen_crossxhtml-pl $t
-  $dirbin/gen_cardhtml-pl
-  exit
-fi
-#
-#  Now collect the events - just to get the banner in fact
-#
-echo "Combining Events" >& status
-echo "Combining Events" 
-$dirbin/gen_crossxhtml-pl $t
-pushd ./Source > /dev/null
-make ../bin/combine_events
-popd > /dev/null
-pushd SubProcesses > /dev/null
-$dirbin/run_combine $p
-mv events.lhe ../Events/
-mv unweighted_events.lhe ../Events/
-popd > /dev/null
-#
-#  do the banner
-#
-cd ./Events
-echo "putting the banner"
-$dirbin/put_banner events.lhe
-$dirbin/put_banner unweighted_events.lhe
-cd ../
-#
-#  Store Events
-#
-echo "Storing Events" >& status
-$dirbin/gen_crossxhtml-pl $t
-$dirbin/store $t
-rm -f RunWeb
-echo " " >& status
-$dirbin/gen_crossxhtml-pl $t
-$dirbin/gen_cardhtml-pl
-date

=== removed file 'Template/bin/clean4grid'
--- Template/bin/clean4grid	2010-10-30 03:26:37 +0000
+++ Template/bin/clean4grid	1970-01-01 00:00:00 +0000
@@ -1,95 +0,0 @@
-#!/bin/bash
-#
-################################################################################
-#                                                                             ##
-#                    MadGraph/MadEvent                                        ##
-#                                                                             ##
-# FILE : clean4grid                                                           ##
-# VERSION : 1.0                                                               ##
-# DATE : 29 Jan 2008                                                          ##
-# AUTHORS : MH & RF - MadGraph team                                           ##
-#                                                                             ##
-# DESCRIPTION : script to clean package for the grid                          ##
-# USAGE : ./clean4grid                                                        ##
-################################################################################
-
-set nonomatch
-
-if [[ ! -d ./bin ]]; then
-    cd ../
-    if [[ ! -d ./bin ]]; then
-        echo "Error: it must be executed from the main, or bin directory"
-        exit
-    fi
-fi
-
-if [[ -d SubProcesses ]]; then
-    cd SubProcesses
-    echo -n "Cleaning  SubProcesses"
-    rm -f *.f  >& /dev/null
-    for i in `cat subproc.mg` ; do
-	cd $i
-	echo -n "."
-	rm -f *.o >& /dev/null
-        rm -f *.f  >& /dev/null
-#	rm -f madevent >& /dev/null
-#	rm -f *ajob* >& /dev/null
-	rm -f gensym >& /dev/null
-        rm -f matrix*.jpg >& /dev/null
-        rm -f matrix.ps >& /dev/null
-	rm -f G*/ftn25 >& /dev/null
-	rm -f G*/ftn26 >& /dev/null
-	rm -f G*/events.lhe >& /dev/null
-	cd ..
-    done
-    echo " "
-    cd ../
-else
-    echo "Error could not find SubProcesses"
-    exit
-fi
-if [[ -d Source ]]; then
-    echo "Removing Source:"
-    cd Source
-    rm -f *.o >& /dev/null
-    rm -f *.f >& /dev/null
-    rm -f *.inc >& /dev/null
-    rm -f DHELAS/*.F >& /dev/null
-    rm -f DHELAS/*.f >& /dev/null
-    rm -f DHELAS/*.o >& /dev/null
-    rm -f DHELAS/*.inc >& /dev/null
-    rm -rf DHELAS/Doc >& /dev/null
-    rm -f MODEL/*.f >& /dev/null
-    rm -f MODEL/*.o >& /dev/null
-    rm -f MODEL/*.dat >& /dev/null
-    rm -f MODEL/*.html >& /dev/null
-    rm -f PDF/*.f >& /dev/null
-    rm -f PDF/*.o >& /dev/null
-    rm -f PDF/*.inc >& /dev/null
-    cd ..
-fi
-if [[ -d HTML ]]; then
-    echo "Removing HTML files:"
-    rm -f HTML/*  >& /dev/null
-fi
-#if [[ -d lib ]]; then
-#    cd lib
-#    echo "Cleaning lib:"
-#    rm -f *.a >& /dev/null
-#    cd ../
-#else
-#    echo "Error could not find lib"
-#    exit
-#fi
-#if [[ -d bin ]]; then
-#    cd bin
-#    echo "Cleaning bin:"
-#    for i in gen_ximprove scale_events select_events sum_html combine_events; do
-#	rm -f $i >& /dev/null
-#    done
-#    cd ../
-#else
-#    echo "Error could not find bin"
-#    exit
-#fi
-

=== modified file 'Template/bin/cleanall'
--- Template/bin/cleanall	2011-07-06 17:23:07 +0000
+++ Template/bin/cleanall	2012-02-04 00:00:25 +0000
@@ -5,7 +5,7 @@
 #     executable files
 #     events.dat   (Except in Events)
 #
-#  Usage: clean
+#  Usage: cleanall
 #
 #  
 #   First we need to get into the main directory
@@ -17,6 +17,15 @@
 	exit
     fi
 fi
+if [[ -d Source ]]; then
+    cd Source
+    echo "Cleaning Source:"
+    make clean >& /dev/null
+    cd ../
+else
+    echo "Error could not find Source"
+    exit
+fi
 if [[ -d SubProcesses ]]; then
     cd SubProcesses
     echo -n "Cleaning  SubProcesses"
@@ -35,23 +44,3 @@
     echo "Error could not find SubProcesses"
     exit
 fi
-if [[ -d Source ]]; then
-    cd Source
-    echo "Cleaning Source:"
-    make clean >& /dev/null
-    cd ../
-else
-    echo "Error could not find Source"
-    exit
-fi
-if [[ -d bin ]]; then
-    cd bin
-    echo "Cleaning bin:"
-    for i in gen_ximprove scale_events select_events sum_html combine_runs combine_events; do
-	rm -f $i >& /dev/null
-    done
-    cd ../
-else
-    echo "Error could not find bin"
-    exit
-fi

=== removed file 'Template/bin/compile'
--- Template/bin/compile	2011-07-11 18:32:26 +0000
+++ Template/bin/compile	1970-01-01 00:00:00 +0000
@@ -1,89 +0,0 @@
-#!/bin/bash
-#
-################################################################################
-#                                                                             ##
-#                    MadGraph/MadEvent                                        ##
-#                                                                             ##
-# FILE : compile                                                              ##
-# VERSION : 1.0                                                               ##
-# DATE : 23 Septembre 2007                                                    ##
-# AUTHOR : MH - MadGraph team                                                 ##
-#                                                                             ##
-# DESCRIPTION : script to compile madevent                                    ##
-# USAGE : ./make_package                                                      ##
-################################################################################
-
-# set nonomatch
-
-if [[ ! -d ./bin ]]; then
-    cd ../
-    if [[ ! -d ./bin ]]; then
-        echo "Error: it must be executed from the main, or bin directory"
-        exit
-    fi
-fi
-
-
-#
-# If argument is equal to 'd' use dynamic libraries
-#
-echo $PWD
-if [[ "$1" == "dynamic" ]]; then
-    echo "Using dynamic libraries (might not work under MacOsX)"
-    export dynamic=true
-else
-    echo "Using static libraries"
-    unset dynamic
-fi
-
-# Check for LHAPDF
-c=`awk '/^[^#].*=.*pdlabel/{print $1}' Cards/run_card.dat`
-if [[ $c == "'lhapdf'" ]]; then
-    echo Using LHAPDF interface!
-    export lhapdf=true
-else
-    unset lhapdf
-fi
-
-#
-#  Now let shell know where to find important executables
-#
-
-bin/compile_Source
-
-if [[ -e error ]];then
-   exit
-fi
-
-if [[ -d SubProcesses ]]; then
-    cd SubProcesses
-    for i in `cat subproc.mg` ; do
-	cd $i
-	echo $i
-	rm -f ajob* >& /dev/null
-	rm -f wait.ajob* >& /dev/null
-	rm -f run.ajob* >& /dev/null
-	rm -f done.ajob* >& /dev/null
-        make madevent > /dev/null
-	if [[ $? -ne 0 ]];then
-	    # Make didn't exit successfully
-	    echo Error make madevent not successful > error
-	    cat error
-	    exit
-        fi
-	cd ..
-     done
-    cd ..
-else
-    echo "Error could not find SubProcesses"
-    exit
-fi
-
-if [[ -d ../DECAY ]]; then
-    echo "DECAY directory found, compiling..."
-    cd ../DECAY
-    make
-fi
-
-echo ""
-

=== removed file 'Template/bin/compile_Source'
--- Template/bin/compile_Source	2011-07-06 22:40:17 +0000
+++ Template/bin/compile_Source	1970-01-01 00:00:00 +0000
@@ -1,32 +0,0 @@
-#!/bin/bash
-#
-################################################################################
-#                                                                             ##
-#                    MadGraph/MadEvent                                        ##
-#                                                                             ##
-# FILE : compile_Source                                                       ##
-# VERSION : 1.0                                                               ##
-# DATE : 6 July 2011                                                          ##
-# AUTHOR : JA - MadGraph team                                                 ##
-#                                                                             ##
-# DESCRIPTION : script to compile Source                                      ##
-################################################################################
-
-if [[ -d Source ]]; then
-    cd Source
-    for i in ../bin/sum_html ../bin/gen_ximprove all ../bin/combine_events ../bin/combine_runs ../bin/sumall; do
-	make $i > /dev/null
-	if [[ $? -ne 0 ]];then
-	    # Make didn't exit successfully
-	    echo Error make $i in Source not successful > ../error
-	    cat ../error
-	    exit
-	fi
-    done
-    cd ..
-else
-    echo 'Error Source directory not found' > error
-    cat error
-    exit
-fi
-

=== removed file 'Template/bin/cpall'
--- Template/bin/cpall	2010-10-30 03:26:37 +0000
+++ Template/bin/cpall	1970-01-01 00:00:00 +0000
@@ -1,4 +0,0 @@
-#!/bin/bash
-for i in `cat subproc.mg` ; do
-   cp -p $1 $i/$1
-done

=== removed file 'Template/bin/gen_crossxhtml-pl'
--- Template/bin/gen_crossxhtml-pl	2011-04-21 15:56:29 +0000
+++ Template/bin/gen_crossxhtml-pl	1970-01-01 00:00:00 +0000
@@ -1,399 +0,0 @@
-#!/usr/bin/perl -w
-
-$server = $ENV{'SERVER_NAME'};
-
-
-
-&writeHtmlProcessPage;
-
-
-sub writeHtmlProcessPage {
- 
-#
-#    Make sure we are in the main directory, not bin
-#
-  if (-d "bin") {
-#       print "We are in main directory \n";
-#       print qx(pwd); 	
-  }
-  else{
-#      print "Not in correct directory \n";
-      chdir("../");
-      if (-d "bin") {
-#	  print "We are in main directory \n";
-      }
-      else{
-	  print "Error gen_crossxhtml-pl must be run from main or bin directory \n";
-	  exit;
-      }	  
-  }  
-  $directory=qx(pwd);
-
-#
-# name of the output file
-#
-  $htfile =  "HTML/crossx" . "\.html"; 
-  open(PAGE,"> $htfile");
-#
-# The name of the current run
-#
-  if($#ARGV < 0) {
-      $RunName="Web";
-  } else {
-      $RunName="$ARGV[0]";
-  }
-
-#
-# Write out the current running situation
-#
-
-#--------
-# Getting info
-#--------
-  if (-e "SubProcesses/procdef_mg5.dat"){
-    open (INCARD,"SubProcesses/procdef_mg5.dat") || die "Error reading file proc_card.dat";
-    @incard=<INCARD>;
-    close (INCARD);
-  } else {
-    open (INCARD,"Cards/proc_card.dat") || exit;
-    @incard=<INCARD>;
-    close (INCARD);
-  }
-
-  # process
-  $listpos = 0;