Merge lp:~mg5core2/mg5amcnlo/2.5.5 into lp:mg5amcnlo/lts

[Olivier Mattelaer]

Hi guys,

I have found (and fix) a quite important bug today.
It is related to the LO gridpack where the grid generated during the setup of the gridpack were not saved inside the gridpack (killing the efficiency of such gridpack).

I therefore plan to freeze this version on Monday and release it asap after that.
As soon as the release will be done, we will need to send an email to CMS and Atlas such that they are aware of such problem.

Cheers,

Olivier

[Rikkert Frederix (frederix)]

Hi Olivier,

Hasn't this always been the case? I thought that the gridpacks became too large when all the grids are saved. I thought that the only really important information to safe was the relative contributions of all the channels, so that one exactly knows how many events each channel should give.

Best,
Rikkert

[Olivier Mattelaer (olivier-mattelaer)]

Hi Rik,

> Hasn't this always been the case?

That was my first reaction when I realised that they were missing. So I checked that point explicitly and indeed this is something new (since 2.4.3).

> I thought that the only really important information to safe was the relative contributions of all the channels, so that one exactly knows how many events each channel should give.

I agree that this is the crucial information. Now discarding the grid is wasting a lot of information.
I also have to double check this but in the last run I did, it was written that the grid was kept fixed.
If this is the case, starting from the default grid sounds a bad idea.

Cheers,

Olivier

> On 20 May 2017, at 02:47, Rikkert Frederix <email address hidden> wrote:
>
> Hi Olivier,
>
> Hasn't this always been the case? I thought that the gridpacks became too large when all the grids are saved. I thought that the only really important information to safe was the relative contributions of all the channels, so that one exactly knows how many events each channel should give.
>
> Best,
> Rikkert
>
> --
> https://code.launchpad.net/~mg5core2/mg5amcnlo/2.5.5/+merge/324347
> Your team mg5core2 is subscribed to branch lp:~mg5core2/mg5amcnlo/2.5.5.

[lp:~mg5core2/mg5amcnlo/2.5.5]

297. By olivier-mattelaer on 2017-05-22

fixing problem with complex mass scheme and single core running of Delphes

298. By olivier-mattelaer on 2017-05-22

add an acceptance test related to complex-mass-scheme and madspin

299. By olivier-mattelaer on 2017-05-22

fixing one acceptance test

300. By olivier-mattelaer on 2017-05-22

update Update Notes

[Rikkert Frederix (frederix)]

Hi Olivier,

I saw that in your fix of the grid-pack creation, you also changed all(?) the default Delphes cards. If this is indeed intended, there should be a mention of this in the UpdateNotes.txt.

For the rest: approve.

Best,
Rikkert

[lp:~mg5core2/mg5amcnlo/2.5.5]

301. By Rikkert Frederix on 2017-05-23

Fix in FxFx merging in case there are diagrams with 1->3 decays. See
question #631090.

302. By Valentin Hirschi on 2017-05-24

1. Fixed a subtle side-effect of having removed the recycling of the FD structure
   attached to the loops, which was necessary to solve a fermion flow problem.
   The problem was that as a result, some of the external wavefunctions were not
   recycled as they should have. This could be fixed by creating the wavefunctions
   of the FD structures with *all* the different lorentz structure *at once* instead
   of progressively adding them.

[Valentin Hirschi (valentin-hirschi)]

Hi Olivier, quick last minute review:

MadSpin/interface_madspin.py:
two left-over misc.sprint

This change seems a bit hacky, is it really safe?
For example, what could happen if the first operator used is exponentiation, i.e. ** ?

=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py 2017-03-27 08:37:47 +0000
+++ aloha/aloha_writers.py 2017-05-24 17:51:30 +0000
@@ -802,7 +802,13 @@
             for ind in numerator.listindices():
                 formatted = self.write_obj(numerator.get_rep(ind))
                 if formatted.startswith(('+','-')):
- formatted = '(%s)*%s' % tuple(formatted.split('*',1))
+ if '*' in formatted:
+ formatted = '(%s)*%s' % tuple(formatted.split('*',1))
+ else:
+ if formatted.startswith('+'):
+ formatted = formatted[1:]
+ else:
+ formatted = '(-1)*%s' % formatted[1:]
                 to_order[self.pass_to_HELAS(ind)] = \
                         ' %s(%d)= %s%s\n' % (self.outname, self.pass_to_HELAS(ind)+1,
                         coeff, formatted)

Going through the other changes nothing struck me.
However, I don't know well enough MadFKS to be able to assess the impact of Rik's clean-up of the code.

So it's all ok for me.

review: Approve

[Olivier Mattelaer (olivier-mattelaer)]

Hi,

Let assume that formatted returned "+23**2"
Then it will be transformed in (+23)**2
So this should not be a problem.

Cheers,

Olivier

PS: actually exponent are always associated to paranthesis so the problem will never occur.

> On 24 May 2017, at 23:26, Valentin Hirschi <email address hidden> wrote:
>
> Review: Approve
>
> Hi Olivier, quick last minute review:
>
> MadSpin/interface_madspin.py:
> two left-over misc.sprint
>
> This change seems a bit hacky, is it really safe?
> For example, what could happen if the first operator used is exponentiation, i.e. ** ?
>
> === modified file 'aloha/aloha_writers.py'
> --- aloha/aloha_writers.py 2017-03-27 08:37:47 +0000
> +++ aloha/aloha_writers.py 2017-05-24 17:51:30 +0000
> @@ -802,7 +802,13 @@
> for ind in numerator.listindices():
> formatted = self.write_obj(numerator.get_rep(ind))
> if formatted.startswith(('+','-')):
> - formatted = '(%s)*%s' % tuple(formatted.split('*',1))
> + if '*' in formatted:
> + formatted = '(%s)*%s' % tuple(formatted.split('*',1))
> + else:
> + if formatted.startswith('+'):
> + formatted = formatted[1:]
> + else:
> + formatted = '(-1)*%s' % formatted[1:]
> to_order[self.pass_to_HELAS(ind)] = \
> ' %s(%d)= %s%s\n' % (self.outname, self.pass_to_HELAS(ind)+1,
> coeff, formatted)
>
> Going through the other changes nothing struck me.
> However, I don't know well enough MadFKS to be able to assess the impact of Rik's clean-up of the code.
>
> So it's all ok for me.
>
>
>
>
> --
> https://code.launchpad.net/~mg5core2/mg5amcnlo/2.5.5/+merge/324347
> Your team mg5core2 is subscribed to branch lp:~mg5core2/mg5amcnlo/2.5.5.

[lp:~mg5core2/mg5amcnlo/2.5.5]

303. By olivier-mattelaer on 2017-05-24

fixing last tests + ma5 related bug when latex is not available

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2017-03-17 19:58:41 +0000
+++ MadSpin/decay.py	2017-05-24 22:17:46 +0000
@@ -1638,6 +1638,12 @@
                  ' '.join('--%s=%s' % (key,value) for (key,value) in opts.items()
                         if key not in ['model', 'force', 'particles'] and value))
         cmd.exec_cmd('import model %s' % model.get('modelpath+restriction'))
+
+        #pattern for checking complex mass scheme.
+        has_cms = re.compile(r'''set\s+complex_mass_scheme\s*(True|T|1|true|$|;)''', re.M)
+        force_CMS =  has_cms.search(self.banner['mg5proccard'])
+        if force_CMS:
+            cmd.exec_cmd('set complex_mass_scheme')
 #        cmd._curr_model = model
 #        cmd._curr_fortran_model = helas_call_writers.FortranUFOHelasCallWriter(model)
         cmd.exec_cmd(line)
@@ -3210,7 +3216,7 @@
                     logger.warning( 'no events for %s' % decay_tag)
                     continue
                 weights.sort(reverse=True)
-                assert weights[0] >= weights[1]
+                assert len(weights) == 1 or weights[0] >= weights[1]
                 ave_weight, std_weight = decay_tools.get_mean_sd(weights)
                 base_max_weight = 1.05 * (ave_weight+self.options['nb_sigma']*std_weight)
 

=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2017-03-23 16:11:59 +0000
+++ MadSpin/interface_madspin.py	2017-05-24 22:17:46 +0000
@@ -1412,6 +1412,8 @@
             info = self.generate_all.all_me[tag]
         except:
             misc.sprint(self.generate_all.all_me)
+            misc.sprint(production)
+            misc.sprint(decays)
             raise
         import copy
         

=== modified file 'Template/Common/Cards/delphes_card_ATLAS.dat'
--- Template/Common/Cards/delphes_card_ATLAS.dat	2016-08-29 22:10:09 +0000
+++ Template/Common/Cards/delphes_card_ATLAS.dat	2017-05-24 22:17:46 +0000
@@ -3,7 +3,7 @@
 #######################################
 
 set ExecutionPath {
-  ParticlePropagator
+ ParticlePropagator
 
   ChargedHadronTrackingEfficiency
   ElectronTrackingEfficiency
@@ -14,9 +14,14 @@
   MuonMomentumSmearing
 
   TrackMerger
+
+  ECal
+  HCal
+
   Calorimeter
   EFlowMerger
-
+  EFlowFilter
+  
   PhotonEfficiency
   PhotonIsolation
 
@@ -24,6 +29,8 @@
   ElectronEfficiency
   ElectronIsolation
 
+  ChargedHadronFilter
+
   MuonEfficiency
   MuonIsolation
 
@@ -32,11 +39,11 @@
   NeutrinoFilter
   GenJetFinder
   GenMissingET
-  
+
   FastJetFinder
 
   JetEnergyScale
-  
+
   JetFlavorAssociation
 
   BTagging
@@ -140,9 +147,9 @@
   # set ResolutionFormula {resolution formula as a function of eta and pt}
 
   # resolution formula for charged hadrons
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
 }
 
 ###################################
@@ -188,26 +195,23 @@
   set OutputArray tracks
 }
 
-#############
-# Calorimeter
-#############
-
-module Calorimeter Calorimeter {
+
+#############
+#   ECAL
+#############
+
+module SimpleCalorimeter ECal {
   set ParticleInputArray ParticlePropagator/stableParticles
   set TrackInputArray TrackMerger/tracks
 
-  set TowerOutputArray towers
-  set PhotonOutputArray photons
-
+  set TowerOutputArray ecalTowers
   set EFlowTrackOutputArray eflowTracks
-  set EFlowPhotonOutputArray eflowPhotons
-  set EFlowNeutralHadronOutputArray eflowNeutralHadrons
-
-  set ECalEnergyMin 0.5
-  set HCalEnergyMin 1.0
-
-  set ECalEnergySignificanceMin 1.0
-  set HCalEnergySignificanceMin 1.0
+  set EFlowTowerOutputArray eflowPhotons
+
+  set IsEcal true
+
+  set EnergyMin 0.5
+  set EnergySignificanceMin 2.0
 
   set SmearTowerCenter true
 
@@ -217,6 +221,107 @@
   # each list starts with the lower edge of the first tower
   # the list ends with the higher edged of the last tower
 
+  # assume 0.02 x 0.02 resolution in eta,phi in the barrel |eta| < 1.5
+
+  set PhiBins {}
+  for {set i -180} {$i <= 180} {incr i} {
+    add PhiBins [expr {$i * $pi/180.0}]
+  }
+
+  # 0.02 unit in eta up to eta = 1.5 (barrel)
+  for {set i -85} {$i <= 86} {incr i} {
+    set eta [expr {$i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  # assume 0.02 x 0.02 resolution in eta,phi in the endcaps 1.5 < |eta| < 3.0
+  set PhiBins {}
+  for {set i -180} {$i <= 180} {incr i} {
+    add PhiBins [expr {$i * $pi/180.0}]
+  }
+
+  # 0.02 unit in eta up to eta = 3
+  for {set i 1} {$i <= 84} {incr i} {
+    set eta [expr { -2.958 + $i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  for {set i 1} {$i <= 84} {incr i} {
+    set eta [expr { 1.4964 + $i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  # take present CMS granularity for HF
+
+  # 0.175 x (0.175 - 0.35) resolution in eta,phi in the HF 3.0 < |eta| < 5.0
+  set PhiBins {}
+  for {set i -18} {$i <= 18} {incr i} {
+    add PhiBins [expr {$i * $pi/18.0}]
+  }
+
+  foreach eta {-5 -4.7 -4.525 -4.35 -4.175 -4 -3.825 -3.65 -3.475 -3.3 -3.125 -2.958 3.125 3.3 3.475 3.65 3.825 4 4.175 4.35 4.525 4.7 5} {
+    add EtaPhiBins $eta $PhiBins
+  }
+
+
+  add EnergyFraction {0} {0.0}
+  # energy fractions for e, gamma and pi0
+  add EnergyFraction {11} {1.0}
+  add EnergyFraction {22} {1.0}
+  add EnergyFraction {111} {1.0}
+  # energy fractions for muon, neutrinos and neutralinos
+  add EnergyFraction {12} {0.0}
+  add EnergyFraction {13} {0.0}
+  add EnergyFraction {14} {0.0}
+  add EnergyFraction {16} {0.0}
+  add EnergyFraction {1000022} {0.0}
+  add EnergyFraction {1000023} {0.0}
+  add EnergyFraction {1000025} {0.0}
+  add EnergyFraction {1000035} {0.0}
+  add EnergyFraction {1000045} {0.0}
+  # energy fractions for K0short and Lambda
+  add EnergyFraction {310} {0.3}
+  add EnergyFraction {3122} {0.3}
+
+  # set ResolutionFormula {resolution formula as a function of eta and energy}
+
+  # set ECalResolutionFormula {resolution formula as a function of eta and energy}
+  # http://arxiv.org/pdf/physics/0608012v1 jinst8_08_s08003
+  # http://villaolmo.mib.infn.it/ICATPP9th_2005/Calorimetry/Schram.p.pdf
+  # http://www.physics.utoronto.ca/~krieger/procs/ComoProceedings.pdf
+  set ResolutionFormula {                      (abs(eta) <= 3.2) * sqrt(energy^2*0.0017^2 + energy*0.101^2) +
+                             (abs(eta) > 3.2 && abs(eta) <= 4.9) * sqrt(energy^2*0.0350^2 + energy*0.285^2)}
+
+
+}
+
+
+
+#############
+#   HCAL
+#############
+
+module SimpleCalorimeter HCal {
+  set ParticleInputArray ParticlePropagator/stableParticles
+  set TrackInputArray ECal/eflowTracks
+
+  set TowerOutputArray hcalTowers
+  set EFlowTrackOutputArray eflowTracks
+  set EFlowTowerOutputArray eflowNeutralHadrons
+
+  set IsEcal false
+
+  set EnergyMin 1.0
+  set EnergySignificanceMin 2.0
+
+  set SmearTowerCenter true
+
+ set pi [expr {acos(-1)}]
+
+  # lists of the edges of each tower in eta and phi
+  # each list starts with the lower edge of the first tower
+  # the list ends with the higher edged of the last tower
+
   # 10 degrees towers
   set PhiBins {}
   for {set i -18} {$i <= 18} {incr i} {
@@ -236,50 +341,99 @@
   }
 
   # default energy fractions {abs(PDG code)} {Fecal Fhcal}
-  add EnergyFraction {0} {0.0 1.0}
+  add EnergyFraction {0} {1.0}
   # energy fractions for e, gamma and pi0
-  add EnergyFraction {11} {1.0 0.0}
-  add EnergyFraction {22} {1.0 0.0}
-  add EnergyFraction {111} {1.0 0.0}
+  add EnergyFraction {11} {0.0}
+  add EnergyFraction {22} {0.0}
+  add EnergyFraction {111} {0.0}
   # energy fractions for muon, neutrinos and neutralinos
-  add EnergyFraction {12} {0.0 0.0}
-  add EnergyFraction {13} {0.0 0.0}
-  add EnergyFraction {14} {0.0 0.0}
-  add EnergyFraction {16} {0.0 0.0}
-  add EnergyFraction {1000022} {0.0 0.0}
-  add EnergyFraction {1000023} {0.0 0.0}
-  add EnergyFraction {1000025} {0.0 0.0}
-  add EnergyFraction {1000035} {0.0 0.0}
-  add EnergyFraction {1000045} {0.0 0.0}
+  add EnergyFraction {12} {0.0}
+  add EnergyFraction {13} {0.0}
+  add EnergyFraction {14} {0.0}
+  add EnergyFraction {16} {0.0}
+  add EnergyFraction {1000022} {0.0}
+  add EnergyFraction {1000023} {0.0}
+  add EnergyFraction {1000025} {0.0}
+  add EnergyFraction {1000035} {0.0}
+  add EnergyFraction {1000045} {0.0}
   # energy fractions for K0short and Lambda
-  add EnergyFraction {310} {0.3 0.7}
-  add EnergyFraction {3122} {0.3 0.7}
+  add EnergyFraction {310} {0.7}
+  add EnergyFraction {3122} {0.7}
 
-  # set ECalResolutionFormula {resolution formula as a function of eta and energy}
-  # http://arxiv.org/pdf/physics/0608012v1 jinst8_08_s08003
+  # http://arxiv.org/pdf/hep-ex/0004009v1
   # http://villaolmo.mib.infn.it/ICATPP9th_2005/Calorimetry/Schram.p.pdf
-  # http://www.physics.utoronto.ca/~krieger/procs/ComoProceedings.pdf
-  set ECalResolutionFormula {                  (abs(eta) <= 3.2) * sqrt(energy^2*0.0017^2 + energy*0.101^2) +
-                             (abs(eta) > 3.2 && abs(eta) <= 4.9) * sqrt(energy^2*0.0350^2 + energy*0.285^2)}
-
   # set HCalResolutionFormula {resolution formula as a function of eta and energy}
-  # http://arxiv.org/pdf/hep-ex/0004009v1
-  # http://villaolmo.mib.infn.it/ICATPP9th_2005/Calorimetry/Schram.p.pdf
-  set HCalResolutionFormula {                  (abs(eta) <= 1.7) * sqrt(energy^2*0.0302^2 + energy*0.5205^2 + 1.59^2) +
+  set ResolutionFormula {                      (abs(eta) <= 1.7) * sqrt(energy^2*0.0302^2 + energy*0.5205^2 + 1.59^2) +
                              (abs(eta) > 1.7 && abs(eta) <= 3.2) * sqrt(energy^2*0.0500^2 + energy*0.706^2) +
                              (abs(eta) > 3.2 && abs(eta) <= 4.9) * sqrt(energy^2*0.09420^2 + energy*1.00^2)}
 }
 
+
+#################
+# Electron filter
+#################
+
+module PdgCodeFilter ElectronFilter {
+  set InputArray HCal/eflowTracks
+  set OutputArray electrons
+  set Invert true
+  add PdgCode {11}
+  add PdgCode {-11}
+}
+
+######################
+# ChargedHadronFilter
+######################
+
+module PdgCodeFilter ChargedHadronFilter {
+  set InputArray HCal/eflowTracks
+  set OutputArray chargedHadrons
+  
+  add PdgCode {11}
+  add PdgCode {-11}
+  add PdgCode {13}
+  add PdgCode {-13}
+}
+
+
+
+###################################################
+# Tower Merger (in case not using e-flow algorithm)
+###################################################
+
+module Merger Calorimeter {
+# add InputArray InputArray
+  add InputArray ECal/ecalTowers
+  add InputArray HCal/hcalTowers
+  set OutputArray towers
+}
+
+
+
 ####################
 # Energy flow merger
 ####################
 
 module Merger EFlowMerger {
 # add InputArray InputArray
-  add InputArray Calorimeter/eflowTracks
-  add InputArray Calorimeter/eflowPhotons
-  add InputArray Calorimeter/eflowNeutralHadrons
-  set OutputArray eflow
+  add InputArray HCal/eflowTracks
+  add InputArray ECal/eflowPhotons
+  add InputArray HCal/eflowNeutralHadrons
+  set OutputArray eflow
+}
+
+######################
+# EFlowFilter
+######################
+
+module PdgCodeFilter EFlowFilter {
+  set InputArray EFlowMerger/eflow
+  set OutputArray eflow
+  
+  add PdgCode {11}
+  add PdgCode {-11}
+  add PdgCode {13}
+  add PdgCode {-13}
 }
 
 ###################
@@ -287,7 +441,7 @@
 ###################
 
 module Efficiency PhotonEfficiency {
-  set InputArray Calorimeter/eflowPhotons
+  set InputArray ECal/eflowPhotons
   set OutputArray photons
 
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
@@ -305,7 +459,7 @@
 
 module Isolation PhotonIsolation {
   set CandidateInputArray PhotonEfficiency/photons
-  set IsolationInputArray EFlowMerger/eflow
+  set IsolationInputArray EFlowFilter/eflow
 
   set OutputArray photons
 
@@ -316,17 +470,6 @@
   set PTRatioMax 0.12
 }
 
-#################
-# Electron filter
-#################
-
-module PdgCodeFilter ElectronFilter {
-  set InputArray Calorimeter/eflowTracks
-  set OutputArray electrons
-  set Invert true
-  add PdgCode {11}
-  add PdgCode {-11}
-}
 
 #####################
 # Electron efficiency
@@ -351,7 +494,7 @@
 
 module Isolation ElectronIsolation {
   set CandidateInputArray ElectronEfficiency/electrons
-  set IsolationInputArray EFlowMerger/eflow
+  set IsolationInputArray EFlowFilter/eflow
 
   set OutputArray electrons
 
@@ -385,7 +528,7 @@
 
 module Isolation MuonIsolation {
   set CandidateInputArray MuonEfficiency/muons
-  set IsolationInputArray EFlowMerger/eflow
+  set IsolationInputArray EFlowFilter/eflow
 
   set OutputArray muons
 
@@ -402,7 +545,7 @@
 
 module Merger MissingET {
 # add InputArray InputArray
-  add InputArray EFlowMerger/eflow
+  add InputArray Calorimeter/towers
   set MomentumOutputArray momentum
 }
 
@@ -502,12 +645,12 @@
 ########################
 
 module JetFlavorAssociation JetFlavorAssociation {
-  
+
   set PartonInputArray Delphes/partons
   set ParticleInputArray Delphes/allParticles
   set ParticleLHEFInputArray Delphes/allParticlesLHEF
   set JetInputArray JetEnergyScale/jets
-  
+
   set DeltaR 0.5
   set PartonPTMin 1.0
   set PartonEtaMax 2.5
@@ -544,7 +687,7 @@
 #############
 
 module TrackCountingTauTagging TauTagging {
- 
+
   set ParticleInputArray Delphes/allParticles
   set PartonInputArray Delphes/partons
   set TrackInputArray TrackMerger/tracks
@@ -554,19 +697,19 @@
   set DeltaRTrack 0.2
 
   set TrackPTMin 1.0
- 
+
   set TauPTMin 1.0
   set TauEtaMax 2.5
 
-  # instructions: {n-prongs} {eff} 
-  
+  # instructions: {n-prongs} {eff}
+
   # 1 - one prong efficiency
   # 2 - two or more efficiency
   # -1 - one prong mistag rate
   # -2 - two or more mistag rate
- 
+
   set BitNumber 0
- 
+
   # taken from ATL-PHYS-PUB-2015-045 (medium working point)
   add EfficiencyFormula {1} {0.70}
   add EfficiencyFormula {2} {0.60}
@@ -603,9 +746,9 @@
   add Branch TrackMerger/tracks Track Track
   add Branch Calorimeter/towers Tower Tower
 
-  add Branch Calorimeter/eflowTracks EFlowTrack Track
-  add Branch Calorimeter/eflowPhotons EFlowPhoton Tower
-  add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
+  add Branch HCal/eflowTracks EFlowTrack Track
+  add Branch ECal/eflowPhotons EFlowPhoton Tower
+  add Branch HCal/eflowNeutralHadrons EFlowNeutralHadron Tower
 
   add Branch GenJetFinder/jets GenJet Jet
   add Branch GenMissingET/momentum GenMissingET MissingET

=== modified file 'Template/Common/Cards/delphes_card_CMS.dat'
--- Template/Common/Cards/delphes_card_CMS.dat	2016-08-29 22:10:09 +0000
+++ Template/Common/Cards/delphes_card_CMS.dat	2017-05-24 22:17:46 +0000
@@ -14,9 +14,14 @@
   MuonMomentumSmearing
 
   TrackMerger
+ 
+  ECal
+  HCal
+ 
   Calorimeter
   EFlowMerger
-
+  EFlowFilter
+  
   PhotonEfficiency
   PhotonIsolation
 
@@ -24,6 +29,8 @@
   ElectronEfficiency
   ElectronIsolation
 
+  ChargedHadronFilter
+
   MuonEfficiency
   MuonIsolation
 
@@ -34,6 +41,7 @@
   GenMissingET
   
   FastJetFinder
+  FatJetFinder
 
   JetEnergyScale
 
@@ -123,9 +131,12 @@
   # tracking efficiency formula for muons
   set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
                                            (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) +
-                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e3) * (0.99) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0e3 )               * (0.99 * exp(0.5 - pt*5.0e-4)) +
+
                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e3) * (0.98) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e3)                * (0.98 * exp(0.5 - pt*5.0e-4)) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
 
@@ -141,9 +152,9 @@
 
   # resolution formula for charged hadrons
   # based on arXiv:1405.6569
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
 }
 
 ###################################
@@ -191,26 +202,125 @@
   set OutputArray tracks
 }
 
-#############
-# Calorimeter
-#############
-
-module Calorimeter Calorimeter {
+
+
+#############
+#   ECAL
+#############
+
+module SimpleCalorimeter ECal {
   set ParticleInputArray ParticlePropagator/stableParticles
   set TrackInputArray TrackMerger/tracks
 
-  set TowerOutputArray towers
-  set PhotonOutputArray photons
-
-  set EFlowTrackOutputArray eflowTracks
-  set EFlowPhotonOutputArray eflowPhotons
-  set EFlowNeutralHadronOutputArray eflowNeutralHadrons
-
-  set ECalEnergyMin 0.5
-  set HCalEnergyMin 1.0
-
-  set ECalEnergySignificanceMin 1.0
-  set HCalEnergySignificanceMin 1.0
+  set TowerOutputArray ecalTowers
+  set EFlowTrackOutputArray eflowTracks
+  set EFlowTowerOutputArray eflowPhotons
+
+  set IsEcal true
+
+  set EnergyMin 0.5
+  set EnergySignificanceMin 2.0
+
+  set SmearTowerCenter true
+
+  set pi [expr {acos(-1)}]
+
+  # lists of the edges of each tower in eta and phi
+  # each list starts with the lower edge of the first tower
+  # the list ends with the higher edged of the last tower
+
+  # assume 0.02 x 0.02 resolution in eta,phi in the barrel |eta| < 1.5
+
+  set PhiBins {}
+  for {set i -180} {$i <= 180} {incr i} {
+    add PhiBins [expr {$i * $pi/180.0}]
+  }
+
+  # 0.02 unit in eta up to eta = 1.5 (barrel)
+  for {set i -85} {$i <= 86} {incr i} {
+    set eta [expr {$i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  # assume 0.02 x 0.02 resolution in eta,phi in the endcaps 1.5 < |eta| < 3.0 (HGCAL- ECAL)
+
+  set PhiBins {}
+  for {set i -180} {$i <= 180} {incr i} {
+    add PhiBins [expr {$i * $pi/180.0}]
+  }
+
+  # 0.02 unit in eta up to eta = 3
+  for {set i 1} {$i <= 84} {incr i} {
+    set eta [expr { -2.958 + $i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  for {set i 1} {$i <= 84} {incr i} {
+    set eta [expr { 1.4964 + $i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  # take present CMS granularity for HF
+
+  # 0.175 x (0.175 - 0.35) resolution in eta,phi in the HF 3.0 < |eta| < 5.0
+  set PhiBins {}
+  for {set i -18} {$i <= 18} {incr i} {
+    add PhiBins [expr {$i * $pi/18.0}]
+  }
+
+  foreach eta {-5 -4.7 -4.525 -4.35 -4.175 -4 -3.825 -3.65 -3.475 -3.3 -3.125 -2.958 3.125 3.3 3.475 3.65 3.825 4 4.175 4.35 4.525 4.7 5} {
+    add EtaPhiBins $eta $PhiBins
+  }
+
+
+  add EnergyFraction {0} {0.0}
+  # energy fractions for e, gamma and pi0
+  add EnergyFraction {11} {1.0}
+  add EnergyFraction {22} {1.0}
+  add EnergyFraction {111} {1.0}
+  # energy fractions for muon, neutrinos and neutralinos
+  add EnergyFraction {12} {0.0}
+  add EnergyFraction {13} {0.0}
+  add EnergyFraction {14} {0.0}
+  add EnergyFraction {16} {0.0}
+  add EnergyFraction {1000022} {0.0}
+  add EnergyFraction {1000023} {0.0}
+  add EnergyFraction {1000025} {0.0}
+  add EnergyFraction {1000035} {0.0}
+  add EnergyFraction {1000045} {0.0}
+  # energy fractions for K0short and Lambda
+  add EnergyFraction {310} {0.3}
+  add EnergyFraction {3122} {0.3}
+
+  # set ResolutionFormula {resolution formula as a function of eta and energy}
+
+  # for the ECAL barrel (|eta| < 1.5), see hep-ex/1306.2016 and 1502.02701
+
+  # set ECalResolutionFormula {resolution formula as a function of eta and energy}
+  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
+  set ResolutionFormula {                      (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
+
+}
+
+
+#############
+#   HCAL
+#############
+
+module SimpleCalorimeter HCal {
+  set ParticleInputArray ParticlePropagator/stableParticles
+  set TrackInputArray ECal/eflowTracks
+
+  set TowerOutputArray hcalTowers
+  set EFlowTrackOutputArray eflowTracks
+  set EFlowTowerOutputArray eflowNeutralHadrons
+
+  set IsEcal false
+
+  set EnergyMin 1.0
+  set EnergySignificanceMin 1.0
 
   set SmearTowerCenter true
 
@@ -248,35 +358,71 @@
   }
 
   # default energy fractions {abs(PDG code)} {Fecal Fhcal}
-  add EnergyFraction {0} {0.0 1.0}
+  add EnergyFraction {0} {1.0}
   # energy fractions for e, gamma and pi0
-  add EnergyFraction {11} {1.0 0.0}
-  add EnergyFraction {22} {1.0 0.0}
-  add EnergyFraction {111} {1.0 0.0}
+  add EnergyFraction {11} {0.0}
+  add EnergyFraction {22} {0.0}
+  add EnergyFraction {111} {0.0}
   # energy fractions for muon, neutrinos and neutralinos
-  add EnergyFraction {12} {0.0 0.0}
-  add EnergyFraction {13} {0.0 0.0}
-  add EnergyFraction {14} {0.0 0.0}
-  add EnergyFraction {16} {0.0 0.0}
-  add EnergyFraction {1000022} {0.0 0.0}
-  add EnergyFraction {1000023} {0.0 0.0}
-  add EnergyFraction {1000025} {0.0 0.0}
-  add EnergyFraction {1000035} {0.0 0.0}
-  add EnergyFraction {1000045} {0.0 0.0}
+  add EnergyFraction {12} {0.0}
+  add EnergyFraction {13} {0.0}
+  add EnergyFraction {14} {0.0}
+  add EnergyFraction {16} {0.0}
+  add EnergyFraction {1000022} {0.0}
+  add EnergyFraction {1000023} {0.0}
+  add EnergyFraction {1000025} {0.0}
+  add EnergyFraction {1000035} {0.0}
+  add EnergyFraction {1000045} {0.0}
   # energy fractions for K0short and Lambda
-  add EnergyFraction {310} {0.3 0.7}
-  add EnergyFraction {3122} {0.3 0.7}
-
-  # set ECalResolutionFormula {resolution formula as a function of eta and energy}
-  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
-  set ECalResolutionFormula {                  (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
-                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
-                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
+  add EnergyFraction {310} {0.7}
+  add EnergyFraction {3122} {0.7}
 
   # set HCalResolutionFormula {resolution formula as a function of eta and energy}
-  set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
+  set ResolutionFormula {                      (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
                              (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.130^2 + energy*2.70^2)}
-}
+
+}
+
+
+#################
+# Electron filter
+#################
+
+module PdgCodeFilter ElectronFilter {
+  set InputArray HCal/eflowTracks
+  set OutputArray electrons
+  set Invert true
+  add PdgCode {11}
+  add PdgCode {-11}
+}
+
+######################
+# ChargedHadronFilter
+######################
+
+module PdgCodeFilter ChargedHadronFilter {
+  set InputArray HCal/eflowTracks
+  set OutputArray chargedHadrons
+  
+  add PdgCode {11}
+  add PdgCode {-11}
+  add PdgCode {13}
+  add PdgCode {-13}
+}
+
+
+###################################################
+# Tower Merger (in case not using e-flow algorithm)
+###################################################
+
+module Merger Calorimeter {
+# add InputArray InputArray
+  add InputArray ECal/ecalTowers
+  add InputArray HCal/hcalTowers
+  set OutputArray towers
+}
+
+
 
 ####################
 # Energy flow merger
@@ -284,18 +430,33 @@
 
 module Merger EFlowMerger {
 # add InputArray InputArray
-  add InputArray Calorimeter/eflowTracks
-  add InputArray Calorimeter/eflowPhotons
-  add InputArray Calorimeter/eflowNeutralHadrons
-  set OutputArray eflow
-}
+  add InputArray HCal/eflowTracks
+  add InputArray ECal/eflowPhotons
+  add InputArray HCal/eflowNeutralHadrons
+  set OutputArray eflow
+}
+
+######################
+# EFlowFilter
+######################
+
+module PdgCodeFilter EFlowFilter {
+  set InputArray EFlowMerger/eflow
+  set OutputArray eflow
+  
+  add PdgCode {11}
+  add PdgCode {-11}
+  add PdgCode {13}
+  add PdgCode {-13}
+}
+
 
 ###################
 # Photon efficiency
 ###################
 
 module Efficiency PhotonEfficiency {
-  set InputArray Calorimeter/eflowPhotons
+  set InputArray ECal/eflowPhotons
   set OutputArray photons
 
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
@@ -313,7 +474,7 @@
 
 module Isolation PhotonIsolation {
   set CandidateInputArray PhotonEfficiency/photons
-  set IsolationInputArray EFlowMerger/eflow
+  set IsolationInputArray EFlowFilter/eflow
 
   set OutputArray photons
 
@@ -324,17 +485,6 @@
   set PTRatioMax 0.12
 }
 
-#################
-# Electron filter
-#################
-
-module PdgCodeFilter ElectronFilter {
-  set InputArray Calorimeter/eflowTracks
-  set OutputArray electrons
-  set Invert true
-  add PdgCode {11}
-  add PdgCode {-11}
-}
 
 #####################
 # Electron efficiency
@@ -359,7 +509,7 @@
 
 module Isolation ElectronIsolation {
   set CandidateInputArray ElectronEfficiency/electrons
-  set IsolationInputArray EFlowMerger/eflow
+  set IsolationInputArray EFlowFilter/eflow
 
   set OutputArray electrons
 
@@ -381,11 +531,9 @@
   # set EfficiencyFormula {efficiency as a function of eta and pt}
 
   # efficiency formula for muons
-  set EfficiencyFormula {                                      (pt <= 10.0)               * (0.00) +
-                                           (abs(eta) <= 1.5) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
-                                           (abs(eta) <= 1.5) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
+  set EfficiencyFormula {                                     (pt <= 10.0)                * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)                * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0)                * (0.95) +
                          (abs(eta) > 2.4)                                                 * (0.00)}
 }
 
@@ -395,7 +543,7 @@
 
 module Isolation MuonIsolation {
   set CandidateInputArray MuonEfficiency/muons
-  set IsolationInputArray EFlowMerger/eflow
+  set IsolationInputArray EFlowFilter/eflow
 
   set OutputArray muons
 
@@ -497,6 +645,43 @@
 }
 
 ##################
+# Fat Jet finder
+##################
+
+module FastJetFinder FatJetFinder {
+  set InputArray EFlowMerger/eflow
+
+  set OutputArray jets
+
+  # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
+  set JetAlgorithm 6
+  set ParameterR 0.8
+
+  set ComputeNsubjettiness 1
+  set Beta 1.0
+  set AxisMode 4
+
+  set ComputeTrimming 1
+  set RTrim 0.2
+  set PtFracTrim 0.05
+
+  set ComputePruning 1
+  set ZcutPrun 0.1
+  set RcutPrun 0.5
+  set RPrun 0.8
+
+  set ComputeSoftDrop 1
+  set BetaSoftDrop 0.0
+  set SymmetryCutSoftDrop 0.1
+  set R0SoftDrop 0.8
+
+  set JetPTMin 200.0
+}
+
+
+
+
+##################
 # Jet Energy Scale
 ##################
 
@@ -601,17 +786,20 @@
   add Branch TrackMerger/tracks Track Track
   add Branch Calorimeter/towers Tower Tower
 
-  add Branch Calorimeter/eflowTracks EFlowTrack Track
-  add Branch Calorimeter/eflowPhotons EFlowPhoton Tower
-  add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
+  add Branch HCal/eflowTracks EFlowTrack Track
+  add Branch ECal/eflowPhotons EFlowPhoton Tower
+  add Branch HCal/eflowNeutralHadrons EFlowNeutralHadron Tower
 
   add Branch GenJetFinder/jets GenJet Jet
   add Branch GenMissingET/momentum GenMissingET MissingET
-
+ 
   add Branch UniqueObjectFinder/jets Jet Jet
   add Branch UniqueObjectFinder/electrons Electron Electron
   add Branch UniqueObjectFinder/photons Photon Photon
   add Branch UniqueObjectFinder/muons Muon Muon
+
+  add Branch FatJetFinder/jets FatJet Jet
+
   add Branch MissingET/momentum MissingET MissingET
   add Branch ScalarHT/energy ScalarHT ScalarHT
 }

=== modified file 'Template/Common/Cards/delphes_card_default.dat'
--- Template/Common/Cards/delphes_card_default.dat	2016-08-29 22:10:09 +0000
+++ Template/Common/Cards/delphes_card_default.dat	2017-05-24 22:17:46 +0000
@@ -14,9 +14,14 @@
   MuonMomentumSmearing
 
   TrackMerger
+ 
+  ECal
+  HCal
+ 
   Calorimeter
   EFlowMerger
-
+  EFlowFilter
+  
   PhotonEfficiency
   PhotonIsolation
 
@@ -24,6 +29,8 @@
   ElectronEfficiency
   ElectronIsolation
 
+  ChargedHadronFilter
+
   MuonEfficiency
   MuonIsolation
 
@@ -34,6 +41,7 @@
   GenMissingET
   
   FastJetFinder
+  FatJetFinder
 
   JetEnergyScale
 
@@ -123,9 +131,12 @@
   # tracking efficiency formula for muons
   set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
                                            (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) +
-                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e3) * (0.99) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0e3 )               * (0.99 * exp(0.5 - pt*5.0e-4)) +
+
                          (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e3) * (0.98) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e3)                * (0.98 * exp(0.5 - pt*5.0e-4)) +
                          (abs(eta) > 2.5)                                                  * (0.00)}
 }
 
@@ -141,9 +152,9 @@
 
   # resolution formula for charged hadrons
   # based on arXiv:1405.6569
-  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
-                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.06^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.10^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.25^2 + pt^2*3.1e-3^2)}
 }
 
 ###################################
@@ -191,26 +202,125 @@
   set OutputArray tracks
 }
 
-#############
-# Calorimeter
-#############
-
-module Calorimeter Calorimeter {
+
+
+#############
+#   ECAL
+#############
+
+module SimpleCalorimeter ECal {
   set ParticleInputArray ParticlePropagator/stableParticles
   set TrackInputArray TrackMerger/tracks
 
-  set TowerOutputArray towers
-  set PhotonOutputArray photons
-
-  set EFlowTrackOutputArray eflowTracks
-  set EFlowPhotonOutputArray eflowPhotons
-  set EFlowNeutralHadronOutputArray eflowNeutralHadrons
-
-  set ECalEnergyMin 0.5
-  set HCalEnergyMin 1.0
-
-  set ECalEnergySignificanceMin 1.0
-  set HCalEnergySignificanceMin 1.0
+  set TowerOutputArray ecalTowers
+  set EFlowTrackOutputArray eflowTracks
+  set EFlowTowerOutputArray eflowPhotons
+
+  set IsEcal true
+
+  set EnergyMin 0.5
+  set EnergySignificanceMin 2.0
+
+  set SmearTowerCenter true
+
+  set pi [expr {acos(-1)}]
+
+  # lists of the edges of each tower in eta and phi
+  # each list starts with the lower edge of the first tower
+  # the list ends with the higher edged of the last tower
+
+  # assume 0.02 x 0.02 resolution in eta,phi in the barrel |eta| < 1.5
+
+  set PhiBins {}
+  for {set i -180} {$i <= 180} {incr i} {
+    add PhiBins [expr {$i * $pi/180.0}]
+  }
+
+  # 0.02 unit in eta up to eta = 1.5 (barrel)
+  for {set i -85} {$i <= 86} {incr i} {
+    set eta [expr {$i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  # assume 0.02 x 0.02 resolution in eta,phi in the endcaps 1.5 < |eta| < 3.0 (HGCAL- ECAL)
+
+  set PhiBins {}
+  for {set i -180} {$i <= 180} {incr i} {
+    add PhiBins [expr {$i * $pi/180.0}]
+  }
+
+  # 0.02 unit in eta up to eta = 3
+  for {set i 1} {$i <= 84} {incr i} {
+    set eta [expr { -2.958 + $i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  for {set i 1} {$i <= 84} {incr i} {
+    set eta [expr { 1.4964 + $i * 0.0174}]
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  # take present CMS granularity for HF
+
+  # 0.175 x (0.175 - 0.35) resolution in eta,phi in the HF 3.0 < |eta| < 5.0
+  set PhiBins {}
+  for {set i -18} {$i <= 18} {incr i} {
+    add PhiBins [expr {$i * $pi/18.0}]
+  }
+
+  foreach eta {-5 -4.7 -4.525 -4.35 -4.175 -4 -3.825 -3.65 -3.475 -3.3 -3.125 -2.958 3.125 3.3 3.475 3.65 3.825 4 4.175 4.35 4.525 4.7 5} {
+    add EtaPhiBins $eta $PhiBins
+  }
+
+
+  add EnergyFraction {0} {0.0}
+  # energy fractions for e, gamma and pi0
+  add EnergyFraction {11} {1.0}
+  add EnergyFraction {22} {1.0}
+  add EnergyFraction {111} {1.0}
+  # energy fractions for muon, neutrinos and neutralinos
+  add EnergyFraction {12} {0.0}
+  add EnergyFraction {13} {0.0}
+  add EnergyFraction {14} {0.0}
+  add EnergyFraction {16} {0.0}
+  add EnergyFraction {1000022} {0.0}
+  add EnergyFraction {1000023} {0.0}
+  add EnergyFraction {1000025} {0.0}
+  add EnergyFraction {1000035} {0.0}
+  add EnergyFraction {1000045} {0.0}
+  # energy fractions for K0short and Lambda
+  add EnergyFraction {310} {0.3}
+  add EnergyFraction {3122} {0.3}
+
+  # set ResolutionFormula {resolution formula as a function of eta and energy}
+
+  # for the ECAL barrel (|eta| < 1.5), see hep-ex/1306.2016 and 1502.02701
+
+  # set ECalResolutionFormula {resolution formula as a function of eta and energy}
+  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
+  set ResolutionFormula {                      (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
+
+}
+
+
+#############
+#   HCAL
+#############
+
+module SimpleCalorimeter HCal {
+  set ParticleInputArray ParticlePropagator/stableParticles
+  set TrackInputArray ECal/eflowTracks
+
+  set TowerOutputArray hcalTowers
+  set EFlowTrackOutputArray eflowTracks
+  set EFlowTowerOutputArray eflowNeutralHadrons
+
+  set IsEcal false
+
+  set EnergyMin 1.0
+  set EnergySignificanceMin 1.0
 
   set SmearTowerCenter true
 
@@ -248,35 +358,71 @@
   }
 
   # default energy fractions {abs(PDG code)} {Fecal Fhcal}
-  add EnergyFraction {0} {0.0 1.0}
+  add EnergyFraction {0} {1.0}
   # energy fractions for e, gamma and pi0
-  add EnergyFraction {11} {1.0 0.0}
-  add EnergyFraction {22} {1.0 0.0}
-  add EnergyFraction {111} {1.0 0.0}
+  add EnergyFraction {11} {0.0}
+  add EnergyFraction {22} {0.0}
+  add EnergyFraction {111} {0.0}
   # energy fractions for muon, neutrinos and neutralinos
-  add EnergyFraction {12} {0.0 0.0}
-  add EnergyFraction {13} {0.0 0.0}
-  add EnergyFraction {14} {0.0 0.0}
-  add EnergyFraction {16} {0.0 0.0}
-  add EnergyFraction {1000022} {0.0 0.0}
-  add EnergyFraction {1000023} {0.0 0.0}
-  add EnergyFraction {1000025} {0.0 0.0}
-  add EnergyFraction {1000035} {0.0 0.0}
-  add EnergyFraction {1000045} {0.0 0.0}
+  add EnergyFraction {12} {0.0}
+  add EnergyFraction {13} {0.0}
+  add EnergyFraction {14} {0.0}
+  add EnergyFraction {16} {0.0}
+  add EnergyFraction {1000022} {0.0}
+  add EnergyFraction {1000023} {0.0}
+  add EnergyFraction {1000025} {0.0}
+  add EnergyFraction {1000035} {0.0}
+  add EnergyFraction {1000045} {0.0}
   # energy fractions for K0short and Lambda
-  add EnergyFraction {310} {0.3 0.7}
-  add EnergyFraction {3122} {0.3 0.7}
-
-  # set ECalResolutionFormula {resolution formula as a function of eta and energy}
-  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
-  set ECalResolutionFormula {                  (abs(eta) <= 1.5) * (1+0.64*eta^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
-                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * (2.16 + 5.6*(abs(eta)-2)^2) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
-                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
+  add EnergyFraction {310} {0.7}
+  add EnergyFraction {3122} {0.7}
 
   # set HCalResolutionFormula {resolution formula as a function of eta and energy}
-  set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
+  set ResolutionFormula {                      (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
                              (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.130^2 + energy*2.70^2)}
-}
+
+}
+
+
+#################
+# Electron filter
+#################
+
+module PdgCodeFilter ElectronFilter {
+  set InputArray HCal/eflowTracks
+  set OutputArray electrons
+  set Invert true
+  add PdgCode {11}
+  add PdgCode {-11}
+}
+
+######################
+# ChargedHadronFilter
+######################
+
+module PdgCodeFilter ChargedHadronFilter {
+  set InputArray HCal/eflowTracks
+  set OutputArray chargedHadrons
+  
+  add PdgCode {11}
+  add PdgCode {-11}
+  add PdgCode {13}
+  add PdgCode {-13}
+}
+
+
+###################################################
+# Tower Merger (in case not using e-flow algorithm)
+###################################################
+
+module Merger Calorimeter {
+# add InputArray InputArray
+  add InputArray ECal/ecalTowers
+  add InputArray HCal/hcalTowers
+  set OutputArray towers
+}
+
+
 
 ####################
 # Energy flow merger
@@ -284,18 +430,33 @@
 
 module Merger EFlowMerger {
 # add InputArray InputArray
-  add InputArray Calorimeter/eflowTracks
-  add InputArray Calorimeter/eflowPhotons
-  add InputArray Calorimeter/eflowNeutralHadrons
-  set OutputArray eflow
-}
+  add InputArray HCal/eflowTracks
+  add InputArray ECal/eflowPhotons
+  add InputArray HCal/eflowNeutralHadrons
+  set OutputArray eflow
+}
+
+######################
+# EFlowFilter
+######################
+
+module PdgCodeFilter EFlowFilter {
+  set InputArray EFlowMerger/eflow
+  set OutputArray eflow
+  
+  add PdgCode {11}
+  add PdgCode {-11}
+  add PdgCode {13}
+  add PdgCode {-13}
+}
+
 
 ###################
 # Photon efficiency
 ###################
 
 module Efficiency PhotonEfficiency {
-  set InputArray Calorimeter/eflowPhotons
+  set InputArray ECal/eflowPhotons
   set OutputArray photons
 
   # set EfficiencyFormula {efficiency formula as a function of eta and pt}
@@ -313,7 +474,7 @@
 
 module Isolation PhotonIsolation {
   set CandidateInputArray PhotonEfficiency/photons
-  set IsolationInputArray EFlowMerger/eflow
+  set IsolationInputArray EFlowFilter/eflow
 
   set OutputArray photons
 
@@ -324,17 +485,6 @@
   set PTRatioMax 0.12
 }
 
-#################
-# Electron filter
-#################
-
-module PdgCodeFilter ElectronFilter {
-  set InputArray Calorimeter/eflowTracks
-  set OutputArray electrons
-  set Invert true
-  add PdgCode {11}
-  add PdgCode {-11}
-}
 
 #####################
 # Electron efficiency
@@ -359,7 +509,7 @@
 
 module Isolation ElectronIsolation {
   set CandidateInputArray ElectronEfficiency/electrons
-  set IsolationInputArray EFlowMerger/eflow
+  set IsolationInputArray EFlowFilter/eflow
 
   set OutputArray electrons
 
@@ -381,11 +531,9 @@
   # set EfficiencyFormula {efficiency as a function of eta and pt}
 
   # efficiency formula for muons
-  set EfficiencyFormula {                                      (pt <= 10.0)               * (0.00) +
-                                           (abs(eta) <= 1.5) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
-                                           (abs(eta) <= 1.5) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
-                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
+  set EfficiencyFormula {                                     (pt <= 10.0)                * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)                * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0)                * (0.95) +
                          (abs(eta) > 2.4)                                                 * (0.00)}
 }
 
@@ -395,7 +543,7 @@
 
 module Isolation MuonIsolation {
   set CandidateInputArray MuonEfficiency/muons
-  set IsolationInputArray EFlowMerger/eflow
+  set IsolationInputArray EFlowFilter/eflow
 
   set OutputArray muons
 
@@ -497,6 +645,43 @@
 }
 
 ##################
+# Fat Jet finder
+##################
+
+module FastJetFinder FatJetFinder {
+  set InputArray EFlowMerger/eflow
+
+  set OutputArray jets
+
+  # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
+  set JetAlgorithm 6
+  set ParameterR 0.8
+
+  set ComputeNsubjettiness 1
+  set Beta 1.0
+  set AxisMode 4
+
+  set ComputeTrimming 1
+  set RTrim 0.2
+  set PtFracTrim 0.05
+
+  set ComputePruning 1
+  set ZcutPrun 0.1
+  set RcutPrun 0.5
+  set RPrun 0.8
+
+  set ComputeSoftDrop 1
+  set BetaSoftDrop 0.0
+  set SymmetryCutSoftDrop 0.1
+  set R0SoftDrop 0.8
+
+  set JetPTMin 200.0
+}
+
+
+
+
+##################
 # Jet Energy Scale
 ##################
 
@@ -601,17 +786,20 @@
   add Branch TrackMerger/tracks Track Track
   add Branch Calorimeter/towers Tower Tower
 
-  add Branch Calorimeter/eflowTracks EFlowTrack Track
-  add Branch Calorimeter/eflowPhotons EFlowPhoton Tower
-  add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
+  add Branch HCal/eflowTracks EFlowTrack Track
+  add Branch ECal/eflowPhotons EFlowPhoton Tower
+  add Branch HCal/eflowNeutralHadrons EFlowNeutralHadron Tower
 
   add Branch GenJetFinder/jets GenJet Jet
   add Branch GenMissingET/momentum GenMissingET MissingET
-
+ 
   add Branch UniqueObjectFinder/jets Jet Jet
   add Branch UniqueObjectFinder/electrons Electron Electron
   add Branch UniqueObjectFinder/photons Photon Photon
   add Branch UniqueObjectFinder/muons Muon Muon
+
+  add Branch FatJetFinder/jets FatJet Jet
+
   add Branch MissingET/momentum MissingET MissingET
   add Branch ScalarHT/energy ScalarHT ScalarHT
 }

=== modified file 'Template/LO/bin/internal/store4grid'
--- Template/LO/bin/internal/store4grid	2016-07-29 15:27:30 +0000
+++ Template/LO/bin/internal/store4grid	2017-05-24 22:17:46 +0000
@@ -53,6 +53,11 @@
 	for j in ftn26.gz ; do
 	    if [[ -e $1_$j ]]; then
 		mv -f $1_$j default_$j
+	    elif [[ -e $j ]]; then
+		mv -f $j default_$j
+	    elif [[ -e ftn26 ]]; then
+		gzip ftn26
+		mv -f $j default_$j
 	    fi
 	done
 	cd ../

=== removed file 'Template/NLO/Source/basecode.f'
--- Template/NLO/Source/basecode.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/basecode.f	1970-01-01 00:00:00 +0000
@@ -1,127 +0,0 @@
-      subroutine basecode_test
-      implicit none
-      integer imax
-      parameter (imax = 8)
-      integer icode,iarray(imax),ibase,i,j
-      logical done
-
-      ibase = 3
-c      do i=0,ibase**3-1
-c         call decode(i,iarray,ibase,imax)
-c         call encode(icode,iarray,ibase,imax)
-c         write(*,*) i,icode,"=",(iarray(j),j=1,imax)
-c      enddo
-      icode = 0
-      call decode(icode,iarray,ibase,imax)
-      iarray(2)=1
-      iarray(4)=1
-      iarray(5)=1
-      iarray(7)=1
-      done = .false.
-      write(*,*) (iarray(j),j=1,imax)
-      do while (.not. done)
-         write(*,*) (iarray(j),j=1,imax)
-         call increment_array(iarray,imax,ibase,done)
-      enddo
-      end
-      
-
-      subroutine EnCode(icode,iarray,ibase,imax)
-c******************************************************************************
-c     Turns array of integers (iarray) values range (0,ibase-1) into a single
-c     integer icode. icode = Sum[ iarray(k) * ibase^k]
-c******************************************************************************
-      implicit none
-c
-c     Arguments
-c
-      integer imax           !Number of integers to encode
-      integer icode          !Output encoded value of iarray
-      integer iarray(imax)   !Input values to be encoded
-      integer ibase          !Base for encoding
-
-c
-c     Local
-c     
-      integer i
-c-----
-c  Begin Code
-c-----
-      icode = 0
-      do i = 1, imax
-         if (iarray(i) .ge. 0 .and. iarray(i) .lt. ibase) then 
-            icode = icode + iarray(i)*ibase**(i-1)
-         else
-            write(*,*) 'Error invalid number to be encoded',i,iarray(i)
-         endif
-      enddo
-      end
-
-      subroutine DeCode(icode,iarray,ibase,imax)
-c******************************************************************************
-c     Decodes icode, into base integers used to create it.
-c     integer icode. icode = Sum[ iarray(k) * ibase^k]
-c******************************************************************************
-      implicit none
-c
-c     Arguments
-c
-      integer imax           !Number of integers to encode
-      integer icode          !Input encoded value of iarray
-      integer iarray(imax)   !Output decoded values icode
-      integer ibase          !Base for encoding
-
-c
-c     Local
-c     
-      integer i, jcode
-c-----
-c  Begin Code
-c-----
-      jcode = icode          !create copy for use
-      do i =  imax, 1, -1
-         iarray(i) = 0
-         do while (jcode .ge. ibase**(i-1) .and. iarray(i) .lt. ibase)
-            jcode = jcode-ibase**(i-1)
-            iarray(i)=iarray(i)+1
-         enddo
-      enddo
-      end
-
-      subroutine increment_array(iarray,imax,ibase,done)
-c************************************************************************
-c     Increments iarray     
-c************************************************************************
-      implicit none
-c
-c     Arguments
-c
-      integer imax          !Input, number of elements in iarray
-      integer ibase         !Base for incrementing, 0 is skipped
-      integer iarray(imax)  !Output:Array of values being incremented
-      logical done          !Output:Set when no more incrementing
-c
-c     Local
-c
-      integer i,j
-      logical found
-c-----
-c  Begin Code
-c-----
-      found = .false.
-      i = 1
-      do while (i .le. imax .and. .not. found)
-         if (iarray(i) .eq. 0) then    !don't increment this
-            i=i+1
-         elseif (iarray(i) .lt. ibase-1) then
-            found = .true.
-            iarray(i)=iarray(i)+1
-         else
-            iarray(i)=1
-            i=i+1
-         endif
-      enddo
-      done = .not. found
-      end
-
-

=== removed file 'Template/NLO/Source/combine_events.f'
--- Template/NLO/Source/combine_events.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/combine_events.f	1970-01-01 00:00:00 +0000
@@ -1,728 +0,0 @@
-      program test
-c*****************************************************************
-c     tests traversing directories to find all events
-c****************************************************************
-      implicit none
-c
-c     Constants
-c
-      include 'run_config.inc'
-      integer    maxsubprocesses
-      parameter (maxsubprocesses=9999)
-      integer    cmax_events
-      parameter (cmax_events=500000)
-      integer    sfnum
-      parameter (sfnum=17)   !Unit number for scratch file
-      integer    maxexternal
-      parameter (maxexternal=17)
-      integer maxpara
-      parameter (maxpara=1000)
-c     
-c     Local
-c
-      character*30 subname(maxsubprocesses)
-      integer i,j,m,ns,nreq,ievent
-      integer kevent,revent,iarray(cmax_events)
-      double precision sum, xsec, xerr, goal_wgt,xarray(cmax_events)
-      integer i4,r8,record_length
-      integer iseed
-      real xran1
-      double precision wgt,maxwgt
-      double precision p(0:4,maxexternal)
-      integer ic(7,maxexternal),n
-      double precision scale,aqcd,aqed
-      character*20 param(maxpara),value(maxpara)
-      integer npara,nunwgt
-      double precision xtrunc, min_goal,max_goal
-      logical keep(cmax_events),done
-      integer ntry
-      logical gridrun,gridpack
-
-      character*140 buff
-
-      data iseed/-1/
-c-----
-c  Begin Code
-c-----
-c
-c     Get requested number of events
-c
-      call load_para(npara,param,value)
-      call get_logical(npara,param,value," gridrun ",gridrun,.false.)
-      call get_logical(npara,param,value," gridpack ",gridpack,.false.)
-      if (gridrun.and.gridpack) then
-         call get_integer(npara,param,value," gevents "  ,nreq  ,2000   )
-      else
-         call get_integer(npara,param,value," nevents "  ,nreq  ,10000   )
-      endif
-
-c   Get information for the <init> block
-      call setrun
-
-c      nreq = 10000
-c
-c     Get total cross section
-c
-      xsec = 0d0
-      xerr = 0d0
-      open(unit=15,file='results.dat',status='old',err=21)
-      read(15,*,err=20) xsec,xerr
-      write(*,*) "Results.dat xsec = ",xsec
- 20   close(15)
- 21   if (nreq .gt. 0 .and. xsec .gt. 0) then
-         goal_wgt = xsec/nreq/4d0   !Extra factor of 4 for weighted events
-      else
-         goal_wgt = 0d0    !Write out everything
-      endif
-c
-c     Get list of subprocesses
-c
-      call get_subprocess(subname,ns)
-
-c
-c     Create scratch file to hold events
-c
-      I4 = 4
-      R8 = 8
-      record_length = 4*I4+maxexternal*I4*7+maxexternal*5*R8+3*R8+
-     &   140
-      open(unit=sfnum,access='direct',file='scratch',err=999,
-     &     recl=record_length)
-c
-c     Loop through subprocesses filling up scratch file with events
-c
-      sum=0d0
-      kevent=0
-      revent=0
-      maxwgt=0d0
-      write(*,*) 'SubProcess/Channel     kept   read   xsec '
-      do i=1,ns
-c         write(*,*) 'Subprocess: ',subname(ns)
-         call read_channels(subname(i),sum,kevent,revent,goal_wgt,maxwgt)
-      enddo 
-c
-c     Get Random order for events
-c
-      do i=1,kevent
-         iarray(i)=i
-         xarray(i)=xran1(iseed)
-      enddo
-      call sortO3(xarray,iarray,kevent)
-c
-c     Write out the events in iarray order
-c
-      open(unit=15,file='events.lhe',status='unknown',err=98)
-      call writebanner(15,kevent,sum,maxwgt,xerr)
-      do i=1,kevent
-         read(sfnum,rec=iarray(i)) wgt,n,
-     &        ((ic(m,j),j=1,maxexternal),m=1,7),ievent,
-     &        ((p(m,j),m=0,4),j=1,maxexternal),scale,aqcd,aqed,
-     &     buff
-         call write_event(15,P,wgt,n,ic,ievent,scale,aqcd,aqed,buff)
-      enddo
-      close(15)
-c
-c     Now select unweighted events.
-c
-      goal_wgt = sum/(nreq*1.03)
-      min_goal = goal_wgt/5d0
-      max_goal = goal_wgt*5d0
-      ntry = 1
-c
-c     Loop to refine guess for goal_wgt while keeping xtrunc<0.01
-c
-      done=.false.
-      do while(.not. done)
-         done=.true.
-         nunwgt=0
-         xtrunc=0d0
-         do i=1,kevent
-            read(sfnum,rec=iarray(i)) wgt,n,
-     &           ((ic(m,j),j=1,maxexternal),m=1,7),ievent,
-     &           ((p(m,j),m=0,4),j=1,maxexternal),scale,aqcd,aqed,
-     &        buff
-            if (wgt .gt. goal_wgt*xran1(iseed)) then
-               keep(i) = .true.
-               nunwgt=nunwgt+1
-               if (wgt .gt. goal_wgt) then
-                  xtrunc=xtrunc+wgt-goal_wgt
-               endif
-            else
-               keep(i)=.false.
-            endif
-         enddo
-         if (xtrunc .gt. 0.01d0*sum) then
-            done=.false.
-            min_goal = max(goal_wgt,min_goal)
-            goal_wgt = goal_wgt*1.3d0            
-            write(*,*) 'Iteration ',ntry, ' too large truncation ',xtrunc/sum,nunwgt
-c            write(*,*) min_goal,goal_wgt,max_goal
-         elseif (nunwgt .lt. nreq) then
-            done=.false.
-            max_goal = min(goal_wgt,max_goal)
-            goal_wgt = goal_wgt*0.95d0
-            write(*,*) 'Iteration ',ntry, ' too few events ',xtrunc/sum,nunwgt
-c            write(*,*) min_goal,goal_wgt,max_goal
-            if (goal_wgt .lt. min_goal) then
-               done=.true.
-               write(*,*) 'Failed to find requested number of unweighted events',nreq,nunwgt
-            endif
-         endif
-         ntry=ntry+1
-         if (ntry .gt. 20) done=.true.
-      enddo
-      if (nunwgt .lt. nreq) then
-         write(*,*) 'Unable to get ',nreq,' events. Writing ',nunwgt
-         nreq = nunwgt
-      else
-         write(*,*) 'Found ',nunwgt,' events writing first ',nreq
-      endif
-      write(*,*) 'Unweighting selected ',nreq, ' events.'
-      write(*,'(a,f5.2,a)') 'Truncated ',xtrunc*100./sum, '% of cross section'
-      open(unit=15,file='unweighted_events.lhe',status='unknown',err=99)
-      call writebanner_u(15,nreq,xsec,xtrunc,xerr)
-      ntry = 0
-      do i=1,kevent
-         if (keep(i) .and. ntry .lt. nreq) then
-            read(sfnum,rec=iarray(i)) wgt,n,
-     &           ((ic(m,j),j=1,maxexternal),m=1,7),ievent,
-     &           ((p(m,j),m=0,4),j=1,maxexternal),scale,aqcd,aqed,
-     $        buff
-            call write_event(15,P,xsec/nreq,n,ic,ievent,scale,aqcd,aqed,buff)
-            ntry=ntry+1
-         endif
-      enddo
-      close(15)
-      close(sfnum)
-      return
- 98   write(*,*) 'Error writing events.dat' 
-      return
- 99   write(*,*) 'Error writing unweighted_events.dat' 
-      return
- 999  write(*,*) 'Error opening scratch file'
-      end
-
-
-      subroutine writebanner(lunw,nevent,sum,maxwgt,xerr)
-c**************************************************************************************
-c     Writes out banner information at top of event file
-c**************************************************************************************
-      implicit none
-c
-c     Arguments
-c     
-      integer lunw,nevent
-      double precision sum,maxwgt,xerr
-c
-c     Local
-c
-      integer i,j
-
-c
-c     Les Houches init block (for the <init> info)
-c
-      integer maxpup
-      parameter(maxpup=100)
-      integer idbmup,pdfgup,pdfsup,idwtup,nprup,lprup
-      double precision ebmup,xsecup,xerrup,xmaxup
-      common /heprup/ idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
-     &     idwtup,nprup,xsecup(maxpup),xerrup(maxpup),
-     &     xmaxup(maxpup),lprup(maxpup)
-
-c
-c     Global
-c
-c      double precision etmin(3:nexternal),etamax(3:nexternal)
-c      double precision                    r2min(3:nexternal,3:nexternal)
-c      double precision s_min(nexternal,nexternal)
-c      common/to_cuts/  etmin     ,etamax     , r2min, s_min
-
-c-----
-c  Begin Code
-c-----
-c
-c     gather the info
-c
-c      call setpara('param_card.dat')
-c      call setcuts
-c
-c     write it out
-c
-c      call write_para(lunw)
-c      write(lunw,'(a70)') '##                                                                    '
-c      write(lunw,'(a70)') '##-------------------                                                 '
-c      write(lunw,'(a70)') '## Run-time options                                                   '
-c      write(lunw,'(a70)') '##-------------------                                                 '
-c      write(lunw,'(a70)') '##                                                                    '
-c     write(lunw,'(a70)') '##********************************************************************'     
-c     write(lunw,'(a70)') '## Standard Cuts                                                     *'
-c     write(lunw,'(a70)') '##********************************************************************'    
-c      write(lunw,'(a13,8i8)')   '## Particle  ',(i,i=3,nexternal)
-c      write(lunw,'(a13,8f8.1)') '## Et       >',(etmin(i),i=3,nexternal)
-c      write(lunw,'(a13,8f8.1)') '## Eta      <',(etamax(i),i=3,nexternal)
-c      do j=3,nexternal-1
-c         write(lunw,'(a,i2,a,8f8.1)') '## d R #',j,'  >',(-0.0,i=3,j),
-c     &        (r2min(i,j),i=j+1,nexternal)
-c         do i=j+1,nexternal
-c            r2min(i,j)=r2min(i,j)**2 !Since r2 returns distance squared
-c         enddo
-c      enddo
-c      do j=3,nexternal-1
-c         write(lunw,'(a,i2,a,8f8.1)') '## s min #',j,'>',
-c     &        (s_min(i,j),i=3,nexternal)
-c      enddo
-c      write(lunw,'(a70)') '#********************************************************************'    
-c
-c     Now write out specific information on the event set
-c
-c
-      write(lunw,'(a)') '<MGGenerationInfo>'
-      write(lunw,'(a30,i10)')   '#  Number of Events        :  ',nevent
-      write(lunw,'(a30,e10.5)') '#  Integrated weight (pb)  :  ',sum
-      write(lunw,'(a30,e10.5)') '#  Max wgt                 :  ',maxwgt
-      write(lunw,'(a30,e10.5)') '#  Average wgt             :  ',sum/nevent
-      write(lunw,'(a)') '</MGGenerationInfo>'
-   
-    
-
-C   Write out compulsory init info
-      write(lunw,'(a)') '</header>'
-      write(lunw,'(a)') '<init>'
-      write(lunw,90) (idbmup(i),i=1,2),(ebmup(i),i=1,2),(pdfgup(i),i=1,2),
-     $   (pdfsup(i),i=1,2),2,nprup
-      do i=1,nprup
-         write(lunw,91) xsecup(i),xerr*xsecup(i)/sum,maxwgt,lprup(i) ! FACTOR OF nevts for maxwgt and wgt? error?
-      enddo
-      write(lunw,'(a)') '</init>'
- 90   FORMAT(2i6,2e19.11,2i2,2i6,i2,i3)
- 91   FORMAT(3e19.11,i4)
-      end
-
-
-      subroutine writebanner_u(lunw,nevent,sum,maxwgt,xerr)
-c**************************************************************************************
-c     Writes out banner information at top of event file
-c**************************************************************************************
-      implicit none
-c
-c     Arguments
-c     
-      integer lunw,nevent
-      double precision sum,maxwgt,xerr
-c
-c     Local
-c
-      integer i,j
-c
-c     Les Houches init block (for the <init> info)
-c
-      integer maxpup
-      parameter(maxpup=100)
-      integer idbmup,pdfgup,pdfsup,idwtup,nprup,lprup
-      double precision ebmup,xsecup,xerrup,xmaxup
-      common /heprup/ idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
-     &     idwtup,nprup,xsecup(maxpup),xerrup(maxpup),
-     &     xmaxup(maxpup),lprup(maxpup)
-c
-c     Global
-c
-c      double precision etmin(3:nexternal),etamax(3:nexternal)
-c      double precision                    r2min(3:nexternal,3:nexternal)
-c      double precision s_min(nexternal,nexternal)
-c      common/to_cuts/  etmin     ,etamax     , r2min, s_min
-
-c-----
-c  Begin Code
-c-----
-c
-c     gather the info
-c
-c      call setpara('param_card.dat')
-c      call setcuts
-c
-c     write it out
-c
-c      call write_para(lunw)
-c      write(lunw,'(a70)') '##                                                                    '
-c      write(lunw,'(a70)') '##-------------------                                                 '
-c      write(lunw,'(a70)') '## Run-time options                                                   '
-c      write(lunw,'(a70)') '##-------------------                                                 '
-c      write(lunw,'(a70)') '##                                                                    '
-c     write(lunw,'(a70)') '##********************************************************************'     
-c     write(lunw,'(a70)') '## Standard Cuts                                                     *'
-c     write(lunw,'(a70)') '##********************************************************************'    
-c      write(lunw,'(a13,8i8)')   '## Particle  ',(i,i=3,nexternal)
-c      write(lunw,'(a13,8f8.1)') '## Et       >',(etmin(i),i=3,nexternal)
-c      write(lunw,'(a13,8f8.1)') '## Eta      <',(etamax(i),i=3,nexternal)
-c      do j=3,nexternal-1
-c         write(lunw,'(a,i2,a,8f8.1)') '## d R #',j,'  >',(-0.0,i=3,j),
-c     &        (r2min(i,j),i=j+1,nexternal)
-c         do i=j+1,nexternal
-c            r2min(i,j)=r2min(i,j)**2 !Since r2 returns distance squared
-c         enddo
-c      enddo
-c      do j=3,nexternal-1
-c         write(lunw,'(a,i2,a,8f8.1)') '## s min #',j,'>',
-c     &        (s_min(i,j),i=3,nexternal)
-c      enddo
-c      write(lunw,'(a70)') '##********************************************************************'    
-c
-c     Now write out specific information on the event set
-c
-
-      write(lunw,'(a)') '<MGGenerationInfo>'
-      write(lunw,'(a30,i10)')   '#  Number of Events        :  ',nevent
-      write(lunw,'(a30,e10.5)') '#  Integrated weight (pb)  :  ',sum
-      write(lunw,'(a30,e10.5)') '#  Truncated wgt (pb)      :  ',maxwgt
-      write(lunw,'(a30,e10.5)') '#  Unit wgt                :  ',sum/nevent
-      write(lunw,'(a)') '</MGGenerationInfo>'
-
-C   Write out compulsory init info
-      write(lunw,'(a)') '</header>'
-      write(lunw,'(a)') '<init>'
-      write(lunw,90) (idbmup(i),i=1,2),(ebmup(i),i=1,2),(pdfgup(i),i=1,2),
-     $   (pdfsup(i),i=1,2),3,nprup
-      do i=1,nprup
-         write(lunw,91) xsecup(i),xerr*xsecup(i)/sum,sum/nevent,lprup(i) ! FACTOR OF nevts for maxwgt and wgt? error?
-      enddo
-      write(lunw,'(a)') '</init>'
- 90   FORMAT(2i6,2e19.11,2i2,2i6,i2,i3)
- 91   FORMAT(3e19.11,i4)
-
-      end
-
-
-      subroutine read_channels(dir,sum,kevent,revent,goal_wgt,maxwgt)
-c*****************************************************************
-c     tests traversing directories to find all events
-c****************************************************************
-      implicit none
-c
-c     Constants
-c
-      character*(*) symfile
-      parameter (symfile='symfact.dat')
-c
-c     Arguments
-c
-      character*(30) dir
-      integer kevent,revent
-      double precision sum,goal_wgt,maxwgt
-c
-c     Local
-c
-      integer i,j, k, ip
-      double precision xi
-      character*50 dirname,dname,channame
-c-----
-c  Begin Code
-c-----
-       i = index(dir," ")
-      dname = dir(1:i-1)// "/" // symfile
-      open(unit=35, file=dname ,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,f8.3,a,a)') 'G',xi,'/'
-               endif
-            endif     
-            ip = index(dirname,'/')
-            channame = dname(1:i-1)// "/" //dirname(1:ip)
-            call read_dir_events(channame(1:i+ip),j,kevent,revent,sum,goal_wgt,maxwgt)
-            write(*,'(a,2i8,e10.3)') channame(1:i+ip),kevent,revent,sum
-         endif
- 98   enddo
- 99   close(35)
-      return
-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
-      channame = dirname(1:ip)
-      call read_dir_events(channame,j,kevent,revent,sum,goal_wgt,maxwgt)
-      write(*,'(a30,i8,e10.3)') channame(1:i+ip),kevent,sum
-      return
-      end
-
-      subroutine read_dir_events(channame,nj,kevent,revent,sum,goal_wgt,maxwgt)
-c********************************************************************
-c********************************************************************
-      implicit none
-c
-c     parameters
-c     
-      integer    sfnum
-      parameter (sfnum=17)   !Unit number for scratch file
-      character*(*) scaled_file
-      parameter (scaled_file='events.lhe')
-      integer maxexternal
-      parameter (maxexternal=17)
-      include 'run_config.inc'
-      integer    max_read
-      parameter (max_read = 2000000)
-c
-c     Arguments
-c
-      character*(*) channame
-      integer nj,kevent,revent
-      double precision sum,goal_wgt,maxwgt
-c
-c     Local
-c
-      double precision wgt
-      double precision p(0:4,maxexternal)
-      double precision gsfact
-      real xwgt(max_read),xtot
-      integer i,j,k,m, ic(7,maxexternal),n
-      double precision scale,aqcd,aqed
-      integer ievent,iseed
-      logical done,found
-      character*140 buff
-      character*60 fullname
-c
-c     Les Houches init block (for the <init> info)
-c
-      integer maxpup
-      parameter(maxpup=100)
-      integer idbmup,pdfgup,pdfsup,idwtup,nprup,lprup
-      double precision ebmup,xsecup,xerrup,xmaxup
-      common /heprup/ idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
-     &     idwtup,nprup,xsecup(maxpup),xerrup(maxpup),
-     &     xmaxup(maxpup),lprup(maxpup)
-      data nprup/0/
-c
-c     external
-c
-      real xran1
-c
-c     data
-c
-      data iseed/-1/
-c-----
-c  Begin Code
-c-----     
-      fullname = channame // "gscalefact.dat"
-      gsfact = 1d0
-      open (unit=15,file=fullname,status='old',err=12)
-      read(15,*) gsfact    !Scale factor for grid runs that only use some channels
- 12   close(15)
-      if (gsfact .eq. 0d0) return
-      fullname = channame // scaled_file      
-      open(unit=15,file=fullname, status='old',err=999)
-      done=.false.
-c
-c     Start by initializing all event variables to zero (not really necessary)
-c
-      do j=1,maxexternal
-         do i=1,7
-            ic(i,j)=0
-         enddo
-         do i=0,4
-            p(i,j) = 0d0
-         enddo
-      enddo
-c
-c     Now loop through events
-c
-      do while (.not. done)
-         call read_event(15,P,wgt,n,ic,ievent,scale,aqcd,aqed,buff,done)
-         if (.not. done) then
-            revent = revent+1
-            wgt = wgt*nj*gsfact                 !symmetry factor * grid factor
-            if (wgt .gt. maxwgt) maxwgt=wgt
-            if (wgt .ge. goal_wgt*xran1(iseed)) then
-               kevent=kevent+1
-               if (wgt .lt. goal_wgt) wgt = goal_wgt
-               write(sfnum,rec=kevent) wgt,n,
-     &              ((ic(i,j),j=1,maxexternal),i=1,7),ievent,
-     &              ((p(i,j),i=0,4),j=1,maxexternal),scale,aqcd,aqed,buff
-               sum=sum+wgt
-               found=.false.
-               do i=1,nprup
-                  if(ievent.eq.lprup(i))then
-                     xsecup(i)=xsecup(i)+wgt
-                     found=.true.
-                  endif
-               enddo
-               if(.not.found)then
-                  nprup=nprup+1
-                  lprup(nprup)=ievent
-                  xsecup(nprup)=wgt
-               endif
-            endif
-         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)         
-c      write(*,*) 'Found ',kevent,' events'
-c      write(*,*) 'Integrated weight',sum
-      return
- 999  write(*,*) 'Error opening file ',channame,scaled_file
-
-      end
-
-
-
-      subroutine get_subprocess(subname,ns)
-c*****************************************************************
-c     tests traversing directories to find all events
-c****************************************************************
-      implicit none
-c
-c     Constants
-c
-      character*(*) plist
-      parameter (plist='subproc.mg')
-      integer    maxsubprocesses
-      parameter (maxsubprocesses=999)
-c
-c     Arguments
-c
-      character*30 subname(maxsubprocesses)
-      integer ns
-c-----
-c  Begin Code
-c-----
-      ns = 1
-      open(unit=15, file=plist,status='old',err=99)
-      do while (.true.)
-         read(15,*,err=999,end=999) subname(ns)
-         ns=ns+1
-      enddo
- 99   subname(ns) = './'
-      write(*,*) "Did not find ", plist
-      return
- 999  ns = ns-1
-      write(*,*) "Found ", ns," subprocesses"
-      close(15)
-      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
-         STOP
-      endif
-      xran1=r(j)
-      r(j)=(float(ix1)+float(ix2)*rm2)*rm1
-      return
-      end
-
-
-      subroutine sort2(array,aux1,n)
-      implicit none
-! Arguments
-      integer n
-      integer aux1(n)
-      double precision array(n)
-!  Local Variables
-      integer i,k
-      double precision temp
-      logical done
-
-!-----------
-! Begin Code
-!-----------
-      do i=n-1,1,-1
-         done = .true.
-         do k=1,i
-            if (array(k) .lt. array(k+1)) then
-               temp = array(k)
-               array(k) = array(k+1)
-               array(k+1) = temp
-               temp = aux1(k)
-               aux1(k) = aux1(k+1)
-               aux1(k+1) = temp
-               done = .false.
-            end if
-         end do
-         if (done) return
-      end do
-      end 
-
-      subroutine sortO3(array,aux1,n)
-
-c   O-Sort Version 3, Sorting routine by Erik Oosterwal
-c   http://www.geocities.com/oosterwal/computer/sortroutines.html
-
-      implicit none
-
-! Arguments
-      integer n
-      integer aux1(n)
-      double precision array(n)
-!  Local Variables
-      integer step,i,itemp
-      double precision SngPhi,SngFib
-
-      SngPhi = 0.78             ! Define phi value
-      SngFib = n * SngPhi       ! Set initial real step size
-      step = int(SngFib)        ! set initial integer step size
-
-      do while (step > 0)
-        do i = 1,n-step         ! Set the range of the lower search cells
-          if (array(aux1(i))<array(aux1(i+step))) then ! Compare cells
-            itemp = aux1(i)     !                       \
-            aux1(i) = aux1(i+step) !                     | Swap cells
-            aux1(i+step) = itemp !                      /
-          end if
-        enddo
-        
-        SngFib = SngFib * SngPhi ! Decrease the Real step size
-        Step = Int(SngFib)      ! Set the integer step value
-
-      enddo
-
-      end

=== modified file 'Template/NLO/Source/genps.inc'
--- Template/NLO/Source/genps.inc	2013-01-11 08:13:45 +0000
+++ Template/NLO/Source/genps.inc	2017-05-24 22:17:46 +0000
@@ -1,44 +1,4 @@
-c*************************************************************************
-c	Parameters used by genps
-c*************************************************************************
-c*************************************************************************
-c       Parameters used by genps and dsample, you must recompile
-c	dsample if you change anything below
-c*************************************************************************
       include 'maxparticles.inc'
-      integer    max_tree    , max_dim
-      parameter (max_tree=200, max_dim=max_tree*(max_branch+1))
-      integer    ng     , maxdim     , maxinvar     , maxconfigs
-c      parameter (ng = 96, maxdim = 25, maxinvar= 50 , maxconfigs=10)
-      parameter (ng = 96, maxdim = 3*(max_particles-2)-1, maxinvar= 4*max_particles, maxconfigs=10)
-
       include 'maxconfigs.inc'
-
-        
-      double precision  xgmin,     xgmax
-      parameter        (xgmin=-1d0, xgmax=1d0)
-
-      integer    maxproc,     maxamps
-      parameter (maxproc=100, maxamps=7500)	
-
-c      integer    maxproc
-c      parameter (maxproc=100)	
-
-      integer    maxevents         !Requires about 1K/event
-      parameter (maxevents=100000)  !Maximum # events to write to disk
-
-c*************************************************************************
-c	Parameters used for parrallel running
-c*************************************************************************
-      integer    max_host   ,maxplace    ,maxpoints    ,maxans
-      parameter (max_host=9,maxplace=9,maxpoints=10,maxans=5)
-      integer    maxprb
-      parameter (maxprb = maxconfigs*maxplace*maxpoints)
-      integer    maxfprb
-      parameter (maxfprb = maxinvar*maxplace*maxpoints)
-
-
-
-
-
-
+      integer    maxproc,     maxamps     ,maxflow
+      parameter (maxproc=100, maxamps=7500,maxflow=999)	

=== removed file 'Template/NLO/Source/hbook.inc'
--- Template/NLO/Source/hbook.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/hbook.inc	1970-01-01 00:00:00 +0000
@@ -1,17 +0,0 @@
-C		Internal common blocks for pheno/hbook routines
-C
-C		LABELS(i) = label for histogram i
-C		nhist	  = number of histograms (starts as 0, max is 20)
-C		idnumber(i) = code number to identify histograms in HFILL
-C		pointer(i)  = index of beginning of data for histo # i
-C		single dim(i) = .true. if single variable, .false. if double
-C
-	parameter (nhistmax=20,nhistmax1=nhistmax+1)
-	real data(10000),error(10000),npoints(10000)
-	integer pointer(nhistmax1),id number(nhistmax)
-	logical single dim(nhistmax)
-	character*40 label(nhistmax)
-
-       common /hbooklabel/ label
-       common /hbookarrays/nhist,id number,pointer,
-     $       single dim,data,error,npoints

=== removed file 'Template/NLO/Source/hcurve.f'
--- Template/NLO/Source/hcurve.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/hcurve.f	1970-01-01 00:00:00 +0000
@@ -1,74 +0,0 @@
-C--------------------------------------------
-C
-C		Routine to dump histogram data to a file
-C
-	subroutine hcurve(id,filename)	
-C
-C		Dumps current histogram number id to file 'filename' and
-C		clears histogram id.
-C
-	include 'hbook.inc'
-	character*(*) filename
-        real sum,npts
-
-	if (nhist .eq. 0) return
-	open (unit=69,name=filename,status='unknown')
-	do i = 1, nhist
-		if (id .eq. idnumber(i)) go to 10
-		end do
-	return
-10	continue
-	k = pointer(i)
-	nx = int(data(k)+.1)
-	xmin = data(k+1)
-	xmax = data(k+2)
-	xbinsize = (xmax-xmin)/nx
-	if (single dim(i)) then
-           sum=0
-           npts=0
-           do m=1,nx
-              sum=sum+data(k+2+m)
-              npts=npts+npoints(k+2+m)
-           enddo
-           write (69,300) label(i)(1:labelleng(label(i)))
-           write (69,700) (xmin+(m-.5)*xbinsize,
-     $          data(k+2+m),sqrt(abs(error(k+2+m))),
-     $          npoints(k+2+m)/(npts*sum+1e-23),m=1,nx)
-	else
-		ny = int(data(k+3) + .1)
-		ymin = data(k+4)
-		ymax = data(k+5)
-		ybinsize = (ymax-ymin)/ny
-		write (69,300) label(i)(1:labelleng(label(i)))
-		k = k + 5
-		do n=1,ny
-                   fixed y =  ymin + (n-.5)*ybinsize
-                   write (69,500) (xmin+(m-.5)*xbinsize,fixed y,
-     $                  data(k+m),m=1,nx)
-		   write(69,*) 
-                   k = k + nx
-                end do
-		end if
-	close (unit=69)
-	return
-300	format ('# Histogram ',a)
-400	format (1x,2g15.6)
-500	format (1x,3g15.6)
-700	format (1x,4g15.6)
-	end
-C
-C
-C
-C
-	function labelleng(string)
-	character*(*) string
-
-	do i=len(string),1,-1
-		if (string(i:i) .ne. ' ') then
-			labelleng=i
-			return
-			end if
-		end do
-	labelleng=1
-	return
-	end

=== removed file 'Template/NLO/Source/hfill.f'
--- Template/NLO/Source/hfill.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/hfill.f	1970-01-01 00:00:00 +0000
@@ -1,37 +0,0 @@
-C----------------------------------------------
-C
-C		Routine to add zincrement to a bin in a histogram
-C
-	subroutine hfill(id,x,y,zincrement)
-C
-C		id = integer associated with the histogram
-C		x  = x value to locate bin (real)
-C		y  = y value to locate bin (real) [ignored for 1-dim histo]
-C		zincrement = value to be added to bin specified by (x,y)
-C
-	include 'hbook.inc'
-	data nhist/0/,pointer(1)/1/
-
-	do i=1,nhist
-		if (id number(i) .eq. id) go to 10
-		end do
-	print*,' id number ',id,' does not belong to a current histogram'
-	return
-10	continue
-	k = pointer(i)
-	nx=data(k)+.1
-	ixoff = (x-data(k+1))/(data(k+2)-data(k+1))*data(k)+1
-	if (ixoff .le. 0 .or. ixoff .gt. nx) return
-	if (single dim(i)) then
-		data(k+2+ixoff)=data(k+2+ixoff)+zincrement
-		error(k+2+ixoff)=error(k+2+ixoff)+zincrement**2
-		npoints(k+2+ixoff)=npoints(k+2+ixoff)+1.
-	else
-		ny=data(k+3)+.1
-		iyoff = (y-data(k+4))/(data(k+5)-data(k+4))*data(k+3)+1
-		if (iyoff .le. 0 .or. iyoff .gt. ny) return
-		ioff = nx*(iyoff-1)+ixoff
-		data(k+5+ioff)=data(k+5+ioff)+zincrement
-		end if
-	return
-	end

=== removed file 'Template/NLO/Source/invarients.f'
--- Template/NLO/Source/invarients.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/invarients.f	1970-01-01 00:00:00 +0000
@@ -1,302 +0,0 @@
-      subroutine set_invarients(nfinal,ninvar)
-c***************************************************************************
-c     Calculates all of the invarients for a 2->n process
-c***************************************************************************
-      implicit none
-c     
-c     Constants
-c
-      include 'genps.inc'
-c
-c     Arguments
-c
-      integer nfinal,ninvar
-c
-c     Local
-c     
-      integer ip1,ip2,ipstart,ipstop,np,i
-      integer ncycle
-      character*10 buff
-c
-c     Global
-c
-      integer              imom(maxinvar),ninvarients
-      common/to_invarients/imom          ,ninvarients
-c-----
-c  Begin Code
-c-----
-
-      do i=1,nfinal
-         imom(i)=i
-      enddo
-      ipstart=1
-      ipstop =nfinal
-      np     =nfinal
-c
-c     First do all the s-channel
-c
-      do ncycle=2,nfinal-1
-         do ip1 = ipstart,ipstop-1
-            do ip2=int((real(imom(ip1))/10.-imom(ip1)/10)*10+.1)+1,
-     $           nfinal
-               np=np+1
-               if (np .gt. maxinvar) then
-                  print*,'Sorry too many invarients',np,ip1,ip2,ncycle
-                  stop
-               endif
-               imom(np)=imom(ip1)*10+imom(ip2)
-               if (imom(np) .lt. 10) then
-                  write(buff,'(a2,i1)') 'S?',imom(np)
-               elseif (imom(np) .lt. 100) then
-                  write(buff,'(a2,i2)') 'S?',imom(np)
-               elseif (imom(np) .lt. 1000) then
-                  write(buff,'(a2,i3)') 'S?',imom(np)
-               elseif (imom(np) .lt. 10000) then
-                  write(buff,'(a2,i4)') 'S?',imom(np)
-               elseif (imom(np) .lt. 100000) then
-                  write(buff,'(a2,i5)') 'S?',imom(np)
-               else
-                  write(buff,'(a2,i6)') 'S?',imom(ip1)
-               endif
-c               call hbook1(100+np-nfinal,buff,100,0.,1.,0.)
-c               write(*,'(i4,i6)') np-nfinal,imom(np)
-               write(*,'(i4,a1,a6)') np-nfinal,'=',buff
-               if ((np-nfinal)/7 .eq. real(np-nfinal)/7.) write(*,*)' '
-            enddo
-         enddo
-         ipstart=ipstop+1
-         ipstop = np
-      enddo
-c
-c     Now do the t-channel
-c
-      ipstop = np
-      do ip1 = 1,ipstop
-c         write(*,'(i4,a2,i6)') np-nfinal+ip1,'a-',imom(ip1)
-         if (imom(ip1) .lt. 10) then
-            write(buff,'(a2,i1)') 'T?',imom(ip1)
-         elseif (imom(ip1) .lt. 100) then
-            write(buff,'(a2,i2)') 'T?',imom(ip1)
-         elseif (imom(ip1) .lt. 1000) then
-            write(buff,'(a2,i3)') 'T?',imom(ip1)
-         elseif (imom(ip1) .lt. 10000) then
-            write(buff,'(a2,i4)') 'T?',imom(ip1)
-         elseif (imom(ip1) .lt. 100000) then
-            write(buff,'(a2,i5)') 'T?',imom(ip1)
-         else
-            write(buff,'(a2,i6)') 'T?',imom(ip1)
-         endif
-c         call hbook1(100+np-nfinal+ip1,buff,100,0.,1.,0.)
-c         write(*,*) np-nfinal+ip1,buff
-         write(*,'(i4,a1,a6)') np-nfinal+ip1,'=',buff
-         if ((np-nfinal+ip1)/7 .eq. real(np-nfinal+ip1)/7.) write(*,*)
-      enddo
-      write(*,*)
-      print*,'Particles, Invarients',nfinal,np-nfinal+np
-      ninvarients=np-nfinal+np
-      ninvar=ninvarients
-      if (ninvarients .gt. maxinvar) then
-         print*,'Error too many invarients to map'
-c         stop
-      endif
-      end
-
-
-      subroutine fill_invarients(nfinal,p1,s,xx)
-c***************************************************************************
-c     Calculates all of the invarients for a 2->n process
-c***************************************************************************
-      implicit none
-c     
-c     Constants
-c     
-      include 'genps.inc'
-c
-c     Arguments
-c
-      integer nfinal
-      double precision p1(0:3,nfinal+2),s,xx(55)
-c
-c     Local
-c     
-      integer ip1,ip2,ipstart,ipstop,np,i,j
-      integer imom(maxinvar)
-      integer ncycle
-      character*10 buff
-      double precision p(0:3,maxinvar),ptemp(0:3)
-c
-c     External
-c
-      double precision dot
-      external         dot
-c-----
-c  Begin Code
-c-----
-
-      do i=1,nfinal
-         imom(i) = i
-         do j=0,3
-            p(j,i)=p1(j,i+2)
-         enddo
-c         write(*,'(i3,4f17.8)') i,(p(j,i),j=0,3)
-      enddo
-      ipstart=1
-      ipstop =nfinal
-      np     =nfinal
-c
-c     First do all the s-channel
-c
-      do ncycle=2,nfinal-1
-         do ip1 = ipstart,ipstop-1
-            do ip2=int((real(imom(ip1))/10.-imom(ip1)/10)*10+.1)+1
-     $           ,nfinal
-               np=np+1
-               if (np .gt. maxinvar) then
-                  print*,'Sorry too many invarients',np,ip1,ip2,ncycle
-                  stop
-               endif
-               imom(np)=imom(ip1)*10+imom(ip2)
-               do j=0,3
-                  p(j,np) = p(j,ip1)+p(j,ip2)
-               enddo
-               xx(np-nfinal) = dot(p(0,np),p(0,np))/s
-c               call hfill(100+np-nfinal,
-c     &              real(dot(p(0,np),p(0,np))/s),0.,wgt)
-c               write(*,'(i4,3f20.8)') np-nfinal,
-c     &              real(dot(p(0,np),p(0,np))/s)
-            enddo
-         enddo
-         ipstart=ipstop+1
-         ipstop = np
-      enddo
-c
-c     Now do the t-channel
-c
-      ipstop = np
-      do ip1 = 1,ipstop
-         do j = 0,3
-            ptemp(j)=p1(j,1)-p(j,ip1)
-         enddo
-         xx(np-nfinal+ip1)= .5d0*(dot(ptemp,ptemp)/s+1d0)
-c         call hfill(100+np-nfinal+ip1,real(-dot(ptemp,ptemp)/s),0.,wgt)
-c         write(*,'(i4,3f20.8)') np-nfinal+ip1,
-c     &         real(-dot(ptemp,ptemp)/s)
-      enddo
-      end
-
-
-      subroutine map_invarients(Minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
-c****************************************************************************
-c     Determines mappings for each structure of invarients onto integration
-c     variables.  Input: Ninvar, iforest.  Output: Minvar, ninvar
-c****************************************************************************
-      implicit none
-c
-c     Constants
-c
-      include 'genps.inc'
-c
-c     Arguments
-c
-      integer Minvar(maxdim,lmaxconfigs),nconfigs,ninvar,nexternal,nincoming
-      integer mincfig,maxcfig
-c
-c     Local
-c
-      integer iconfig, jgrid,j, nbranch
-      logical found,tchannel
-c
-c     Global
-c
-      integer              imom(maxinvar),ninvarients
-      common/to_invarients/imom          ,ninvarients
-      integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
-
-c-----
-c  Begin Code
-c----
-
-      nbranch = nexternal-2
-      jgrid=0
-c
-c     
-c     Try simple mapping if nconfigs = 1
-c
-
-      if (nconfigs .eq. 1) then
-         do j=1,3*nbranch-4+2
-            minvar(j,mincfig)=j
-         enddo
-         jgrid=j-1
-      else
-c      if (ep) jgrid=1
-c      if (pp) jgrid=2
-      do iconfig=mincfig,maxcfig
-         tchannel = .false.         
-         do j=1,nbranch-1
-            if (iforest(1,-j,iconfig) .eq. 1) then
-               tchannel=.true.
-            endif
-            jgrid=jgrid+1
-            minvar(j,iconfig) = jgrid
-            if (tchannel .and. j .lt. nbranch-1) then
-               jgrid=jgrid+1            
-               minvar(nbranch-1+2*j,iconfig)=jgrid
-            endif
-         enddo             !Each Branch
-         if (.not. tchannel .and. nincoming.eq.2) then          !Don't need last s-channel
-            jgrid=jgrid-1
-            minvar(nbranch-1,iconfig)=0
-         endif
-c         if (pp) then
-c            jgrid=jgrid+1
-c            minvar(3*nbranch-3,iconfig)=jgrid
-c            jgrid=jgrid+1
-c            minvar(3*nbranch-2,iconfig)=jgrid
-c         elseif (ep) then
-c            jgrid=jgrid+1
-c            minvar(3*nbranch-3,iconfig)=jgrid
-c         endif
-      enddo  !Each configurations
-      endif
-      ninvar = jgrid
-      end
-
-      subroutine sortint(n,ra)
-      integer 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
-
-

=== removed file 'Template/NLO/Source/leshouche.inc'
--- Template/NLO/Source/leshouche.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/leshouche.inc	1970-01-01 00:00:00 +0000
@@ -1,8 +0,0 @@
-      DATA (IDUP(I,1),I=1,5)/-1,2,-11,12,21/
-      DATA (MOTHUP(1,I,  1),I=1, 5)/  0,  0,  1,  1,  1/
-      DATA (MOTHUP(2,I,  1),I=1, 5)/  0,  0,  2,  2,  2/
-      DATA (ICOLUP(1,I,  1),I=1, 5)/  0,502,  0,  0,502/
-      DATA (ICOLUP(2,I,  1),I=1, 5)/501,  0,  0,  0,501/
-      DATA (IDUP(I,2),I=1,5)/-3,4,-11,12,21/
-      DATA (MOTHUP(1,I,  2),I=1, 5)/  0,  0,  1,  1,  1/
-      DATA (MOTHUP(2,I,  2),I=1, 5)/  0,  0,  2,  2,  2/

=== modified file 'Template/NLO/Source/makefile'
--- Template/NLO/Source/makefile	2015-07-14 18:21:29 +0000
+++ Template/NLO/Source/makefile	2017-05-24 22:17:46 +0000
@@ -8,33 +8,19 @@
 IREGIDIR= ./IREGI/src/
 STDHEPDIR= ./StdHEP/
 
-PROCESS= hfill.o matrix.o myamp.o 
-
-HBOOK = hfill.o hcurve.o hbook1.o hbook2.o
-
-GENERIC = $(alfas_functions).o invarients.o hfill.o pawgraphs.o ran1.o rw_events.o rw_routines.o kin_functions.o \
-	  open_file.o basecode.o setrun.o run_printout.o dgauss.o ranmar.o getissud.o
-
-INCLUDEF= coupl.inc genps.inc hbook.inc psample.inc pmass.inc nexternal.inc cluster.inc
-
-BANNER  =  write_banner.o  rw_events.o ranmar.o  kin_functions.o open_file.o  rw_routines.o alfas_functions.o
-
-COMBINE = combine_events.o  rw_events.o ranmar.o  kin_functions.o open_file.o rw_routines.o alfas_functions.o setrun.o
-
-
+PROCESS= matrix.o myamp.o 
+
+GENERIC = $(alfas_functions).o rw_routines.o kin_functions.o setrun.o	\
+	  run_printout.o dgauss.o ranmar.o getissud.o
+
+INCLUDEF = coupl.inc pmass.inc cluster.inc
 
 .f.o: ; $(FC) $(FFLAGS) -c $*.f
 
-all: $(LIBDIR)libdhelas.a $(LIBDIR)libgeneric.a $(LIBDIR)libpdf.a  $(LIBDIR)libmodel.a  $(LIBDIR)libcernlib.a \
-    param_card.inc
+all: $(LIBDIR)libdhelas.a $(LIBDIR)libgeneric.a $(LIBDIR)libpdf.a	\
+    $(LIBDIR)libmodel.a $(LIBDIR)libcernlib.a param_card.inc
 	rm -f PDF/*.o
 
-scaled.dat:
-	combine_events
-
-$(BINDIR)sum_html: sum_html.o 
-	$(FC) $(FFLAGS) -o sum_html sum_html.o   
-	mv sum_html $(BINDIR)
 
 $(LIBDIR)libdhelas.a: DHELAS
 	rm -f $(LIBDIR)libdhelas.a
@@ -62,23 +48,12 @@
 	rm -f $(LIBDIR)libmodel.a
 	cd MODEL; make	
 
-genps.inc: nexternal.inc
-	touch genps.inc
-
 run_card.inc: ../Cards/run_card.dat
 	../bin/aMCatNLO treatcards run
 
 param_card.inc: ../Cards/param_card.dat
 	../bin/aMCatNLO treatcards param
 
-$(BINDIR)gen_ximprove: gen_ximprove.o ranmar.o rw_routines.o open_file.o
-	$(FC) $(FFLAGS) -o gen_ximprove gen_ximprove.o ranmar.o rw_routines.o open_file.o
-	mv gen_ximprove $(BINDIR)
-
-$(BINDIR)combine_events: $(COMBINE)
-	$(FC) $(FFLAGS) -o combine_events $(COMBINE) $(LIBDIR)/libmodel.a $(LIBDIR)/libpdf.a 
-	mv combine_events $(BINDIR)
-
 CutTools: $(LIBDIR)libcts.a
 libcuttools: $(LIBDIR)libcts.a
 
@@ -128,4 +103,3 @@
 	if [ -d $(CUTTOOLSDIR) ]; then cd $(CUTTOOLSDIR); make clean; cd ..; fi
 	if [ -d $(IREGIDIR) ]; then cd $(IREGIDIR); make clean; cd ..; fi
 	if [ -d $(STDHEPDIR) ]; then cd $(STDHEPDIR); make clean; cd ..; fi
-	rm -f $(BINDIR)/combine_events $(BINDIR)/gen_ximprove

=== modified file 'Template/NLO/Source/nexternal.inc'
--- Template/NLO/Source/nexternal.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/nexternal.inc	2017-05-24 22:17:46 +0000
@@ -1,4 +0,0 @@
-      INTEGER    NEXTERNAL
-      PARAMETER (NEXTERNAL=5)
-      INTEGER    NINCOMING
-      PARAMETER (NINCOMING=2)

=== removed file 'Template/NLO/Source/open_file.f'
--- Template/NLO/Source/open_file.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/open_file.f	1970-01-01 00:00:00 +0000
@@ -1,46 +0,0 @@
-      subroutine open_file(lun,filename,fopened)
-c***********************************************************************
-c     opens file input-card.dat in current directory or above
-c***********************************************************************
-      implicit none
-c
-c     Arguments
-c
-      integer lun
-      logical fopened
-      character*(*) filename
-      character*90  tempname
-      integer fine
-      integer i
-      
-c-----
-c  Begin Code
-c-----
-c
-c     first check that we will end in the main directory
-c
-
-c
-c 	  if I have to read a card
-c
-
-      tempname=filename 	 
-      fine=index(tempname,' ') 	 
-      if(fine.eq.0) fine=len(tempname)
-
-      if(index(filename,"_card").gt.0) then
-         tempname='Cards/'//tempname(1:fine)
-         fine=fine+6
-      endif
-      
-      fopened=.false.
-      do i=0,5
-         open(unit=lun,file=tempname,status='old',ERR=30)
-         fopened=.true.
-         exit
- 30      tempname='../'//tempname
-         if (i.eq.5)then
-            write(*,*) 'Warning: file ',filename,' not found'
-         endif
-      enddo
-      end

=== removed file 'Template/NLO/Source/pawgraphs.f'
--- Template/NLO/Source/pawgraphs.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/pawgraphs.f	1970-01-01 00:00:00 +0000
@@ -1,83 +0,0 @@
-      subroutine graph_init
-c*************************************************************************
-c     Set up graphing
-c*************************************************************************
-      implicit none
-c
-c     Local
-c
-      real xmin,xmax
-c
-c     Global
-c
-      real        h(80000)
-      common/pawc/h
-c-----
-c  Begin Code
-c-----
-c      call hlimit(80000)
-c
-c     Total
-c
-c      print*,'Setting up graphs'
-c      call hbook1(1,'s hat',100,0.,500.,0.)
-      end
-
-      subroutine graph_point2(x,y)
-      double precision x,y
-      end
-
-
-      subroutine graph_point(p,dwgt)
-c***************************************************************************
-c     fill historgrams
-c***************************************************************************
-      implicit none
-c
-c     Constants
-c
-      double precision  pi              , to_deg
-      parameter        (pi = 3.1415927d0, to_deg=180d0/pi)
-c
-c     Arguments
-c
-      double precision dwgt
-      REAL*8 P(0:3,7)
-c
-c     Local
-c
-      real*4 wgt
-      real*8 ptot(0:3),maxamp, shat
-      integer i,iconfig, imax
-c
-c     Global
-c
-      include 'run.inc'
-
-c
-c     External
-c
-      double precision dot,et,eta,r2
-c-----
-c  Begin Code
-c-----
-      wgt=dwgt
-c      call hfill(1,real(et(p(0,4))),0.,wgt)
-      end
-
-      subroutine graph_store
-c*************************************************************************
-c     Stores graphs
-c*************************************************************************
-      implicit none
-
-c-----
-c  Begin Code
-c-----
-c      call hcurve(1,'shat.dat')
-c      call hrput(0,'wg.paw','N')
-      end
-
-
-
-

=== removed file 'Template/NLO/Source/psample.inc'
--- Template/NLO/Source/psample.inc	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/psample.inc	1970-01-01 00:00:00 +0000
@@ -1,9 +0,0 @@
-c
-c Global variables used by psample
-c
-      integer                ihost(max_host),npnts,nans,icpu(max_host)
-      integer                 ierror(max_host)
-      character*30                   hostname(max_host),program
-      common /sample_machine/ ihost,icpu,ierror,npnts,nans, 
-     &                        hostname,program
-

=== removed file 'Template/NLO/Source/ran1.f'
--- Template/NLO/Source/ran1.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/ran1.f	1970-01-01 00:00:00 +0000
@@ -1,33 +0,0 @@
-      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
-         STOP
-      endif
-      xran1=r(j)
-      r(j)=(float(ix1)+float(ix2)*rm2)*rm1
-      return
-      end

=== removed file 'Template/NLO/Source/rw_events.f'
--- Template/NLO/Source/rw_events.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/rw_events.f	1970-01-01 00:00:00 +0000
@@ -1,159 +0,0 @@
-      subroutine read_event(lun,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
-c********************************************************************
-c     Reads one event from data file #lun
-c     ic(*,1) = Particle ID
-c     ic(*,2) = Mothup(1)
-c     ic(*,3) = Mothup(2)
-c     ic(*,4) = ICOLUP(1)
-c     ic(*,5) = ICOLUP(2)
-c     ic(*,6) = ISTUP   -1=initial state +1=final  +2=decayed
-c     ic(*,7) = Helicity
-c********************************************************************
-      implicit none
-
-      double precision pi
-      parameter (pi = 3.1415926d0)
-c
-c     Arguments
-c      
-      integer lun
-      integer nexternal, ic(7,*)
-      logical done
-      double precision P(0:4,*),wgt,aqcd,aqed,scale
-      integer ievent
-      character*140 buff2
-c
-c     Local
-c
-      integer i,j,k
-      character*(256) buff
-      double precision xdum1,xdum2
-c
-c     Global
-c
-      logical banner_open
-      integer lun_ban
-      common/to_banner/banner_open, lun_ban
-
-      data lun_ban/37/
-      data banner_open/.false./
-c-----
-c  Begin Code
-c-----     
-      buff2=' '
-      done=.false.
-      if (.not. banner_open) then
-         open (unit=lun_ban, status='scratch')
-         banner_open=.true.
-      endif
- 11   read(lun,'(a256)',end=99,err=99) buff
-      do while(index(buff,"<event") .eq. 0)
-         write(lun_ban,'(a)') buff
-         read(lun,'(a256)',end=99,err=99) buff
-      enddo
-      read(lun,*,err=11, end=11) nexternal,ievent,wgt,scale,aqed,aqcd
-      do i=1,nexternal
-         read(lun,*,err=99,end=99) ic(1,i),ic(6,i),(ic(j,i),j=2,5),
-     $     (p(j,i),j=1,3),p(0,i),p(4,i),xdum1,xdum2
-         ic(7,i)=xdum2
-      enddo
-      do while(index(buff,"</event") .eq. 0)
-         read(lun,'(a256)',end=99,err=99) buff
-         if(buff(1:1).eq.'#') buff2=buff(1:140)
-      enddo
-c      gal(1) = sqrt(4d0*pi*aqed)
-c      g      = sqrt(4d0*pi*aqcd)
-      return
- 99   done=.true.
-      return
- 55   format(i3,5e19.11)         
-      end
-
-      subroutine write_event(lun,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff)
-c********************************************************************
-c     Writes one event from data file #lun according to LesHouches
-c     ic(1,*) = Particle ID
-c     ic(2.*) = Mothup(1)
-c     ic(3,*) = Mothup(2)
-c     ic(4,*) = ICOLUP(1)
-c     ic(5,*) = ICOLUP(2)
-c     ic(6,*) = ISTUP   -1=initial state +1=final  +2=decayed
-c     ic(7,*) = Helicity
-c********************************************************************
-      implicit none
-c
-c     parameters
-c
-      double precision pi
-      parameter (pi = 3.1415926d0)
-c
-c     Arguments
-c      
-      integer lun, ievent
-      integer nexternal, ic(7,*)
-      double precision P(0:4,*),wgt
-      double precision aqcd, aqed, scale
-      character*140 buff
-c
-c     Local
-c
-      integer i,j,k
-c
-c     Global
-c
-
-c-----
-c  Begin Code
-c-----     
-c      aqed= gal(1)*gal(1)/4d0/pi
-c      aqcd = g*g/4d0/pi
-
-      write(lun,'(a)') '<event>'
-      write(lun,'(i2,i4,4e15.7)') nexternal,ievent,wgt,scale,aqed,aqcd
-      do i=1,nexternal
-         write(lun,51) ic(1,i),ic(6,i),(ic(j,i),j=2,5),
-     $     (p(j,i),j=1,3),p(0,i),p(4,i),0.,real(ic(7,i))
-      enddo
-      if(buff(1:1).eq.'#') write(lun,'(a)') buff(1:len_trim(buff))
-      write(lun,'(a)') '</event>'
-      return
- 51   format(i9,5i5,5e19.11,f3.0,f4.0)
-      end
-
-      subroutine write_comments(lun)
-c********************************************************************
-c     Outputs all of the banner comment lines back at the top of
-c     the file lun.
-c********************************************************************
-      implicit none
-c
-c     Arguments
-c
-      integer lun
-c
-c     Local
-c
-      character*(80) buff
-c
-c     Global
-c
-      logical banner_open
-      integer lun_ban
-      common/to_banner/banner_open, lun_ban
-
-c-----
-c  Begin Code
-c-----     
-c      write(*,*) 'Writing comments'
-      if (banner_open) then
-         rewind(lun_ban)
-         do while (.true.) 
-            read(lun_ban,'(a79)',end=99,err=99) buff
-            write(lun,'(a79)') buff
-c            write(*,*) buff
-         enddo
- 99      close(lun_ban)
-         banner_open = .false.
-      endif
-      end
-

=== removed file 'Template/NLO/Source/rw_events.short.f'
--- Template/NLO/Source/rw_events.short.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/rw_events.short.f	1970-01-01 00:00:00 +0000
@@ -1,162 +0,0 @@
-      subroutine read_event(lun,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,done)
-c********************************************************************
-c     Reads one event from data file #lun
-c     ic(*,1) = Particle ID
-c     ic(*,2) = Mothup(1)
-c     ic(*,3) = Mothup(2)
-c     ic(*,4) = ICOLUP(1)
-c     ic(*,5) = ICOLUP(2)
-c     ic(*,6) = ISTUP   -1=initial state +1=final  +2=decayed
-c     ic(*,7) = Helicity
-c********************************************************************
-      implicit none
-c
-c     parameters
-c
-      integer    MaxParticles
-      parameter (MaxParticles=15)
-      double precision pi
-      parameter (pi = 3.1415926d0)
-c
-c     Arguments
-c      
-      integer lun
-      integer nexternal, ic(7,MaxParticles)
-      logical done
-      double precision P(0:3,MaxParticles),wgt,aqcd,aqed,scale
-      integer ievent
-c
-c     Local
-c
-      integer i,j,k
-      character*(132) buff
-c
-c     Global
-c
-c      include 'coupl.inc'
-c      real*8          scale
-c      logical               fixed_ren_scale,fixed_fac_scale
-c      common/to_scale/scale,fixed_ren_scale,fixed_fac_scale
-
-      logical banner_open
-      integer lun_ban
-      common/to_banner/banner_open, lun_ban
-
-      data lun_ban/37/
-      data banner_open/.false./
-c-----
-c  Begin Code
-c-----     
-      done=.false.
-      if (.not. banner_open) then
-         open (unit=lun_ban, status='scratch')
-         banner_open=.true.
-      endif
- 11   read(lun,'(a132)',end=99,err=99) buff
-      do while(index(buff,"#") .ne. 0)
-         write(lun_ban,'(a)') buff
-         read(lun,'(a132)',end=99,err=99) buff
-      enddo
-      read(buff,*,err=11, end=11) nexternal,k,wgt,scale,aqed,aqcd
-      do j=1,7
-         read(lun,*,err=99,end=99) (ic(j,i),i=1,nexternal)!This is info
-      enddo      
-      do j=1,nexternal
-         read(lun,55,err=99,end=99) k,(p(i,j),i=0,3)
-      enddo
-c      gal(1) = sqrt(4d0*pi*aqed)
-c      g      = sqrt(4d0*pi*aqcd)
-      return
- 99   done=.true.
-      return
- 55   format(i3,4e19.11)         
-      end
-
-      subroutine write_event(lun,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed)
-c********************************************************************
-c     Writes one event from data file #lun according to LesHouches
-c     ic(*,1) = Particle ID
-c     ic(*,2) = Mothup(1)
-c     ic(*,3) = Mothup(2)
-c     ic(*,4) = ICOLUP(1)
-c     ic(*,5) = ICOLUP(2)
-c     ic(*,6) = ISTUP   -1=initial state +1=final  +2=decayed
-c     ic(*,7) = Helicity
-c********************************************************************
-      implicit none
-c
-c     parameters
-c
-      integer    MaxParticles
-      parameter (MaxParticles=15)
-      double precision pi
-      parameter (pi = 3.1415926d0)
-c
-c     Arguments
-c      
-      integer lun, ievent
-      integer nexternal, ic(7,MaxParticles)
-      double precision P(0:3,MaxParticles),wgt
-      double precision aqcd, aqed, scale
-c
-c     Local
-c
-      integer i,j,k
-c
-c     Global
-c
-
-c-----
-c  Begin Code
-c-----     
-c      aqed= gal(1)*gal(1)/4d0/pi
-c      aqcd = g*g/4d0/pi
-      write(lun,'(2i8,4e15.7)') nexternal,ievent,wgt,scale,aqed,aqcd
-      do j=1,7
-         write(lun,51) (ic(j,i),i=1,nexternal)  !This is info
-      enddo
-      do j=1,nexternal
-         write(lun,55) j,(p(i,j),i=0,3)
-      enddo
-      return
- 51   format(19i5)
- 55   format(i3,4e19.11)         
-      end
-
-      subroutine write_comments(lun)
-c********************************************************************
-c     Outputs all of the banner comment lines back at the top of
-c     the file lun.
-c********************************************************************
-      implicit none
-c
-c     Arguments
-c
-      integer lun
-c
-c     Local
-c
-      character*(80) buff
-c
-c     Global
-c
-      logical banner_open
-      integer lun_ban
-      common/to_banner/banner_open, lun_ban
-
-c-----
-c  Begin Code
-c-----     
-c      write(*,*) 'Writing comments'
-      if (banner_open) then
-         rewind(lun_ban)
-         do while (.true.) 
-            read(lun_ban,'(a79)',end=99,err=99) buff
-            write(lun,'(a79)') buff
-c            write(*,*) buff
-         enddo
- 99      close(lun_ban)
-         banner_open = .false.
-      endif
-      end
-

=== modified file 'Template/NLO/Source/rw_routines.f'
--- Template/NLO/Source/rw_routines.f	2015-01-06 11:33:59 +0000
+++ Template/NLO/Source/rw_routines.f	2017-05-24 22:17:46 +0000
@@ -1,391 +1,3 @@
-      subroutine load_para(npara,param,value)
-c----------------------------------------------------------------------
-c Read the params from the run_card.dat file
-c---------------------------------------------------------------------- 
-      implicit none
-c
-c     arguments
-c     
-      character*20 param(*),value(*)
-      integer npara
-c
-c     local
-c
-      logical fopened,done
-      integer iunit
-      character*20 ctemp
-      integer k,i,l1,l2,iproc
-      character*132 buff
-      data iunit/21/
-c
-c     global
-c
-      integer ngroup
-      common/to_group/ngroup
-c
-c----------
-c     start
-c----------
-c
-c     read the run_card.dat
-c
-      npara=0
-      param(1)=' '
-      value(1)=' '
-c
-c     open file
-c
-      call open_file(iunit,'run_card.dat',fopened)
-      if(.not.fopened) then
-         write(*,*) 'Error: File run_card.dat not found'
-         stop
-      else
-c
-c     first look for process-specific parameters
-c
-      done=.false.
-      do while(.not.done)  
-         read(iunit,'(a132)',end=20,err=20) buff
-         if(buff(1:1).ne.'#' .and. index(buff,"=").gt.0
-     $        .and. index(buff,"@").gt.0) then
-            l1=index(buff,"@")
-            l2=index(buff,"!")
-            if(l2.eq.0) l2=l1+20  !maybe there is no comment...
-            read(buff(l1+1:l2),*,err=11) iproc
-            if(iproc.ne.ngroup) cycle
-
-            l1=index(buff,"=")
-            l2=index(buff,"@")
-            if(l2-l1.lt.0) cycle
-            npara=npara+1
-c
-             value(npara)=buff(1:l1-1)
-             ctemp=value(npara)
-             call case_trap2(ctemp)
-             value(npara)=ctemp
-c
-             param(npara)=" "//buff(l1+1:l2-1)
-             ctemp=param(npara)
-             call case_trap2(ctemp)
-             param(npara)=ctemp
-c
- 11          cycle
-         endif
-      enddo
- 20   rewind(iunit)
-c
-c     read in values
-c
-      done=.false.
-      do while(.not.done)  
-         read(iunit,'(a132)',end=96,err=96) buff
-         if(buff(1:1).ne.'#' .and. index(buff,"=").gt.0
-     $        .and. index(buff,"@").le.0) then
-            l1=index(buff,"=")
-            l2=index(buff,"!")
-            if(l2.eq.0) l2=l1+20  !maybe there is no comment...
-            if(l2-l1.lt.0) cycle
-            npara=npara+1
-c
-             value(npara)=buff(1:l1-1)
-             ctemp=value(npara)
-             call case_trap2(ctemp)
-             value(npara)=ctemp
-c
-             param(npara)=" "//buff(l1+1:l2-1)
-             ctemp=param(npara)
-             call case_trap2(ctemp)
-             param(npara)=ctemp
-c
-         endif
-      enddo
- 96   close(iunit)
-      endif
-c
-c     open file
-c
-c
-c     tjs modified 11-16-07 to include grid_card.dat
-c
-      call open_file(iunit,'grid_card.dat',fopened)
-      if(fopened) then
-c
-c     first look for process-specific parameters
-c
-      done=.false.
-      do while(.not.done)  
-         read(iunit,'(a132)',end=30,err=30) buff
-         if(buff(1:1).ne.'#' .and. index(buff,"=").gt.0
-     $        .and. index(buff,"@").gt.0) then
-            l1=index(buff,"@")
-            l2=index(buff,"!")
-            if(l2.eq.0) l2=l1+20  !maybe there is no comment...
-            read(buff(l1+1:l2),*,err=21) iproc
-            if(iproc.ne.ngroup) cycle
-
-            l1=index(buff,"=")
-            l2=index(buff,"@")
-            if(l2-l1.lt.0) cycle
-            npara=npara+1
-c
-             value(npara)=buff(1:l1-1)
-             ctemp=value(npara)
-             call case_trap2(ctemp)
-             value(npara)=ctemp
-c
-             param(npara)=" "//buff(l1+1:l2-1)
-             ctemp=param(npara)
-             call case_trap2(ctemp)
-             param(npara)=ctemp
-c
- 21          cycle
-         endif
-       enddo
- 30   rewind(iunit)
-c
-c     read in values
-c
-      done=.false.
-      do while(.not.done)  
-         read(iunit,'(a132)',end=99,err=99) buff
-         if(buff(1:1).ne.'#' .and. index(buff,"=").gt.0
-     $        .and. index(buff,"@").le.0) then
-            l1=index(buff,"=")
-            l2=index(buff,"!")
-            if(l2.eq.0) l2=l1+20  !maybe there is no comment...
-            if(l2-l1.lt.0) cycle
-            npara=npara+1
-c
-             value(npara)=buff(1:l1-1)
-             ctemp=value(npara)
-             call case_trap2(ctemp)
-             value(npara)=ctemp
-c
-             param(npara)=" "//buff(l1+1:l2-1)
-c             write (*,*) param(npara),l1,l2
-             ctemp=param(npara)
-             call case_trap2(ctemp)
-             param(npara)=ctemp
-c             write(*,*) "New param:",param(npara)," = ", value(npara)
-c
-         endif
-      enddo
- 99   close(iunit)
-      endif
-
-      return
-      end
-
-
-
-      subroutine get_real(npara,param,value,name,var,def_value)
-c----------------------------------------------------------------------------------
-c     finds the parameter named "name" in param and associate to "value" in value 
-c----------------------------------------------------------------------------------
-      implicit none
-
-c
-c     arguments
-c
-      integer npara
-      character*20 param(*),value(*)
-      character*(*)  name
-      real*8 var,def_value
-      character*20 c_param,c_name
-c
-c     local
-c
-      logical found
-      integer i
-c
-c     start
-c
-      i=1
-      found=.false.
-      do while(.not.found.and.i.le.npara)
-          call firststring(c_param,param(i))
-          call firststring(c_name,name)
-         found = (c_param .eq. c_name)
-         if (found) read(value(i),*) var
-c         if (found) write (*,*) name,var
-         i=i+1
-      enddo
-      if (.not.found) then
-         write (*,*) "Warning: parameter ",name," not found"
-         write (*,*) "         setting it to default value ",def_value
-         var=def_value
-      endif
-      return
-
-      end
-c
-
-      subroutine get_integer(npara,param,value,name,var,def_value)
-c----------------------------------------------------------------------------------
-c     finds the parameter named "name" in param and associate to "value" in value 
-c----------------------------------------------------------------------------------
-      implicit none
-c
-c     arguments
-c
-      integer npara
-      character*20 param(*),value(*)
-      character*(*)  name
-      integer var,def_value
-      character*20 c_param,c_name
-c
-c     local
-c
-      logical found
-      integer i
-c
-c     start
-c
-      i=1
-      found=.false.
-      do while(.not.found.and.i.le.npara)
-         call firststring(c_param,param(i))
-          call firststring(c_name,name)
-         found = (c_param .eq. c_name)
-         if (found) read(value(i),*) var
-c         if (found) write (*,*) name,var
-         i=i+1
-      enddo
-      if (.not.found) then
-         write (*,*) "Warning: parameter ",name," not found"
-         write (*,*) "         setting it to default value ",def_value
-         var=def_value
-      endif
-      return
-
-      end
-c
-      subroutine get_int8(npara,param,value,name,var,def_value)
-c----------------------------------------------------------------------------------
-c     finds the parameter named "name" in param and associate to "value" in value 
-c----------------------------------------------------------------------------------
-      implicit none
-c
-c     arguments
-c
-      integer npara
-      character*20 param(*),value(*)
-      character*(*)  name
-      integer def_value
-      integer*8 var
-      character*20 c_param,c_name
-c
-c     local
-c
-      logical found
-      integer i
-c
-c     start
-c
-      i=1
-      found=.false.
-      do while(.not.found.and.i.le.npara)
-         call firststring(c_param,param(i))
-          call firststring(c_name,name)
-         found = (c_param .eq. c_name)
-         if (found) read(value(i),*) var
-c         if (found) write (*,*) name,var
-         i=i+1
-      enddo
-      if (.not.found) then
-         write (*,*) "Warning: parameter ",name," not found"
-         write (*,*) "         setting it to default value ",def_value
-         var=def_value
-      endif
-      return
-
-      end
-c
-      subroutine get_string(npara,param,value,name,var,def_value)
-c----------------------------------------------------------------------------------
-c     finds the parameter named "name" in param and associate to "value" in value 
-c----------------------------------------------------------------------------------
-      implicit none
-
-c
-c     arguments
-c
-      integer npara
-      character*20 param(*),value(*)
-      character*(*)  name
-      character*(*)  var,def_value
-      character*20 c_param,c_name
-c
-c     local
-c
-      logical found
-      integer i
-c
-c     start
-c
-      i=1
-      found=.false.
-      do while(.not.found.and.i.le.npara)
-         call firststring(c_param,param(i))
-          call firststring(c_name,name)
-         found = (c_param .eq. c_name)
-         if (found) read(value(i),*) var
-c         if (found) write (*,*) name,var
-         i=i+1
-      enddo
-      if (.not.found) then
-         write (*,*) "Warning: parameter ",name," not found"
-         write (*,*) "         setting it to default value ",def_value
-         var=def_value
-      endif
-      return
-
-      end
-c
-      subroutine get_logical(npara,param,value,name,var,def_value)
-c----------------------------------------------------------------------------------
-c     finds the parameter named "name" in param and associate to "value" in value 
-c----------------------------------------------------------------------------------
-      implicit none
-
-c
-c     arguments
-c
-      integer npara
-      character*20 param(*),value(*)
-      character*(*)  name
-      logical  var,def_value
-      character*20 c_param,c_name
-c
-c     local
-c
-      logical found
-      integer i
-c
-c     start
-c
-      i=1
-      found=.false.
-      do while(.not.found.and.i.le.npara)
-         call firststring(c_param,param(i))
-          call firststring(c_name,name)
-         found = (c_param .eq. c_name)
-         if (found) read(value(i),*) var
-c         if (found) write (*,*) name,var
-         i=i+1
-      enddo
-      if (.not.found) then
-         write (*,*) "Warning: parameter ",name," not found"
-         write (*,*) "         setting it to default value ",def_value
-         var=def_value
-      endif
-      return
-
-      end
-c
-
-
-
       subroutine case_trap2(name)
 c**********************************************************    
 c change the string to lowercase if the input is not
@@ -436,27 +48,3 @@
 
       return
       end
-
-c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-c ++
-c ++ firststring -> return the first "word" of string 
-c ++ & remove whitespaces around
-c ++ Needed to correct a bug in "get_" routines
-c ++ Michel Herquet - CP3 - 05-04-2006
-c ++
-c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-      subroutine firststring(first,string)
-
-      implicit none
-      character*(*) string
-      character*20 first
-      character*20 temp
-
-      temp=string
-      do while(temp(1:1) .eq. ' ') 
-	temp=temp(2:len(temp))
-      end do
-      first=temp(1:index(temp,' ')-1)
-
-      end

=== modified file 'Template/NLO/Source/setrun.f'
--- Template/NLO/Source/setrun.f	2016-09-08 23:15:34 +0000
+++ Template/NLO/Source/setrun.f	2017-05-24 22:17:46 +0000
@@ -20,7 +20,6 @@
 c
 c     include
 c
-      include 'genps.inc'
       include 'PDF/pdf.inc'
       include 'run.inc'
       include 'alfas.inc'
@@ -56,7 +55,6 @@
      &     idwtup,nprup,xsecup(maxpup),xerrup(maxpup),
      &     xmaxup(maxpup),lprup(maxpup)
 c
-      include 'nexternal.inc'
       logical gridrun,gridpack
       integer          iseed
       common /to_seed/ iseed
@@ -65,13 +63,8 @@
       common /event_normalisation/event_norm
       integer iappl
       common /for_applgrid/ iappl
+      integer idum
 C      
-      integer    maxflow
-      parameter (maxflow=999)
-      integer idup(nexternal,maxproc)
-      integer mothup(2,nexternal,maxproc)
-      integer icolup(2,nexternal,maxflow)
-      include 'born_leshouche.inc'
 c
 c----------
 c     start
@@ -192,7 +185,10 @@
         elseif(lpp(i).eq.-3) then
           idbmup(i)=-11
         elseif(lpp(i).eq.0) then
-          idbmup(i)=idup(i,1)
+           open (unit=71,status='old',file='initial_states_map.dat')
+           read (71,*,err=100)idum,idum,idbmup(1),idbmup(2)
+           close (71)
+c          idbmup(i)=idup(i,1)
         else
           idbmup(i)=lpp(i)
         endif
@@ -206,7 +202,7 @@
          call numberPDFm(1,nmemPDF(1))
          if (nmemPDF(1).eq.1) then
             nmemPDF(1)=0
-            lpdfvar(1)=0
+            lpdfvar(1)=.False.
          endif
       else
          nmemPDF(1)=0
@@ -215,6 +211,9 @@
       return
  99   write(*,*) 'error in reading'
       return
+ 100  write(*,*) '"initial_states_map.dat" not found (or incorrect'/
+     $     /' format) by "Source/setrun"'
+      stop 1
       end
 
 C-------------------------------------------------

=== removed file 'Template/NLO/Source/sum.f'
--- Template/NLO/Source/sum.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/sum.f	1970-01-01 00:00:00 +0000
@@ -1,55 +0,0 @@
-      program sum
-c********************************************************************************
-c     Program to combine results from all of the different sub amplitudes 
-c     and given total
-c     cross section and error.
-c*****************************************************************************
-      implicit none
-c
-c     Constants
-c
-      character*(*) rfile
-      parameter (rfile='result')
-      character*(*) symfile
-      parameter (symfile='symfact.dat')
-      integer   max_amps
-      parameter (max_amps=999)
-c
-c     local
-c
-      double precision xsec(max_amps), xerr(max_amps)
-      character*80 fname
-      integer i,j
-      double precision xtot,errtot
-c-----
-c  Begin Code
-c-----
-      xtot=0d0
-      errtot=0d0
-      open(unit=15,file=symfile,status='old',err=999)
-      do while (.true.)
-         read(15,*,err=99) i,j
-         if (j .gt. 0) then
-            if (i .lt. 10) then
-               write(fname,'(a,i1,a,a)') 'G',i,'/',rfile
-c               write(*,*) fname
-            else if (i .lt. 100) then
-               write(fname,'(a,i2,a,a)') 'G',i,'/',rfile
-c               write(*,*) fname
-            else if (i .lt. 1000) then
-               write(fname,'(a,i3,a,a)') 'G',i,'/',rfile
-c               write(*,*) fname
-            endif
-            open(unit=25,file=fname,status='old',err=95)
-            read(25,*) xsec(i), xerr(i)
-            xtot = xtot+xsec(i)*j
-            errtot=errtot+j*xerr(i)**2
-            write(*,'(2i4,2e12.4)') i,j, xsec(i),xerr(i)
- 95         close(25)
-         endif
-      enddo
- 99   write(*,*) 'done',xtot,sqrt(errtot)
-      close(15)
-      stop
- 999  write(*,*) 'error'
-      end

=== removed file 'Template/NLO/Source/write_banner.f'
--- Template/NLO/Source/write_banner.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/Source/write_banner.f	1970-01-01 00:00:00 +0000
@@ -1,116 +0,0 @@
-      program write_banner
-      implicit none
-c
-c     parameters
-c
-      integer    MaxParticles
-      parameter (MaxParticles=15)
-c
-c     Local
-c
-      integer lunr, lunw,luni
-      data luni,lunr,lunw/14,15,16/ ! unit numbers for reading and writing
-      integer ic(7,MaxParticles),next
-      double precision P(0:4,MaxParticles),wgt
-      real*8 sum,mxwgt
-      logical done
-      integer i,imax,j
-      integer nevent,nfound
-      character*25 infile,outfile
-      integer iseed
-      data iseed/9999/
-      character*30 process,QED,QCD
-      double precision scale,aqcd,aqed
-      integer ievent
-c--cuts
-c      double precision etmin(3:nexternal),etamax(3:nexternal)
-c      double precision                    r2min(3:nexternal,3:nexternal)
-c      double precision s_min(nexternal,nexternal)
-c      common/to_cuts/  etmin     ,etamax     , r2min, s_min
-c
-c     open the event file
-c
-      write(*,*) 'input the event file (e.g. Events/events.dat)'
-      read(*,'(a)')  infile
-      open(unit=lunr,file=infile,status='old')
-c
-c     open banner file
-c
-      write(*,*) 'input the banner file name (e.g. banner-events.dat)'
-      read(*,'(a)')  outfile
-      open(lunw,file=outfile,status='unknown')
-c
-c     gather the info
-c
-      call setpara('param_card.dat')
-c      call setcuts
-c      call get_seed(iseed)
-c
-c     All the info is gathered. Now start writing it out.
-c
-      call write_para(lunw)
-      write(lunw,'(a70)') '##                                                                    '
-      write(lunw,'(a70)') '##-------------------                                                 '
-      write(lunw,'(a70)') '## Run-time options                                                   '
-      write(lunw,'(a70)') '##-------------------                                                 '
-      write(lunw,'(a70)') '##                                                                    '
-c      write(lunw,'(a70)') '##********************************************************************'     
-c      write(lunw,'(a70)') '## Standard Cuts                                                     *'
-c      write(lunw,'(a70)') '##********************************************************************'    
-c      write(lunw,'(a13,8i8)')   '## Particle  ',(i,i=3,nexternal)
-c      write(lunw,'(a13,8f8.1)') '## Et       >',(etmin(i),i=3,nexternal)
-c      write(lunw,'(a13,8f8.1)') '## Eta      <',(etamax(i),i=3,nexternal)
-c      do j=3,nexternal-1
-c         write(lunw,'(a,i2,a,8f8.1)') '## d R #',j,'  >',(-0.0,i=3,j),
-c     &        (r2min(i,j),i=j+1,nexternal)
-c         do i=j+1,nexternal
-c            r2min(i,j)=r2min(i,j)**2 !Since r2 returns distance squared
-c         enddo
-c      enddo
-c      do j=3,nexternal-1
-c         write(lunw,'(a,i2,a,8f8.1)') '## s min #',j,'>',
-c     &        (s_min(i,j),i=3,nexternal)
-c      enddo
-      write(lunw,'(a70)') '##********************************************************************'    
-c
-c     Now write out specific information on the event set
-c
-      done=.false.
-      nevent=0
-      nfound=0
-      sum=0d0
-      mxwgt=-1d0
-      do while (.not. done)
-         call read_event(lunr,P,wgt,next,ic,ievent,scale,aqcd,aqed,done)
-         sum=sum+wgt
-         mxwgt = max(wgt,mxwgt)
-         nevent=nevent+1
-      enddo
-
-      write(lunw,'(a70)') '##                                                                    '
-      write(lunw,'(a70)') '##-------------------                                                 '
-      write(lunw,'(a70)') '## Event information                                                  '
-      write(lunw,'(a70)') '##-------------------                                                 '
-      write(lunw,'(a70)') '##                                                                    '
-      write(lunw,'(a70)') '##********************************************************************'    
-      write(lunw,'(a30,i10)')   '##  Number of Events       :  ',nevent
-      write(lunw,'(a30,e10.5)') '##  Integrated weight (pb) :  ',sum
-      write(lunw,'(a30,e10.5)') '##  Max wgt                :  ',mxwgt
-      write(lunw,'(a30,e10.5)') '##  Average wgt            :  ',sum/nevent
-      write(lunw,'(a70)') '##********************************************************************'    
-
-      rewind(lunr)
-      done = .false.
-      nevent=0
-      do while (.not. done)
-         call read_event(lunr,P,wgt,next,ic,ievent,scale,aqcd,aqed,done)
-         nevent=nevent+1
-         if (.not. done) then
-            call write_event(lunw,p,wgt,next,ic,nevent,scale,aqcd,aqed)
-         endif
-      enddo
-
-      close(lunw)
-      close(lunr)
-
-      end

=== modified file 'Template/NLO/SubProcesses/add_write_info.f'
--- Template/NLO/SubProcesses/add_write_info.f	2015-12-16 16:13:10 +0000
+++ Template/NLO/SubProcesses/add_write_info.f	2017-05-24 22:17:46 +0000
@@ -39,7 +39,7 @@
       external ran2
 
 c Jamp amplitudes of the Born (to be filled with a call the sborn())
-      double Precision amp2(maxamps), jamp2(0:maxamps)
+      double Precision amp2(ngraphs), jamp2(0:ncolor)
       common/to_amps/  amp2,       jamp2
 
 C iforest and other configuration info. Read once and saved.
@@ -70,8 +70,6 @@
       logical OnBW(-nexternal:0)
 
 c LesHouches info
-      integer maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
 c      include "leshouche.inc"
@@ -540,15 +538,15 @@
       do i=-ns,nexpart
          jpart(4,i)=icolalt(1,i)
          jpart(5,i)=icolalt(2,i)
-         if(i.eq.1.or.i.eq.2) then 
+         if(i.ge.1.and.i.le.nincoming) then 
             ito(i)=i            ! initial state particle
-         elseif(i.ge.3) then 
+         elseif(i.ge.nincoming+1) then 
             ito(i)=i+nres       ! final state particle
          elseif(i.le.-1.and.jpart(6,i).eq.2) then
             ires=ires+1
             ito(i)=2+ires       ! s-channel resonances
          else 
-            ito(i)=0
+            ito(i)=i
             if(i.eq.0) cycle
          endif
          if(jpart(2,i).lt.0.and.jpart(6,jpart(2,i)).ne.2) then
@@ -983,8 +981,6 @@
       implicit none
       include 'genps.inc'
       include 'nexternal.inc'
-      integer    maxflow
-      parameter (maxflow=999)
       integer i
       integer idup(nexternal,maxproc)
       integer mothup(2,nexternal,maxproc)

=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f	2016-09-08 23:15:34 +0000
+++ Template/NLO/SubProcesses/check_poles.f	2017-05-24 22:17:46 +0000
@@ -14,7 +14,7 @@
       integer return_code
       double precision tolerance, tolerance_default
       double precision, allocatable :: accuracies(:)
-      double precision accuracy
+      double precision accuracy2
       double precision ren_scale, energy
       include 'genps.inc'
       include 'nexternal.inc'
@@ -54,6 +54,7 @@
       logical first_time
       data first_time/.TRUE./
       include 'FKSParams.inc'
+      include 'mint.inc'
       
 C-----
 C  BEGIN CODE
@@ -100,6 +101,9 @@
           pmass_rambo(i-nincoming) = pmass(i)
       enddo
 
+      iconfig=1
+      ichan=1
+      iconfigs(1)=iconfig
 c Find the nFKSprocess for which we compute the Born-like contributions,
 c ie. which is a Born+g real-emission process
       do nFKSprocess=1,fks_configs
@@ -182,7 +186,7 @@
           call sborn(p_born, born)
           call sloopmatrix_thres(p_born,virt_wgts,tolerance,
      1 accuracies,return_code)
-          accuracy=accuracies(0)
+          accuracy2=accuracies(0)
 
           finite = virt_wgts(1,0)
           single = virt_wgts(2,0)
@@ -190,7 +194,7 @@
 
 C         If MadLoop was still in initialization mode, then skip this
 C         point for the checks
-          if (accuracy.lt.0.0d0) goto 200
+          if (accuracy2.lt.0.0d0) goto 200
 C         Otherwise, perform the check
           npointsChecked = npointsChecked +1
 

=== modified file 'Template/NLO/SubProcesses/cluster.f'
--- Template/NLO/SubProcesses/cluster.f	2015-12-16 16:13:10 +0000
+++ Template/NLO/SubProcesses/cluster.f	2017-05-24 22:17:46 +0000
@@ -213,7 +213,7 @@
                      if(prwidth(k,ignum).gt.ZERO) then
                      if(btest(mlevel,4))
      $                       write(*,*)'Adding resonance ',ignum,icmp(l)
-                        resmap(icmp(l),ignum)=.true.
+                        resmap(icmp(l),ignum,nFKSprocess)=.true.
                      endif
                   enddo
 c     proceed w/ next table, since there is no possibility,
@@ -278,8 +278,6 @@
       integer tprid(-max_branch:-1,n_max_cg)
       integer mapconfig(0:n_max_cg)
       common/c_configs_inc/iforest,sprop,tprid,mapconfig
-      integer maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
       common /c_leshouche_inc/idup,mothup,icolup,niprocs
@@ -341,8 +339,6 @@
       double precision prwidth(-max_branch:-1,n_max_cg)
       integer prow(-max_branch:-1,n_max_cg)
       common/c_props_inc/prmass,prwidth,prow
-      integer maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
       common /c_leshouche_inc/idup,mothup,icolup,niprocs
@@ -458,7 +454,8 @@
 c     check if we have constraint from onshell resonances
             foundbw=.true.
             do j=1,nbw
-               if(resmap(ibwlist(1,j),id_cl(nFKSprocess,idij,i)))then
+               if(resmap(ibwlist(1,j),id_cl(nFKSprocess,idij,i)
+     $              ,nFKSprocess))then
                   cycle
                endif
                foundbw=.false.
@@ -520,15 +517,18 @@
       include 'message.inc'
       include 'real_from_born_configs.inc'
       include 'run.inc'
+      include 'fks_info.inc'
 c Argument
       double precision p(0:3,nexternal)
 c Local
       integer i,j,k,n,idi,idj,idij,icgs(0:n_max_cg),nleft,iwin,jwin
-     $     ,iwinp,imap(nexternal,2)
+     $     ,iwinp,imap(nexternal,2),i_fks,j_fks
       double precision pcmsp(0:3),p1(0:3),pi(0:3),nr(0:3),nn2,ct,st
      $     ,minpt2ij,pt2ij(n_max_cl),zij(n_max_cl)
       double precision pz(0:3)
       data (pz(i),i=0,3)/1d0,0d0,0d0,1d0/
+      integer igraphscounter
+      data igraphscounter/149/
 c Common
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
@@ -540,6 +540,9 @@
       double precision dj, pydj, djb, pyjb, dot, SumDot, zclus
       external findmt,dj, pydj, djb, pyjb, dot, SumDot, zclus, combid
 
+      i_fks=fks_i_d(nFKSprocess)
+      j_fks=fks_j_d(nFKSprocess)
+      
       if (btest(mlevel,1))
      $   write (*,*)'New event'
 
@@ -626,7 +629,10 @@
                      endif
                   endif
 c     Check if smallest pt2 ("winner")
-                  if (pt2ij(idij).lt.minpt2ij) then
+                  if (pt2ij(idij).lt.minpt2ij .or.
+     &                 (pt2ij(idij).eq.minpt2ij .and.
+     &                 (i.eq.i_fks .or. i.eq.j_fks) .and.
+     &                 (j.eq.i_fks .or. j.eq.j_fks))  ) then
                      iwin=j
                      jwin=i
                      minpt2ij=pt2ij(idij)
@@ -799,6 +805,11 @@
                   exit
                endif
             enddo
+            if (igraphs(0).ge.2) then ! "randomize" the graph taken
+               igraphscounter=igraphscounter+1
+               igraphs(0)=1
+               igraphs(1)=igraphs(mod(igraphscounter,igraphs(0))+1)
+            endif
             zcl(n+1)=1
             cluster=.true.
             clustered=.true.
@@ -866,7 +877,10 @@
                      write(*,*)'     cf. djb: ',djb(pcl(0,idi))
                   endif
                endif
-               if (pt2ij(idij).lt.minpt2ij) then
+               if (pt2ij(idij).lt.minpt2ij .or.
+     &                 (pt2ij(idij).eq.minpt2ij .and.
+     &                 (i.eq.i_fks .or. i.eq.j_fks) .and.
+     &              (j.eq.i_fks .or. j.eq.j_fks))  ) then
                   iwin=j
                   jwin=i
                   minpt2ij=pt2ij(idij)

=== modified file 'Template/NLO/SubProcesses/cluster.inc'
--- Template/NLO/SubProcesses/cluster.inc	2013-03-22 14:03:58 +0000
+++ Template/NLO/SubProcesses/cluster.inc	2017-05-24 22:17:46 +0000
@@ -9,7 +9,7 @@
       logical heavyrad(n_max_cg)
       common/cl_map/id_cl,heavyrad
 c     resmap gives potential resonances for diagrams
-      logical resmap(n_max_cl,n_max_cg)
+      logical resmap(n_max_cl,n_max_cg,fks_configs)
       common/res_map/resmap
       double precision pt2ijcl(nexternal),zcl(nexternal),mt2ij(nexternal),mt2last
       double precision pcl(0:4,n_max_cl) ! 4 is mass**2

=== modified file 'Template/NLO/SubProcesses/cuts.f'
--- Template/NLO/SubProcesses/cuts.f	2015-12-16 16:13:10 +0000
+++ Template/NLO/SubProcesses/cuts.f	2017-05-24 22:17:46 +0000
@@ -428,8 +428,6 @@
       double precision pmass(nexternal)
       common/to_mass/pmass
 c PDG codes of particles
-      integer maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
       common /c_leshouche_inc/idup,mothup,icolup,niprocs
@@ -868,8 +866,6 @@
       implicit none
       include "genps.inc"
       include 'nexternal.inc'
-      integer maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
 c      include 'leshouche.inc'
@@ -887,8 +883,6 @@
       implicit none
       include "genps.inc"
       include 'nexternal.inc'
-      integer    maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc)
       integer mothup(2,nexternal,maxproc)
       integer icolup(2,nexternal,maxflow)

=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2017-03-08 10:26:43 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2017-05-24 22:17:46 +0000
@@ -22,20 +22,10 @@
 c
 c     Global
 c
-      integer                                      nsteps
-      character*40          result_file,where_file
-      common /sample_status/result_file,where_file,nsteps
-      integer ngroup
-      common/to_group/ngroup
-      data ngroup/0/
 cc
       include 'run.inc'
       include 'coupl.inc'
       
-      double precision twgt, maxwgt,swgt(maxevents)
-      integer                             lun, nw
-      common/to_unwgt/twgt, maxwgt, swgt, lun, nw
-
 c Vegas stuff
       common/tosigint/ndim
 
@@ -136,18 +126,6 @@
         n1(i)=0
       enddo
       
-      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
       call setpara('param_card.dat')   !Sets up couplings and masses
       call setcuts               !Sets up cuts and particle masses
@@ -430,7 +408,7 @@
       nbody=.true.
       calculatedBorn=.false.
       call get_born_nFKSprocess(nFKS_picked,nFKS_born)
-      call update_fks_dir(nFKS_born,iconfig)
+      call update_fks_dir(nFKS_born)
       jac=1d0
       call generate_momenta(ndim,iconfig,jac,x,p)
       if (p_born(0,1).lt.0d0) goto 12
@@ -465,7 +443,7 @@
          wgt_me_born=0d0
          wgt_me_real=0d0
          jac=MC_int_wgt
-         call update_fks_dir(iFKS,iconfig)
+         call update_fks_dir(iFKS)
          call generate_momenta(ndim,iconfig,jac,x,p)
          if (p_born(0,1).lt.0d0) cycle
          call compute_prefactors_n1body(vegas_wgt,jac)
@@ -525,16 +503,16 @@
       return
       end
 
-      subroutine update_fks_dir(nFKS,iconfig)
+      subroutine update_fks_dir(nFKS)
       implicit none
-      integer nFKS,iconfig
+      integer nFKS
       integer              nFKSprocess
       common/c_nFKSprocess/nFKSprocess
       nFKSprocess=nFKS
       call fks_inc_chooser()
       call leshouche_inc_chooser()
       call setcuts
-      call setfksfactor(iconfig,.false.)
+      call setfksfactor(.false.)
       return
       end
       

=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2017-01-30 15:18:32 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2017-05-24 22:17:46 +0000
@@ -383,7 +383,7 @@
             endif
 c Randomly pick the contribution that will be written in the event file
             call pick_unweight_contr(iFKS_picked)
-            call update_fks_dir(iFKS_picked,iconfig)
+            call update_fks_dir(iFKS_picked)
             call fill_rwgt_lines
             call finalize_event(x,weight,lunlhe,putonshell)
          enddo
@@ -724,16 +724,12 @@
          write(*,*) 'Not subdividing B.W.'
          lbw(0)=0
       else
-         lbw(0)=1
-         jconfig=dconfig*1000.1
-         write(*,*) 'Using dconfig=',jconfig
-         call DeCode(jconfig,lbw(1),3,nexternal)
-         write(*,*) 'BW Setting ', (lbw(j),j=1,nexternal-2)
+         write(*,*) 'Error BW setting: not supported at NLO'
+         stop 1
       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
 
@@ -830,7 +826,7 @@
             call get_born_nFKSprocess(nFKS_in,nFKS_out)
             nFKS_picked_nbody=nFKS_out
          endif
-         call update_fks_dir(nFKS_picked_nbody,iconfig)
+         call update_fks_dir(nFKS_picked_nbody)
          jac=1d0
          call generate_momenta(ndim,iconfig,jac,x,p)
          if (p_born(0,1).lt.0d0) goto 12
@@ -865,7 +861,7 @@
             wgt_me_real=0d0
             wgt_me_born=0d0
             iFKS=proc_map(proc_map(0,1),i)
-            call update_fks_dir(iFKS,iconfig)
+            call update_fks_dir(iFKS)
             jac=1d0/vol1
             probne=1d0
             gfactsf=1.d0
@@ -1175,8 +1171,6 @@
       common/event_attributes/nattr,npNLO,npLO
       integer              nFKSprocess
       common/c_nFKSprocess/nFKSprocess
-      integer    maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
       common /c_leshouche_inc/idup,mothup,icolup,niprocs
@@ -1211,16 +1205,16 @@
       end
 
 
-      subroutine update_fks_dir(nFKS,iconfig)
+      subroutine update_fks_dir(nFKS)
       implicit none
-      integer nFKS,iconfig
+      integer nFKS
       integer              nFKSprocess
       common/c_nFKSprocess/nFKSprocess
       nFKSprocess=nFKS
       call fks_inc_chooser()
       call leshouche_inc_chooser()
       call setcuts
-      call setfksfactor(iconfig,.true.)
+      call setfksfactor(.true.)
       return
       end
 

=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
--- Template/NLO/SubProcesses/fks_singular.f	2017-01-25 15:41:36 +0000
+++ Template/NLO/SubProcesses/fks_singular.f	2017-05-24 22:17:46 +0000
@@ -649,7 +649,7 @@
      &                         fkssymmetryfactorDeg,ngluons,nquarks
       integer            mapconfig(0:lmaxconfigs), iconfig
       common/to_mconfigs/mapconfig,                iconfig
-      Double Precision amp2(maxamps), jamp2(0:maxamps)
+      Double Precision amp2(ngraphs), jamp2(0:ncolor)
       common/to_amps/  amp2,          jamp2
       double precision   diagramsymmetryfactor
       common /dsymfactor/diagramsymmetryfactor
@@ -797,7 +797,7 @@
      &                         fkssymmetryfactorDeg,ngluons,nquarks
       integer            mapconfig(0:lmaxconfigs), iconfig
       common/to_mconfigs/mapconfig,                iconfig
-      Double Precision amp2(maxamps), jamp2(0:maxamps)
+      Double Precision amp2(ngraphs), jamp2(0:ncolor)
       common/to_amps/  amp2,          jamp2
       double precision   diagramsymmetryfactor
       common /dsymfactor/diagramsymmetryfactor
@@ -1612,8 +1612,6 @@
       include 'fks_info.inc'
       include 'genps.inc'
       integer k,ict,iFKS
-      integer    maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      $     icolup(2,nexternal,maxflow),niprocs
       common /c_leshouche_inc/idup,mothup,icolup,niprocs

=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
--- Template/NLO/SubProcesses/genps_fks.f	2016-09-30 08:53:12 +0000
+++ Template/NLO/SubProcesses/genps_fks.f	2017-05-24 22:17:46 +0000
@@ -1838,6 +1838,7 @@
      $        ,dum3,dum3,jac,s)
       else
          write (*,*) 'ERROR #39 in genps_fks.f',s_mass,smin
+         jac=-1d0
       endif
       tau=s/stot
       jac=jac/stot

=== modified file 'Template/NLO/SubProcesses/iproc_map.f'
--- Template/NLO/SubProcesses/iproc_map.f	2015-12-16 16:13:10 +0000
+++ Template/NLO/SubProcesses/iproc_map.f	2017-05-24 22:17:46 +0000
@@ -12,8 +12,6 @@
       COMMON /SUBPROC/ PD, IPROC
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
-      integer maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
       common /c_leshouche_inc/idup,mothup,icolup,niprocs
@@ -224,8 +222,6 @@
       common/c_flavour_map/flavour_map
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
-      integer maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
       common /c_leshouche_inc/idup,mothup,icolup,niprocs

=== modified file 'Template/NLO/SubProcesses/leshouche_inc_chooser.f'
--- Template/NLO/SubProcesses/leshouche_inc_chooser.f	2015-12-16 16:13:10 +0000
+++ Template/NLO/SubProcesses/leshouche_inc_chooser.f	2017-05-24 22:17:46 +0000
@@ -7,8 +7,6 @@
       integer i,j,k
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
-      integer maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
       common /c_leshouche_inc/idup,mothup,icolup,niprocs

=== modified file 'Template/NLO/SubProcesses/makefile_fks_dir'
--- Template/NLO/SubProcesses/makefile_fks_dir	2016-06-28 14:44:52 +0000
+++ Template/NLO/SubProcesses/makefile_fks_dir	2017-05-24 22:17:46 +0000
@@ -37,25 +37,18 @@
 	 MC_integer.o $(reweight_xsec_events_pdf_dummy)			\
 	 $(applgrid_interface)
 
-# Files needed for vegas, mint & mintMC
+# Files needed for mintFO & mintMC
 RUN= $(FO_ANALYSE) $(FILES) cuts.o pythia_unlops.o recluster.o 		\
 	fill_MC_mshell.o born_hel.o open_output_files.o			\
 	add_write_info.o BinothLHA.o madfks_plot.o
 
-# Files for vegas
-VEGAS= $(RUN) driver_vegas.o vegas2.o handling_lhe_events.o
-
-# Files for vegas
+# Files for mintFO
 MINTFO= $(RUN) driver_mintFO.o mint-integrator2.o handling_lhe_events.o
 
 # Files for mintMC
 MINTMC= $(RUN) mint-integrator2.o driver_mintMC.o write_event.o		\
 	handling_lhe_events.o
 
-# Files for Born reweighting
-BORN_REWEIGHT= $(RUN) driver_reweight.o handling_lhe_events.o		\
-	write_event.o
-
 # Files for check_poles
 POLES= $(FILES) cuts.o pythia_unlops.o recluster.o check_poles.o	\
        BinothLHA.o born_hel.o
@@ -96,10 +89,6 @@
 	rm fks_singular.o
 	strip gensym
 
-test_Sij: symmetry_fks_test_Sij.o $(TEST)
-	$(FC) $(LDFLAGS) -o test_Sij symmetry_fks_test_Sij.o $(TEST) $(APPLLIBS) $(LINKLIBS) $(FJLIBS)
-	strip test_Sij
-
 test_ME: symmetry_fks_test_ME.o $(TEST)
 	$(FC) $(LDFLAGS) -o test_ME symmetry_fks_test_ME.o $(TEST) $(APPLLIBS) $(LINKLIBS) $(FJLIBS)
 	rm symmetry_fks_test_ME.o
@@ -114,11 +103,6 @@
 	$(FC) -o check_poles $(POLES) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(LDFLAGS)
 	strip check_poles
 
-madevent_vegas: $(VEGAS) $(libmadloop) makefile $(LIBS)
-	$(FC) -o madevent_vegas $(VEGAS) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS) $(LDFLAGS)
-	rm handling_lhe_events.o
-	strip madevent_vegas
-
 madevent_mintMC: $(MINTMC) $(libmadloop) makefile $(LIBS)
 	$(FC) -o madevent_mintMC $(MINTMC) $(NLOLIBS) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS) $(LDFLAGS)
 	rm handling_lhe_events.o
@@ -129,10 +113,6 @@
 	rm handling_lhe_events.o
 	strip madevent_mintFO
 
-madevent_reweight: $(BORN_REWEIGHT) $(libmadloop) makefile $(LIBS)
-	$(FC) -o madevent_reweight $(BORN_REWEIGHT) $(APPLLIBS) $(LINKLIBS) $(FJLIBS) $(FO_EXTRAPATHS) $(FO_EXTRALIBS) $(LDFLAGS)
-	rm handling_lhe_events.o
-
 reweight_xsec_events: $(RWGFILES) makefile $(LIBS)
 	$(FC) -o reweight_xsec_events $(RWGFILES) $(LINKLIBS) $(FJLIBS) $(LDFLAGS)
 	rm handling_lhe_events.o

=== modified file 'Template/NLO/SubProcesses/mint-integrator2.f'
--- Template/NLO/SubProcesses/mint-integrator2.f	2017-03-03 23:30:52 +0000
+++ Template/NLO/SubProcesses/mint-integrator2.f	2017-05-24 22:17:46 +0000
@@ -163,7 +163,9 @@
             do kchan=1,nchans
                do kint=0,nint_used
                   xgrid(kint,kdim,kchan)=dble(kint)/nint_used
-                  nhits(kint,kdim,kchan)=0
+                  if(kint.gt.0) then
+                     nhits(kint,kdim,kchan)=0
+                  endif
                enddo
                regridded(kchan)=.true.
             enddo

=== modified file 'Template/NLO/SubProcesses/montecarlocounter.f'
--- Template/NLO/SubProcesses/montecarlocounter.f	2017-03-03 23:36:24 +0000
+++ Template/NLO/SubProcesses/montecarlocounter.f	2017-05-24 22:17:46 +0000
@@ -984,7 +984,7 @@
       double precision p_born(0:3,nexternal-1)
       common/pborn/p_born
 
-      double Precision amp2(maxamps), jamp2(0:maxamps)
+      double Precision amp2(ngraphs), jamp2(0:ncolor)
       common/to_amps/  amp2,       jamp2
 
       integer i_fks,j_fks

=== modified file 'Template/NLO/SubProcesses/pythia_unlops.f'
--- Template/NLO/SubProcesses/pythia_unlops.f	2015-12-16 16:13:10 +0000
+++ Template/NLO/SubProcesses/pythia_unlops.f	2017-05-24 22:17:46 +0000
@@ -18,8 +18,6 @@
       double precision p(0:3,nexternal),eCM
       logical passUNLOPScuts
       INTEGER I, J
-      integer maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
       common /c_leshouche_inc/idup,mothup,icolup,niprocs

=== modified file 'Template/NLO/SubProcesses/reweight.f'
--- Template/NLO/SubProcesses/reweight.f	2015-12-16 16:13:10 +0000
+++ Template/NLO/SubProcesses/reweight.f	2017-05-24 22:17:46 +0000
@@ -73,7 +73,7 @@
          if (iimode.eq.2) then
             gamma=gamma+CF*alphasq0**2*log(Q1/q0)*kappa/(2d0*pi**2) ! A2
          endif
-         if (iipdg.gt.NF) then ! include mass effects
+         if (abs(iipdg).gt.NF) then ! include mass effects
             qom=q0/mass(iipdg)
             gamma=gamma+CF*alphasq0/pi/2d0*( 0.5d0 - qom*atan(1d0/qom) -
      $           (1d0-0.5d0*qom**2)*log(1d0+1d0/qom**2) )
@@ -490,10 +490,10 @@
       if((ipart(1,ida1).ge.1.and.ipart(1,ida1).le.2).or.
      $   (ipart(1,ida2).ge.1.and.ipart(1,ida2).le.2))then
 c     Check if ida1 is outgoing parton or ida2 is outgoing parton
-         if(.not.islast.and.ipart(1,ida2).ge.1
-     $        .and.ipart(1,ida2).le.2.and.isparton(idda1).or.ipart(1
-     $        ,ida1).ge.1.and.ipart(1
-     $        ,ida1).le.2.and.isparton(idda2))then
+         if((.not.islast) .and. ((ipart(1,ida2).ge.1 .and. ipart(1
+     $        ,ida2).le.2 .and. isparton(idda1)) .or. (ipart(1
+     $        ,ida1).ge.1 .and. ipart(1,ida1).le.2 .and.
+     $        isparton(idda2))))then
            ispartonvx=.true.
         else
            ispartonvx=.false.
@@ -1222,8 +1222,6 @@
      $     ,fks_configs,lmaxconfigs),jlast(2)
       real*8 q2bck(2)
       common /to_rw/jlast,njetstore,iqjetstore,q2bck
-      integer maxflow
-      parameter (maxflow=999)
       integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
      &     icolup(2,nexternal,maxflow),niprocs
       common /c_leshouche_inc/idup,mothup,icolup,niprocs

=== removed file 'Template/NLO/SubProcesses/reweight_events.f'
--- Template/NLO/SubProcesses/reweight_events.f	2012-09-18 09:48:33 +0000
+++ Template/NLO/SubProcesses/reweight_events.f	1970-01-01 00:00:00 +0000
@@ -1,375 +0,0 @@
-C***************************************************************
-c Program that reads an Les Houches event file (with unweighted
-c Born events) and calculates the virtual corrections for each
-c event. The weights in the event file should sum to the total
-c cross section.
-c
-c It must be run from a P* directory, because nexternal.inc,
-c coupl.inc, fks.inc, pmass.inc  and run.inc should all be there.
-c
-c Call to virtual is via the BinothLHA() interface subroutine
-c (in the corresponding file.)
-c
-C***************************************************************
-c The virt_wgt coming from BinothLHA() should have the same
-c adjusted identical particle symmetry factor as the Born in
-c the MadFKS framework (i.e. divided by dble(ngluons)).
-c This is also needed for normal inclusion of the virtual, so
-c in general, nothing needs to be changed in BinothLHA().
-C***************************************************************
-c
-      program reweight_events
-      implicit none
-      include "nexternal.inc"
-      include "coupl.inc"
-      include "run.inc"
-      integer MaxParticles
-      parameter (MaxParticles=20)
-      integer lun,lun2,lun3
-      parameter (lun=99,lun2=98,lun3=97)
-      integer npart, ic(7,MaxParticles),i,j,k,l
-      logical done,firsttime(2)
-      data firsttime/.true.,.true./
-      double precision P(0:4,MaxParticles),wgt,aqcd,aqed,xscale,virt_wgt
-      integer ievent,nevent,unwgt,mxevt,unw_eventp,unw_eventn
-      character*140 buff2,filename_in, filename_out,filename_unw
-      double precision sumw,sumw_abs,mxwgt,p_born(0:3,nexternal-1)
-      double precision born_wgt,pdf_wgt
-      external dlum
-      integer ID_save(MaxParticles)
-      double complex wgt1(2)
-      integer iseed
-      data iseed/1/
-      double precision rnd,ran2
-      external ran2
-      logical calculatedBorn
-      common/ccalculatedBorn/calculatedBorn
-      character*1 proceed
-
-      write (*,*) 'Give Les Houches event file with unweighted'//
-     &     ' Born events'     
-      read (*,'(a)') filename_in
-      write (*,*) 'Unweight the reweighted events? (0: no, 1: yes)'
-      read (*,*) unwgt
-      if (unwgt.eq.0) then
-         write (*,*) 'No unweighting'
-      elseif (unwgt.eq.1) then
-         write (*,*) 'Also unweighting the events'
-      else
-         write (*,*) 'Error: did not understand input',unwgt
-      endif
-
-      filename_out="rwgt_events.lhe"
-      filename_unw="unwgt_events.lhe"
-
-      open (unit=lun,file='nbodyonly.fks',status='old',err=96)
-      read (lun,'(a)') proceed
-      close(lun)
-      if (proceed.eq.'N') then
-         write (*,*) 'Can only reweight in a "nbodyonly" directory'
-         stop
-      endif
-
-      call setpara('param_card.dat')   !Sets up couplings and masses
-      call setrun
-
-      open (unit=lun,file=filename_in,status='old',err=99)
-      
-      done = .false.
-      sumw=0d0
-      mxwgt=0d0
-      nevent=0
-      write (*,*) ''
-      write (*,*) 'Parsing events to do some simple tests...'
-      do while (.not.done)
-         call read_event(lun,P,wgt,npart,ic,ievent,
-     &        xscale,aqcd,aqed,buff2,done)
-         if (done) cycle
-         nevent=nevent+1
-         if(mod(nevent,5000).eq.0) print *,'At event ',nevent
-         if(npart.gt.20) then
-            write (*,*) 'Error: Too many particles in event: ',npart
-            stop
-         endif
-         sumw=sumw+wgt
-         mxwgt = max(wgt,mxwgt)
-         if (mxwgt.ne.wgt) then
-            write (*,*) 'Error: not all events have equal weights'
-            stop
-         endif
-         j=0
-         do i=1,npart
-            if (abs(ic(6,i)).eq.1) then
-               j=j+1
-               if (nevent.eq.1) then
-                  id_save(j)=ic(1,i)
-               else
-                  if (ic(1,i).ne.id_save(j) .and. firsttime(1)) then
-                     firsttime(1)=.false.
-                     write (*,*) '* Warning: not all events'// 
-     &                    'have the same final state partons'
-                  endif
-               endif
-            endif
-            if (ic(7,i).ne.0 .and. firsttime(2)) then
-               firsttime(2)=.false.
-               write (*,*) '* Warning: helicity info in Born events'//
-     &              ' will be ignored'
-            endif
-         enddo
-      enddo
-
-      write (*,*) '...',nevent,' events parsed'
-      write (*,*) ''
-      write (*,*) 'Reweighting events...'
-
-      rewind(lun)
-
-      open (unit=lun2,file=filename_out,status='unknown',err=98)
-      
-      sumw=0d0
-      sumw_abs=0d0
-      mxwgt=0d0
-      call fill_common_blocks()
-      
-      do i=1,nevent
-         do j=1,140
-            buff2(j:j)=' '
-         enddo
-         if (nevent.gt.50000) then
-            if(mod(i,10000).eq.0) print *,'At event ',i
-         elseif (nevent.gt.5000) then
-            if(mod(i,1000).eq.0) print *,'At event ',i
-         else
-            if(mod(i,100).eq.0) print *,'At event ',i
-         endif
-c Read event
-         call read_event(lun,P,wgt,npart,ic,ievent,
-     &        xscale,aqcd,aqed,buff2,done)
-         if (done) then
-            write (*,*) 'ERROR, some events are missing '//
-     &           'from event file',i,nevent
-            stop
-         endif
-         j=0
-         do l=1,npart
-            if (abs(ic(6,l)).eq.1) then
-               j=j+1
-               do k=0,3
-                  p_born(k,j)=p(k,l)
-               enddo
-            endif
-         enddo
-         calculatedBorn=.false.
-c Set scale and couplings (on an event-by-event basis)
-         call set_alphaS(p_born)
-         if (abs((aqcd-g**2/(16d0*atan(1d0)))/aqcd).gt.1d-5) then
-            write (*,*) "Error: inconsistency in strong coupling"
-            write (*,*) aqcd,g**2/(16d0*atan(1d0)),scale
-            stop
-         endif
-         call get_helicity_simplified()
-c Compute Born
-         call sborn(p_born,wgt1)
-         born_wgt=dble(wgt1(1))
-c Compute Virtual
-         call BinothLHA(p_born,born_wgt,virt_wgt)
-c Reweight the events
-         wgt=wgt * virt_wgt/born_wgt
-c Write event
-         call write_event(lun2,P,wgt,npart,ic,ievent,
-     &        xscale,aqcd,aqed,buff2)
-         sumw=sumw+wgt
-         sumw_abs=sumw_abs+abs(wgt)
-         if (abs(wgt).gt.abs(mxwgt)) then
-            mxwgt = wgt
-            mxevt=i
-         endif
-      enddo
-
-
-      write (*,*)
-     &     'Reweighted events are written to file ',filename_out
-      write (*,*) 'Cross section:',sumw
-      write (*,*) 'average abs. of weights:',sumw_abs/dble(nevent)
-      write (*,*) 'max wgt:',mxwgt,' for event',mxevt
-      write (lun2,'(a)')'</LesHouchesEvents>'
-c Do the unweighting
-      if (unwgt.eq.1) then
-         write (*,*) ''
-         write (*,*) 'Unweighting...'
-         rewind(lun2)
-         open(unit=lun3,file=filename_unw,status='unknown',err=97)
-
-         unw_eventp=0
-         unw_eventn=0
-         do i=1,nevent
-            call read_event(lun2,P,wgt,npart,ic,ievent,
-     &           xscale,aqcd,aqed,buff2,done)
-            if (done) then
-               write (*,*) 'ERROR, some events are missing '//
-     &              'from event file',i,nevent
-               stop
-            endif
-            rnd=ran2()
-            if (rnd.lt.abs(wgt/mxwgt)) then
-               if (wgt.gt.0) then
-                  wgt=1d0
-                  unw_eventp=unw_eventp+1
-               else
-                  wgt=-1d0
-                  unw_eventn=unw_eventn+1
-               endif
-               call write_event(lun3,P,wgt,npart,ic,ievent,
-     &              xscale,aqcd,aqed,buff2)
-            endif
-         enddo
-         
-         write (*,*)
-     &        'Unweighted events are written to file ',filename_unw
-         write (*,*) 'Number of events written is',unw_eventp+unw_eventn
-         write (*,*) 'positive:',unw_eventp
-         write (*,*) 'negative:',unw_eventn
-         write (*,*) 'Unweighting efficiency (in %):',
-     &        dble(unw_eventp+unw_eventn)/dble(nevent)*100d0
-         write(lun3,'(a)')'</LesHouchesEvents>'
-         close(lun3)
-      endif
-
-
-      close(lun)
-      close(lun2)
-      return