Merge lp:~madnklo/mg5amcnlo/MadNkLO_ISR into lp:~madnklo/mg5amcnlo/MadNkLO

[Valentin Hirschi]

At long last, here is the alpha version of the ISR.
There is still quite some checks to be done and the option "group_subprocesses" still need testing.
That being said however, all is in place and one can already investigate NLO processes.
NNLO processes and beyond should already be structurally OK, but that needs further testing too.

In any case, for now the following can be tested/run:

MG5_aMC>set subtraction_currents_scheme colorful
MG5_aMC>import model loop_sm
MG5_aMC>set group_subprocesses False
MG5_aMC>define p = u d
MG5_aMC>define pbar = u~ d~
MG5_aMC>generate p pbar > z QED=1 --NLO --beam_types=proton@(1,2,3,4,-1,-2,-3,-4,21)
MG5_aMC>display contributions
MG5_aMC>output TEST_ISR_ddx_z
MG5_aMC>display contributions --format=4
MG5_aMC>display contributions --format=2
MG5_aMC>launch TEST_ISR_ddx_z -i
MG5_aMC>display integrands --format=4
ME7@NLO::TEST_ISR_ddx_z > test_IR_poles --integrands=BF1 --counterterms=all --set_PDFs_to_unity=False
ME7@NLO::TEST_ISR_ddx_z > test_IR_limits --limits=F(1) --integrands=BF1 --counterterms=all --set_PDFs_to_unity=False

You can also already integrate the B, V, BF1 and BF2 contributions, while for the R one we're still missing appropriate mappings for 2>1.

I'll update this posts with progress on the remaining issues, but ISR is structurally done now!

[Valentin Hirschi (valentin-hirschi)]

I should add that the process definition is overly complicated just so as to illustrate the new option '--beam_types'.
But it's perfectly possible to simply type:

MG5_aMC>generate p p > z QED=1 --NLO

and things will be set automatically.

[lp:~madnklo/mg5amcnlo/MadNkLO_ISR]

531. By Valentin Hirschi on 2018-08-20

1. Fix to integrands
2. Skeletton for p p > X soft mapping and walker

532. By Nicolas Deutschmann on 2018-08-27

Implemented the Soft vs Initial state mapping to higher/lower multiplicity. All methods raise a NotImplemented error for now as they are *not* tested yet.

533. By Nicolas Deutschmann on 2018-08-27

Soft Vs Initial Mapping corrected to pass invertibility unit_test. Commutativity and associativity tests removed. Still raises NotImplemented error until validated with test_IR_limits

534. By Nicolas Deutschmann on 2018-08-28

Debugged the evaluation of the Initial-final soft-collinears. pp->Z soft limit validated. pp->Z collinear has issues with cuts.

535. By Valentin Hirschi on 2018-08-30

1. Fixed the definition of the mapping_variable `Q` in InitialLorentzOneMapping. @Simone: CHECK.
2. Fixed the soft-initial-collinear current so as to correctly define pC as 'pC = pC + pS' and use the `is_cut` function.
3. Improved the display of plots from test_IR_limits on a grid of subplots on a single page as opposed to multiple windows.
4. Tested that all limits and poles work on 'u u~ > z g'.

536. By Valentin Hirschi on 2018-08-30

1. Merged with the latest version of the trunk.

[Valentin Hirschi (valentin-hirschi)]

New contributions dubbed 'beam-factorization' contributions (BF) are introduced, together with their corresponding PDF counterterms.
These are generated automatically by inspecting existing contributions at lower orders, so that this part of the code is N^kLO ready already in principle.

The workflow of the __call__ and sigma functions of the ME7Integrands has been significantly modified/improved in order to have a more streamlined processing of the flavour mapping and back-propagation in ISR counterterms.
This improvement is realised with the new concept that the 'sigma' function no longer returns a single weight but instead a collection of 'ME7Event', each containing a list o weights for all the flavour structures it can adopt.

The application of the cuts and observables is then much more elegant as it is directly applied to the overall list of events generated for one specific set of input variables in the unit hypercube.

Given this large structural change, it is desirable to merge this branch sooner than later in the MadNkLO trunk. It is however not complete as it must still be tested beyond 2>1, at NNLO too, and the finite part of the integrated soft and soft-collinear counterterms is still missing.

This is however mostly irrelevant w.r.t this merge because the structural changes that were necessary for ISR are now established.

[lp:~madnklo/mg5amcnlo/MadNkLO_ISR]

537. By Valentin Hirschi on 2018-08-31

1. Improved and fixed the interface of MadNkLO to Cuba@Vegas, which is now a viable alternative to Vegas3.
2. Cuba@Vegas supports the same options as Vegas (--save_grids, --load_grids, etc...) plus some other.
3. Integrator options should now be supplied using the new command 'set_integrator_options <integrator_name> options'
4. When considering observables, the status for Cuba@Vegas is the same as for Vegas3: OK during refine only and without parallellization (which must therefore be done by hand)

538. By Valentin Hirschi on 2018-09-03

1. Improved pyCuba@Vegas to better handle saving and loading grids.
2. Improved handling of interrupts during parallel runs (however not satisfactory yet).

539. By Valentin Hirschi on 2018-09-03

1. Slight improvement of the KeyboardInterrupt handling.
2. Added to the pyCuba@Vegas grid output file the last integration result, to easily be able to recover it with max_eval==max_eval_survey==0.

540. By Valentin Hirschi on 2018-09-03

1. Fixed a bug introduced which prevented correct generation of chsi.

541. By Valentin Hirschi on 2018-09-06

1. First attempt at implementing BS contributions

542. By Valentin Hirschi on 2018-09-07

1. Implemented the BS contribution which is correctly generated when using the colorful scheme with a mapping recoiling soft against initial states.
2. Completed the implementation so that the BS integrand and its counterterms are correctly processed. This is running but not tested yet.
3. TODO: Adapt the phase-space generator so as to generate correlated chsi values.

543. By Valentin Hirschi on 2018-09-08

1. Finished the implementation of the ME7Integrand for the BS contribution and tested that test_IR_limits works for it.

544. By Valentin Hirschi on 2018-09-09

1. Finalized the implementation of the ME7_integrand_BS, including test_IR_poles and the integration.
2. Minor fixes for pyCuba@Vegas, still not perfect yet and as quirks.

545. By Valentin Hirschi on 2018-09-09

1. Fixed typo in last push

546. By Valentin Hirschi on 2018-09-09

1. Fix a typo in last push and solved a problem in the integration of the end-point beam-soft integrated CT.

547. By Valentin Hirschi on 2018-09-10

1. Fix on how to determine correlated beam currents for IntegratedBeamCurrents

548. By Valentin Hirschi on 2018-09-10

1. Corrected the spelling chsi->xi

549. By Nicolas Deutschmann on 2018-09-10

Added missing argument (walker) in definitions of scale_configuration

550. By Valentin Hirschi on 2018-09-10

1. Fixed so '\xi' to 'xi' as the slash was sometimes understood as an escape character.

551. By Valentin Hirschi on 2018-09-11

1. Finally, flavor mapping and mirrorring is now supported.
2. The problem of overcounting integrated ISR has been adressed.
3. In theory, all should now be working and 'p p > z QED=1 --NLO' is tested
   to be integrating. Next step is validation.

552. By Valentin Hirschi on 2018-09-12

1. Fixed the passing of options to the Vegas3 integrator

553. By Nicolas Deutschmann on 2018-09-13

Changed output format to JSON, default is now to append result to existing file, new launch option --run_name

554. By Valentin Hirschi on 2018-09-13

Fixed issues that prevented p p > z j to run fine. Now it does.

555. By Valentin Hirschi on 2018-09-13

Fixed some issues to generate NNLO ISR, e.g. p p > z --NNLO with option '--ignore_integrated_counterterms=RR' for the output

556. By Valentin Hirschi on 2018-09-14

Fixed severe bug in the subtraction of EpsilonExpansion

557. By Nicolas Deutschmann on 2018-09-17

added the integrated soft and soft collinear counterterms computed using the Soft vs IS mapping of Vittorio.

558. By Nicolas Deutschmann on 2018-09-17

Corrected the scaling of the matrix element when using Vittorio's Soft vs IS mapping

559. By Nicolas Deutschmann on 2018-09-18

Corrected a typo in colorful.NLO.beam_factorization : need to use filtered_flavor_matrix

560. By Nicolas Deutschmann on 2018-09-18

changes to ME7Event.convolve_flavor:
- doc reflects that a flavor matrix second entry (end_flavors) is actually a tuple to allow multiple flavors at once
- unspecified flavor matrix entries are set to 0, not 1

added ME7Event.is_empty : check if the convolution was with an empty matrix

in ME7Integrand.convolve_event_with_beam_factorization_currents:
- after convolution check if this has made the event empty, then return None as an event. This makes it be ignored

in test_IR_poles_for_process: now weights are summed, not events. This avoids issues with mapped events that have different kinematics (swapped leg orders) which cannot be summed.

561. By Nicolas Deutschmann on 2018-09-18

Refactored flavor masking in a separate function and corrected it (it did not account for end_flavors being tuples)

562. By Valentin Hirschi on 2018-09-21

1. BUG FIX: Real-emission local counterterms were not convoluted with the right PDFs.
2. BUG FIX: Sigma failed to add correctly the counterterms to the collection of events generated because of an indentation error.
3. Removed the ugly loop-induced specific class and refactored MadNkLO so that it can do loop-induced exactly like other processes.
4. Added options at generation stage to overwrite the process definition used for any contribution.
5. Improved MadLoop standalone loop-induced output so as to add the missing F2PY function hooks that the MadLoopAccessors needed.
6. Corrected how the flag 'requires_mirrorring' is set for integrated counterterms.
7. Fixed how the mirrored events are generated, so as to really swap the entire attribute relevant to both initial states.
8. Improved the printouts of test_IR_poles.
9. Fixed when the cuts (higher and lower) are called in the test_* functions depending on the options specified by the user.
Conclusion: This is the first revision that has a chance of yielding correct ISR NLO QCD cross-sections.

563. By Valentin Hirschi on 2018-09-21

1. Small fix in Vegas3 integrator options

564. By Valentin Hirschi on 2018-09-23

1. Fixed how the mirrored events of ISR collinear are generated as they need boosting since the maappinds sends them off the c.o.m. frame.

565. By Valentin Hirschi on 2018-09-24

1. Fixed how the flavour map is applied to the integrated collinear counterterms and [F] currents.

566. By Nicolas Deutschmann on 2018-09-24

Added dispatch in softs recoiling agaist IS according to dipole type (II, IF, FF). Corrected II integrated soft. Added sanity check to avoid handling massive emitters.

567. By Nicolas Deutschmann on 2018-09-24

corrected Bjorken x bound to avoid crashing when doing e+ e- annihilation. (x<=1 instead of x<1)

568. By Valentin Hirschi on 2018-09-24

1. Improved the printout of test_IR_poles
2. FIXED the PS point specifed when creating counter-events didn't have the right kinematics.
3. FIXED PS point of the counter-event *must* be boosted back to the com for distributions to be correct
   (this is because the xi rescalines must *not* apply with the change of variables chosen by MadNkLO to generate them)

569. By Nicolas Deutschmann on 2018-09-25

Corrected the printout when creating a new folder for results and improved comments on the result folder appending/creating logic

570. By Valentin Hirschi on 2018-09-28

1. Fixed support for the option --loop_filter
2. Fixed LO beam type guessing
3. Fixed implementation of cuts at LO
4. Fixed return 0. -> continue in __call__ of the integrand base class

571. By Valentin Hirschi on 2018-10-01

1. Fixed how to implement user-based cuts.

572. By Valentin Hirschi on 2018-10-01

1. Fixed a bug in the function that constructs the flavor cuts

573. By Valentin Hirschi on 2018-10-01

1. Fix bug in handling process options

574. By Valentin Hirschi on 2018-10-01

1. Trivial bug fixed in flavor cut building function

575. By Valentin Hirschi on 2018-10-01

1. Correctly apply Gabor's cut to initial state local collinear CT.
2. Compute Q directly in evaluate_(integrated)_counterterm and pass it to currents (possibly overwriting what was computed in the mappings).

576. By Valentin Hirschi on 2018-10-02

1. Also passed Q to the F-type of currents which is necessary structurally, even though they do not depend on it.

577. By Valentin Hirschi on 2018-10-02

1. Fixed the sum of initial momentum 'Q' passed to the PDF counterterms and integrated initial-state collinear counterterms so that they do not include the xi rescalings.

578. By Valentin Hirschi on 2018-10-02

1. Fixed the fact that the '+' distributions in the integrated ISR collinear CT must also apply on the 'Q^2/(xi1*xi2)' prefactor.
   This brings MadNkLO in full agreement with MadFKS for the process 'p p > z' ! :)

579. By Valentin Hirschi on 2018-10-03

1. Fixed last implementation of the log(Q^2)_plus so that it does not crash on endpoint evaluation.

580. By Valentin Hirschi on 2018-10-03

Typo in last push

581. By Valentin Hirschi on 2018-10-03

1. Fixed the definition of NF in the beam factorization currents.

582. By Nicolas Deutschmann on 2018-10-05

Added a new current scheme to separate colorful and colorful with Vittorio's mapping

583. By Valentin Hirschi on 2018-10-05

1. Improved printouts of the counterterm structure of the events generated, also in test_IR_poles.

584. By Valentin Hirschi on 2018-10-05

1. Safe-guarded the normalization factor in test_IR_poles.

585. By Nicolas Deutschmann on 2018-10-05

Inserted forgotten multiplicity factor in integrated counterterms

586. By Nicolas Deutschmann on 2018-10-05

Added the possibility to have 1,2,3 quarks in protons.
A warning is displayed when using a single quark: debug only

587. By Nicolas Deutschmann on 2018-10-05

Bug correction: identifying colorful_pp as a scheme using BS contributions and correcting selection of initial-initial eikonals

588. By Simone Lionetti <email address hidden> on 2018-10-05

Added safeguard for loading libcuba.so

589. By Valentin Hirschi on 2018-10-05

1. Fixed obvious mistake in the recursion setup for function get_reduced_flavor of the CountertermNode.

590. By Valentin Hirschi on 2018-10-07

1. Fixed issue with mapped processes having non-consecutive leg numbers.
2. Insure to completely skipped the beam factorization treatment when no contribution has an active beam.

591. By Nicolas Deutschmann on 2018-10-11

ProcessKey now always sorts final states PDGs. The sort_PDGs options only applies to the initial state.

592. By Simone Lionetti <email address hidden> on 2018-10-11

Improved error message of get_parent

593. By Simone Lionetti <email address hidden> on 2018-10-11

Partially fixed counterterm handling in test_IR_limits

594. By Simone Lionetti <email address hidden> on 2018-10-12

Improved handling of flavor assignments and cut events in test_IR_limits

595. By Simone Lionetti <email address hidden> on 2018-10-12

Commented out unused momenta from CataniSeymour currents that could cause fails

596. By Simone Lionetti <email address hidden> on 2018-11-26

Removed a printout of the soft PS point during mapping

597. By Simone Lionetti <email address hidden> on 2018-11-26

Implemented currents for the subtraction of 'e+ e- > u d u~ d~'
with distributed softs.

598. By Simone Lionetti <email address hidden> on 2018-11-26

Improvements to ME7 while coding qqQQ subtractions

These include:
- Saving of plots in grid format
- Print statements for debugging updated (and commented)
- ME-def line is not plotted by default (meaningless with reshuffled counterterms)
- Formatting and cleanup

599. By Simone Lionetti <email address hidden> on 2018-11-27

Fixed a bug in the combination of spin-correlations, some formatting

600. By Valentin Hirschi on 2018-11-28

1. Fixed spin-correlated MEs

601. By Valentin Hirschi on 2018-11-29

1. Forced the fastjet clustering algorithm to be genkt

602. By Valentin Hirschi on 2018-11-30

1. Improved some quirks in test_IR_limits

603. By Valentin Hirschi on 2018-11-30

1. Forced the absence of a factor in the initial collinear currents of colorful.

604. By Valentin Hirschi on 2018-11-30

1. Force the massive Q ref. collinear variables for colorful.

605. By Valentin Hirschi on 2018-11-30

1. Added a new attribute FROZEN_DIMENSIONS to the ME7Integrand base class that allows
   to freeze integration variables to a desired values so as to artificially lower the
   dimensionality of the phase-space integral to perform.

606. By Nicolas Deutschmann <email address hidden> on 2018-12-07

removed colorful cut and factors in local soft and SC counterterms

607. By Nicolas Deutschmann <email address hidden> on 2018-12-07

Temporary xi^2 correction factor in colorful_pp BS bulk to crrect flux factor. A more permanent change is warranted

608. By Nicolas Deutschmann <email address hidden> on 2018-12-07

Correct Q value in coloful_pp BS counterterm

609. By Nicolas Deutschmann <email address hidden> on 2018-12-10

Now using the right prefactor for colorful_pp beam_soft

610. By Valentin Hirschi on 2018-12-10

1. Improved the slowdown from the dot functions in the Vector classes.

611. By Valentin Hirschi on 2019-01-07

1. Added a minimal fix for not applying the Bjorken x's xi boundaries on the counterterm part of the beam convolution CTs. Will need to be generalised.

612. By Valentin Hirschi on 2019-01-14

1. Cleanly implemented the upper bound based on xi for the Bjorken x's in the convolution kernels that are not of type 'counterterms'.
2. Fixed the multiplicity of integrated ISR CT when doing flavor grouping.

613. By Valentin Hirschi on 2019-01-14

1. Merged with latest version of the MadNkLO trunk.

614. By Valentin Hirschi on 2019-01-14

1. Fixed all MadNkLO tests

=== modified file 'Template/ME7/Cards/run_card.dat'
--- Template/ME7/Cards/run_card.dat	2018-07-26 20:20:01 +0000
+++ Template/ME7/Cards/run_card.dat	2019-01-14 20:53:09 +0000
@@ -68,6 +68,8 @@
 #*********************************************************************
 # Standard Cuts                                                      *
 #*********************************************************************
+ %(flavor_cuts)s = flavor_cuts ! function string, or keyword: 'hardcoded'
+#*********************************************************************
 # Minimum and maximum pt's (for max, -1 means no cut)                *
 #*********************************************************************
  %(ptj)s  = ptj       ! minimum pt for the jets 

=== modified file 'madgraph/core/accessors.py'
--- madgraph/core/accessors.py	2018-05-25 08:29:19 +0000
+++ madgraph/core/accessors.py	2019-01-14 20:53:09 +0000
@@ -45,6 +45,20 @@
     
     cache_active = False
     
+    #######################################################################################
+    #  WARNING:
+    #  The squared orders and coupling orders (those that correspond to the perturbed couplings
+    #  of the reduced processes of counterterms is not correctly set for now. It is typically
+    # only relevant when performing mixed QCD-QED corrections which MadNkLO does not support
+    # yet. So for now the option below disables the order_matching in the process key for
+    # the orders specified in the variable below. Eventually however, all orders should be
+    # matched correctly, except 'WEIGHTED' maybe, which is typically only used for generation
+    # performance optimizations.
+    _DISABLED_ORDERS_FOR_PROCESS_MATCHING = ['QCD','QED','WEIGHTED']
+    #######################################################################################
+    # Disable process coupling orders-matching
+    
+    
     def get_key_for_cache(self,
                  process,  PDGs, sort_PDGs, vetoed_attributes, allowed_attributes, opts):
         """ Generates the look-up key for the dynamically created dictionary 'process_key_cache'
@@ -69,7 +83,7 @@
                 # intputs from this __init__
                 sort_PDGs = True,
                 # Exclude certain process attributes when creating the key for a certain process.
-                vetoed_attributes = ['model','legs','uid','has_mirror_process','legs_with_decays'],
+                vetoed_attributes = ['model','legs','uid','has_mirror_process','legs_with_decays','beam_factorization'],
                 # Specify only selected attributes to end up in the ProcessKey.
                 # If None, this filter is deactivated.
                 allowed_attributes = None,
@@ -114,7 +128,10 @@
                                      " via the options 'PDGs', 'legs' or 'process' (with precedence in this order).")
             
             if sort_PDGs:
-                self.key_dict['PDGs'] = ( tuple(sorted(self.key_dict['PDGs'][0])), tuple(sorted(self.key_dict['PDGs'][1])) )
+                self.key_dict['PDGs'] = ( tuple(sorted(self.key_dict['PDGs'][0])), tuple(sorted(self.key_dict['PDGs'][1])))
+            else: # The final state still needs to be sorted
+                self.key_dict['PDGs'] = ( self.key_dict['PDGs'][0], tuple(sorted(self.key_dict['PDGs'][1])))
+
         
         # Now make sure to instantiate a default process if the user didn't select one
         if not process:
@@ -175,12 +192,11 @@
             # overlap and cannot be constructed unambiguously in reduced processes for counterterms
             elif proc_attr in ['orders']:
                 self.key_dict[proc_attr] = hash_dict(dict((k, v) for k, v in value.items() 
-                                                             if k!='WEIGHTED'), proc_attr)
+                            if k not in self._DISABLED_ORDERS_FOR_PROCESS_MATCHING), proc_attr)
             elif proc_attr in ['sqorders_types', 'squared_orders']:
                 self.key_dict[proc_attr] = hash_dict(dict((k, v) for k, v in value.items() if \
-                        ( k!='WEIGHTED' and ( isinstance(process, subtraction.Current) or
-                                              (k not in process['sqorders_types']) or 
-                                              process['sqorders_types'][k]=='==') )
+                        ( k not in self._DISABLED_ORDERS_FOR_PROCESS_MATCHING and ( isinstance(process, subtraction.Current) or
+                        (k not in process['sqorders_types']) or process['sqorders_types'][k]=='==') )
                     ), proc_attr)
     
             elif proc_attr == 'parent_subtraction_leg':
@@ -192,7 +208,7 @@
                 self.key_dict[proc_attr] = hash_list(sorted(value), proc_attr)
     
             # Now generically treat other attributes
-            elif isinstance(value, (int, str, tuple, bool)):
+            elif value is None or isinstance(value, (int, str, tuple, bool)):
                 self.key_dict[proc_attr] = value
             elif isinstance(value, list):
                 self.key_dict[proc_attr] = hash_list(value, proc_attr)
@@ -884,7 +900,7 @@
     def call_subtraction_current(self, *args, **opts):
         """ Wrapper around the actual call of the subtraction currents, so as to be
         able to easily time it with a profiler."""
-        
+
         return self.subtraction_current_instance.evaluate_subtraction_current(*args, **opts)
 
     def __call__(self, current, **opts):
@@ -937,7 +953,7 @@
         # Update the cache with the new results produced
         if not cache_key is None:
             all_evaluations = self.cache.add_result(all_evaluations, **cache_key)
-        
+
         # Return both the specific evaluation asked for
         # and all results available for this cache_key
         return self.evaluation_class(evaluation_asked_for), self.result_class(all_evaluations)
@@ -979,14 +995,10 @@
 
         if cls is F2PYMEAccessor:
             target_class = None
-            n_loops = process.get('n_loops')
-            if process.get('perturbation_couplings') and n_loops==1:
+            if process.get('perturbation_couplings') and process.get('NLO_mode') in ['virt','sqrvirt']:
                 target_class = F2PYMEAccessorMadLoop
-            elif n_loops==0:
+            else:
                 target_class = F2PYMEAccessor
-            else:
-                raise MadGraph5Error("Could not determine the type of F2PYMEAccessor suited for "+
-                                     " %s with %d loops."%(process.nice_string(), n_loops))
 
             return super(F2PYMEAccessor, cls).__new__(target_class, 
                                                 process, f2py_module_path, slha_card_path, **opts)
@@ -1109,7 +1121,7 @@
             self.f2py_module = reload(self.f2py_module)
 
         # Now gather various properties about the Matrix Elements
-        # Get all helicities orderin
+        # Get all helicities ordering
         self.helicity_configurations = dict((tuple(hels),i+1) for i, hels in enumerate(
                                               self.get_function('get_helicity_definitions')() ))
         # Reversed map
@@ -1498,25 +1510,30 @@
 
     cache_active = False
 
-    def __init__(self, process, f2py_module_path, slha_card_path, madloop_resources_path=None, **opts):
+    def __init__(self, process, f2py_module_path, slha_card_path, 
+                                        madloop_resources_path=None, root_path='', **opts):
         """ Use the MadLoop resources path for MadLoop outputs """
   
-        super(F2PYMEAccessorMadLoop, self).__init__(process, 
-                                                 f2py_module_path, slha_card_path, **opts)
-        
         # Add the additional inputs to the back-up to be used later for the dump 
-        #  (if not already added by the mother)
-        if not 'madloop_resources_path' in self.initialization_inputs['opts']:
-            self.initialization_inputs['opts']['madloop_resources_path'] = madloop_resources_path
-
+        #  (if not already present). It is important that this path is set before the 
+        # call to the mother's __init__ implementation because this path must have been
+        # already set at the time synchronize() is called.
         if madloop_resources_path:
             # Save the location of MadLoop resources directory path
             if not os.path.isabs(madloop_resources_path):
                 raise MadGraph5Error("The initialization of F2PYMEAccessorMadLoop necessitates an "+
                                                             "absolute path for madloop_resources_path.")
-            self.madloop_resources_path = os.path.relpath(madloop_resources_path, self.root_path)
+            self.madloop_resources_path = os.path.relpath(madloop_resources_path, root_path)
         else:
             self.madloop_resources_path = None
+  
+        super(F2PYMEAccessorMadLoop, self).__init__(process, 
+                            f2py_module_path, slha_card_path, root_path=root_path, **opts)
+
+        # Add the input madloop_resouces_path to the initialization inputs necessary for the
+        # dump (if not already present).
+        if not 'madloop_resources_path' in self.initialization_inputs['opts']:
+            self.initialization_inputs['opts']['madloop_resources_path'] = madloop_resources_path
 
         if self.f2py_module is None:
             # We are at the generation stage and the f2py module has never been compiled yet.
@@ -1528,6 +1545,9 @@
         self.loop_squared_orders = {None: 0}
         if self.has_function('get_nsqso_loop'):
             n_squared_orders = self.get_function('get_nsqso_loop')()
+            # For loop-induced processes, the split order names would not have been set yet, so it is
+            # always OK to re-import it here.
+            self.split_order_names = [''.join(name).strip() for name in self.get_function('get_split_order_names')()]
             for i_sq_order in range(1,n_squared_orders+1):
                 key = tuple(sorted([(self.split_order_names[i], value) for i, value in enumerate(
                                             self.get_function('ml5get_squared_orders_for_soindex')(i_sq_order) )]))
@@ -1553,6 +1573,11 @@
         """ Synchronizes this accessor with the possibly updated value of parameter cards and ME source code.
         Must be defined by daughter classes."""
         
+        # Before the whole synchronization make sure that the current MadLoop resource path
+        # is correctly propagated to fortran.
+        if self.madloop_resources_path:
+            self.f2py_module.initialise_madloop_path(pjoin(self.root_path, self.madloop_resources_path))
+        
         ret_value = super(F2PYMEAccessorMadLoop, self).synchronize(ME7_options, from_init=from_init,**opts)
         
         # Never do anything when called from the __init__ function
@@ -1665,8 +1690,6 @@
         if not self.module_initialized:
             # These functions are not prefixed, so we should not ask to get them via the accessor self.get_function
             self.f2py_module.initialise(pjoin(self.root_path, self.slha_card_path))
-            if self.madloop_resources_path:
-                    self.f2py_module.initialise_madloop_path(pjoin(self.root_path, self.madloop_resources_path))
             self.module_initialized = True
         
         
@@ -1684,6 +1707,11 @@
                        'estimated_accuracies': estimated_accuracies,
                        'return_code': return_code}   
 
+        if __debug__ and return_code%100 ==4:
+            logger.debug("MadLoop library '%s@%s' could not reach target accuracy %.2e for PS_point:\n%s"%(
+                self.f2py_module_path[1], self.f2py_module_path[0],
+                required_accuracy, str(PS_point) ))
+            logger.debug("Accuracies obtained: %s"%str(estimated_accuracies))
         # Gather additional newly generated output_data to be returned and placed in the cache.
         output_datas = self.gather_output_datas(main_output, new_opts)
         
@@ -1936,18 +1964,21 @@
 
         return new_color_correlations
             
-    def format_spin_correlation(self, process, spin_correlations):
-        """ Synchronize the numbers in the spin correlation specifier to the leg_number in the process."""        
+    def format_spin_correlation(self, process, spin_correlation):
+        """Synchronize the leg numbers in the spin correlation specifier
+        to the leg numbers in the _reduced_ process.
+        """
         
         leg_number_to_pos_dict = {}
         for leg_pos, leg in enumerate(process.get_initial_legs()+process.get_final_legs()):
             leg_number_to_pos_dict[leg.get('number')] = leg_pos+1
-        
-        new_spin_correlations = []
-        for spin_correlation in spin_correlations:
-            new_spin_correlations.append( ( leg_number_to_pos_dict[spin_correlation[0]], spin_correlation[1] ) )
-
-        return new_spin_correlations
+
+        new_spin_correlation = []
+        for pol_vector in spin_correlation:
+            new_spin_correlation.append(
+                (leg_number_to_pos_dict[pol_vector[0]], pol_vector[1],))
+
+        return new_spin_correlation
         
     def __call__(self, *args, **opts):
         """ Quick access to directly calling a Matrix Element. The first argument should always be the MEdictionary key
@@ -2030,7 +2061,7 @@
                         args[0].accessor[key] = (ME_accessor, call_options)
                     else:
                         args[0].accessor = {key:(ME_accessor, call_options)}
-
+        
         return ME_accessor(*call_args, **call_options)
     
     def add_MEAccessor(self, ME_accessor, allow_overwrite=False):

=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2018-09-14 10:14:26 +0000
+++ madgraph/core/base_objects.py	2019-01-14 20:53:09 +0000
@@ -56,6 +56,7 @@
 
         for item in init_dict.keys():
             self.set(item, init_dict[item])
+            
         
 
     def __getitem__(self, name):
@@ -2813,6 +2814,15 @@
         # states are interchanged. This is only used for post-processing by ME7
         self['has_mirror_process'] = False
 
+        # The following will store the beam factorization currents that applies
+        # to this particular proces. The elements of the list have the following format
+        # [{ 'beam_one': <BeamCurrent>,
+        #    'beam_two': <BeamCurrent> }]
+        # while None in place of a current indicates no convolution.
+        # It may be necessary to have several instances of BeamCurrent, for example for the
+        # combination F1^(1)F2^(0) and F1^(0)F2^(1)
+        self['beam_factorization'] = [{'beam_one': None, 'beam_two': None},]
+
     def filter(self, name, value):
         """Filter for valid process property values."""
 
@@ -2838,6 +2848,21 @@
                     raise self.PhysicsObjectError, \
                           "%s is not a valid string" % str(value)
 
+        if name == 'beam_factorization':
+            if not isinstance(value, list):
+                raise self.PhysicsObjectError, \
+                        "%s is not a valid list" % str(value)
+            for bf in value:
+                if not isinstance(bf, dict) or \
+                   not set(bf.keys())==set(['beam_one','beam_two']):
+                    raise self.PhysicsObjectError, \
+                        "%s is not a valid list of dictionary with keys ['beam_one','beam_two']" % str(bf)
+                from madgraph.core.subtraction import BeamCurrent
+                for key in ['beam_one','beam_two']:
+                    if (not bf[key] is None) and (not isinstance(bf[key],BeamCurrent)):
+                        raise self.PhysicsObjectError, \
+                        "%s is not None or a valid instance of a BeamCurrent" % str(bf[key])
+
         if name == 'split_orders':
             if not isinstance(value, list):
                 raise self.PhysicsObjectError, \
@@ -2948,10 +2973,6 @@
             if value and isinstance(value, list) and \
                not isinstance(value[0], list):
                 value = [value]
-        
-        if name == 'perturbation_couplings':
-            if value and self['n_loops']==0:
-                self['n_loops']=1
 
         return super(Process, self).set(name, value) # call the mother routine
 
@@ -2989,7 +3010,8 @@
                 'forbidden_onsh_s_channels', 'forbidden_s_channels',
                 'forbidden_particles', 'is_decay_chain', 'decay_chains',
                 'legs_with_decays', 'perturbation_couplings', 'has_born', 
-                'NLO_mode','split_orders','n_loops', 'has_mirror_process']
+                'NLO_mode','split_orders','n_loops', 'has_mirror_process',
+                'beam_factorization']
 
     def nice_string(self, indent = 0, print_weighted = True, prefix = True):
         """Returns a nicely formatted string about current process
@@ -3117,12 +3139,33 @@
         if self['has_mirror_process']:
             mystr = mystr + ' (mirrored)'
         
-        if not self.get('decay_chains'):
-            return mystr
+        if self.get('decay_chains'):
+            for decay in self['decay_chains']:
+                mystr = mystr + '\n' + \
+                        decay.nice_string(indent + 2).replace('Process', 'Decay')         
 
-        for decay in self['decay_chains']:
-            mystr = mystr + '\n' + \
-                    decay.nice_string(indent + 2).replace('Process', 'Decay')         
+        # Finally add beam factorization
+        if len(self['beam_factorization'])>1 or \
+           ( len(self['beam_factorization'])>0 and
+             ( (not self['beam_factorization'][0]['beam_one'] is None) or \
+               (not self['beam_factorization'][0]['beam_two'] is None)) ):
+            beam_factorization_terms_str = []
+            for bft in self['beam_factorization']:
+                if (not bft['beam_one'] is None) and \
+                   (not bft['beam_two'] is None):
+                    parenthesis = '{%s}'
+                else:
+                    # I find it more esthetical to always include the curly brackets
+                    parenthesis = '{%s}'
+                beam_factorization_terms_str.append(parenthesis%('|'.join([
+                    bc.__str__(print_n=True, print_pdg=False, print_state=False,
+                               print_loops=True, print_orders=False) for bc in 
+                    [bft['beam_one'],bft['beam_two']] if not bc is None 
+                ])))
+            if len(self['beam_factorization'])>1:
+                mystr += ' x [%s]'%(' + '.join(beam_factorization_terms_str))
+            else:
+                mystr += ' x %s'%(' + '.join(beam_factorization_terms_str))
 
         return mystr
 
@@ -3411,6 +3454,35 @@
             self.initial_final_pdgs = (initial_pdgs, final_pdgs)
             return self.initial_final_pdgs
 
+    def get_n_loops_in_beam_factorization(self):
+        """ Returns the total number of loops in each of the two beam factorization currents,
+        this sum is expected to be identical for each beam factorization term."""
+        
+        bft = self['beam_factorization'][0]
+        n_loops = 0
+        for beam_name in ['beam_one','beam_two']:
+            if not bft[beam_name] is None:
+                n_loops += bft[beam_name]['n_loops']
+
+        return n_loops
+        
+    def get_squared_orders_in_beam_factorization(self):
+        """ Returns the total squared orders in each of the two beam factorization currents,
+        this sum is expected to be identical for each beam factorization term."""
+        
+        bft = self['beam_factorization'][0]
+        squared_orders = {}
+        
+        for beam_name in ['beam_one','beam_two']:        
+            if not bft[beam_name] is None:
+                for squared_order, value in bft[beam_name]['squared_orders'].items():
+                    if squared_order in squared_orders:
+                        squared_orders[squared_order] += value
+                    else:
+                        squared_orders[squared_order] = value
+
+        return squared_orders
+
     def format_PS_point_for_ME_call(self, PS_point):
         """ From a dictionary formatted PS point and a process, returns the PS point as a flat list, ordered as
         the legs in the process."""
@@ -3542,7 +3614,8 @@
         
         return ( tuple( l.get('id') for l in self.get('legs') if not l.get('state') ), 
                  tuple( l.get('id') for l in self.get('legs') if l.get('state') ),
-                 self.get('id'), self.get('n_loops'), tuple(self.get('perturbation_couplings')) )
+                 self.get('id'), self.get('n_loops'), tuple(self.get('perturbation_couplings')),
+                 tuple( (str(bf['beam_one']),str(bf['beam_two'])) for bf in self.get('beam_factorization')) )
 
     def compare_for_sort(self, other):
         """Sorting routine which allows to sort processes for
@@ -3937,7 +4010,8 @@
             'split_orders': self.get('split_orders'),
             'NLO_mode': self.get('NLO_mode'),
             'n_loops': self.get('n_loops'),
-            'has_mirror_process': self.get('has_mirror_process')
+            'has_mirror_process': self.get('has_mirror_process'),
+            'beam_factorization': self.get('beam_factorization'), 
             })
             
     def get_process(self, initial_state_ids, final_state_ids):
@@ -4016,7 +4090,13 @@
                  correction_couplings       = [], 
                  n_unresolved_particles     = 0,
                  n_loops                    = -1,
-                 squared_orders_constraints = {}
+                 squared_orders_constraints = {},
+                 beam_types                 = (None, None),
+                 # Active simply means the a convolution is in place.
+                 beam_factorization_active  = (False, False),
+                 # Correlated beam convolutions means that both beams are convoluted with
+                 # the same variable xi
+                 correlated_beam_convolution = False
                 ):
         """ Instantiate a contribution definition with all necessary information
         to generate the actual contribution. """
@@ -4028,8 +4108,47 @@
         self.correction_order          = correction_order
         self.correction_couplings      = correction_couplings
         self.n_unresolved_particles    = n_unresolved_particles
-        # Make sure that n_loops matches the one in the process definition. We however keep
-        # this attribute in ContributionDefinition as well for convenience.
+        
+        # Values of the dictionary 'beam_factorization' are dicts of the form:
+        #  {  'beam_type' : beam_type_name_str
+        #     'beam_PDGs' : (beam_PDGs,)
+        #  }
+        self.beam_factorization = {}
+        for i, beam_name in enumerate(['beam_one', 'beam_two']):
+            if beam_types[i] is None:   
+                self.beam_factorization[beam_name] = None
+            else:
+                self.beam_factorization[beam_name] = {
+                    'beam_type'  : beam_types[i][0],
+                    'beam_PDGs'  : tuple(beam_types[i][1]),
+                    # Active means that a convolution with beam factorization is in place.
+                    'active'     : beam_factorization_active[i]
+                }
+        if correlated_beam_convolution and beam_factorization_active != (True, True):
+            raise MadGraph5Error("The beam convolution can be specified to be correlated only if"+
+                                 " both beams are convoluted.")
+        self.correlated_beam_convolution = correlated_beam_convolution
+
+        # Extract n_loops from the one in the process definition. However these need *not*
+        # be the same, for example a Loop-Induced contribution 'RV' would have a process 
+        # definition with n_loops=2 while the contribution itself has n_loops=1 because it
+        # must always match the expected number of loops in such contribution (i.e. RV in this
+        # example) in the case of a non-loop induced process.
+        if n_loops == -1:
+            self.n_loops               = process_definition.get('n_loops')
+        else:
+            self.n_loops = n_loops 
+
+        # Possibly extract n_loops from the one in the process definition. Notice that the
+        # attribute 'n_loops' of the contribution and of the process *must* be the same 
+        # even in the case of loop-induced processes. This attribute 'n_loops' is meant to
+        # be used only for the whole MadNkLO construction where it *must* be equivalent to 
+        # labeling the number of loops one would find in this contribution or the processes
+        # hosted *in the situation where the lowest order is tree-level*.
+        # In order to specify the number of loops the process must really have when its 
+        # ME code is exported on disk, one must either set the 'NLO_mode' of the process
+        # to be 'virt' or 'sqrvirt', or use the squared order constraint 'NLOOP' when using
+        # multi-loop UFO form factors.
         if n_loops == -1:
             self.n_loops               = process_definition.get('n_loops')
         else:
@@ -4038,12 +4157,60 @@
                   "should be instantiated with a ProcessDefinition with the same n_loops.")
             else:
                 self.n_loops = n_loops 
+
         # Squared orders to consider. Note that this is not the same as the attribute
         # of the process definition, since it can be {'QED':[2,4]} to indicates that all
         # squared order for QED between 2 and 4 must be considered. A constraint of 
         # {'QED':[2,2]} means that QED^2 must be exactly equal to 2.
         self.squared_orders_constraints = squared_orders_constraints
     
+    def get_copy(self):
+        """ Returns a partially-shallow copy of self."""
+        return ContributionDefinition(
+                process_definition          = self.process_definition.get_copy(), 
+                overall_correction_order    = self.overall_correction_order,
+                correction_order            = self.correction_order,
+                correction_couplings        = self.correction_couplings, 
+                n_unresolved_particles      = self.n_unresolved_particles,
+                n_loops                     = self.n_loops,
+                squared_orders_constraints  = dict(self.squared_orders_constraints),
+                beam_factorization          = copy.deepcopy(self.beam_factorization),
+                correlated_beam_convolution = self.correlated_beam_convolution
+            )
+    
+    def get_beam_types(self):
+        """ Return the beam types in a tuple format."""
+
+        return ( 
+            ( self.beam_factorization['beam_one']['beam_type'],
+              self.beam_factorization['beam_one']['beam_PDGs']) if
+            (not self.beam_factorization['beam_one'] is None) else None,
+            ( self.beam_factorization['beam_two']['beam_type'],
+              self.beam_factorization['beam_two']['beam_PDGs']) if
+            (not self.beam_factorization['beam_two'] is None) else None )
+
+    def has_beam_factorization(self):
+        """ Checks whether this contribution has as any beam factorization active."""
+
+        return ( ( (not self.beam_factorization['beam_one'] is None) and 
+                   self.beam_factorization['beam_one']['active'] ) or
+                 ( (not self.beam_factorization['beam_two'] is None) and 
+                   self.beam_factorization['beam_two']['active'] ) )
+
+    def is_beam_active(self, beam_name):
+        """ Checks whether factorization of the beam 'beam_name' is active or not."""
+        return (not self.beam_factorization[beam_name] is None) and \
+                                              self.beam_factorization[beam_name]['active']
+
+    def is_loop_induced(self):
+        """ Returns whether this contribution hosts a 'loop-induced' process. This is not
+        so important in the whole MadNkLO construction and for now we simply use the 
+        mode of the process definition, but later we'll also have to investigate the 
+        coupling orders (to see if any is 'NLOOP') as (multi-)loops in Matrix Element can also 
+        be obtained with UFO form factors."""
+        n_process_loops = 1 if (self.process_definition.get('NLO_mode') in ['virt','sqrvirt']) else 0
+        return ((n_process_loops - self.n_loops) > 0)
+
     def nice_string(self):
         """ Nice representation of a contribution definition."""
         res = []
@@ -4054,7 +4221,48 @@
         res.append('%-30s:   %s'%('n_loops',self.n_loops))
         res.append('%-30s:   %s'%('squared_orders_constraints',self.squared_orders_constraints))
         res.append('%-30s:   %s'%('process_definition',self.process_definition.nice_string().replace('Process: ','')))
+        res.append('%-30s:   %s'%('beam_one',self.get_beam_factorization_nice_string(
+                self.process_definition['model'],self.beam_factorization['beam_one'])))
+        res.append('%-30s:   %s'%('beam_two',self.get_beam_factorization_nice_string(
+                self.process_definition['model'],self.beam_factorization['beam_two'])))
         return '\n'.join(res)
+
+    def get_beam_factorization_nice_string(self, model, beam_factorization_specification):
+        """ Returns a nice representation of the beam factorization specification in 
+        argument"""
+        
+        
+        if beam_factorization_specification is None:
+            return ''
+        
+        res = beam_factorization_specification['beam_type']
+
+        if beam_factorization_specification['active']:
+            if self.correlated_beam_convolution:
+                res += ' (correlated convolution)'
+            else:
+                res += ' (convoluted)'
+
+        res += ' [including: %s]'%(' '.join(model.get_particle(pdg).get_name() for pdg 
+                                     in beam_factorization_specification['beam_PDGs']))
+        return res
+    
+    def get_beam_factorization_short_name(self):
+        """ Return a short string identifying the beam factorization applied to this
+        contribution. For instance it can be:
+           F1      --> PDF CT applied to beam one
+           F2      --> PDF CT applied to beam two
+           F1F2    --> PDF CT applied to both beams
+           S       --> Correlated convolution of both beams
+        """
+        if self.correlated_beam_convolution:
+            return 'S'
+        res = ''
+        for i_beam, beam_name in enumerate(['beam_one', 'beam_two']):
+            if not self.beam_factorization[beam_name] is None and \
+                self.beam_factorization[beam_name]['active']:
+                res += 'F%d'%(i_beam+1)
+        return res        
     
     def short_name(self):
         """ Returns a short-hand notation for that contribution."""
@@ -4062,6 +4270,8 @@
                     'V'*self.n_loops
         if shortname=='':
             shortname = 'B'
+        # Add factorization labels
+        shortname += self.get_beam_factorization_short_name()
         return shortname
 
     def get_shell_name(self):
@@ -4079,6 +4289,9 @@
                 shell_name += 'B'
             else:
                 shell_name += 'V'*nV_in_conjugated_amplitude
+            beam_factorization_label = self.get_beam_factorization_short_name()
+            if beam_factorization_label != '':
+                shell_name += '_%s'%beam_factorization_label
         return shell_name
 
 class EpsilonExpansion(dict):
@@ -4152,6 +4365,11 @@
     __truediv__ = __div__
     
     def __radd__(self, other):
+        # The check below is useful for this overloaded __radd__ to work with the Python
+        # intrinsic sum() function.
+        if other==0:
+            return
+        misc.sprint(other)
         for k, v in other.items():
             try:
                 self[k] += v
@@ -4163,7 +4381,7 @@
             try:
                 self[k] -= v
             except KeyError:
-                self[k] = -v        
+                self[k] = -v   
     
     def __rmul__(self, other):
         if not isinstance(other, EpsilonExpansion):
@@ -4191,7 +4409,26 @@
             raise NotImplementedError
         
     __rtruediv__ = __rdiv__
-    
+
+    def get_term(self, term_specifier):
+        """ Get on particular term of the expansion, specified either by a string or by
+        an integer."""
+        dict_key = 0
+        term_specifier = term_specifier.lower()
+        if isinstance(term_specifier, str):
+            if term_specifier == 'sum_all':
+                return sum(v for k,v in self.items() if isinstance(k, int))
+            if term_specifier == 'finite':
+                dict_key = 0
+            elif term_specifier.startswith('eps^'):
+                dict_key = int(term_specifier[4:])
+            else:
+                raise MadGraph5Error('Incorrect specification of the epsilon expansion term to access: %s'%term_specifier)
+        else:
+            dict_key = term_specifier
+
+        return self[dict_key]
+
     def to_human_readable_dict(self):
         """ Transforms this expansion into a dictionary with human readable keys.
         Typically to be used in MEEvaluation results dictionaries. """

=== modified file 'madgraph/core/contributions.py'
--- madgraph/core/contributions.py	2018-06-16 19:41:18 +0000
+++ madgraph/core/contributions.py	2019-01-14 20:53:09 +0000
@@ -50,6 +50,14 @@
 
 logger = logging.getLogger('contributions')
 
+# The initial state collinear counterterms may be combined because the
+# BeamCurrents that implements them return the full backward evolution flavor
+# matrix, so that several initial state collinear counterterms that differ only
+# by the resulting backward evolved flavor (e.g. 'd g <- d' and 'g d~ <- d') could
+# be combined. This option is not implemented as of now, so the only valid value
+# for the global parameter below is False.
+_COMBINE_INITIAL_STATE_INTEGRATED_COLLINEAR_COUNTERTERMS = False
+
 #===============================================================================
 # Contribution mother class
 #===============================================================================
@@ -64,35 +72,39 @@
 
         if cls is Contribution:
             target_type = 'Unknown'
-            if contribution_definition.correction_order == 'LO':
-                if contribution_definition.n_loops == 0 and \
-                  contribution_definition.n_unresolved_particles == 0:
-                    target_type = 'Born'
-                elif contribution_definition.n_loops == 1 and \
-                  contribution_definition.n_unresolved_particles == 0:
-                    target_type = 'LoopInduced_Born'                    
-            elif contribution_definition.correction_order == 'NLO':
-                if contribution_definition.n_loops == 1 and \
-                   contribution_definition.n_unresolved_particles == 0:
+            correction_order = contribution_definition.correction_order.count('N')
+            beam_factorization_count = 0
+            if contribution_definition.is_beam_active('beam_one'):
+                beam_factorization_count += 1
+            if contribution_definition.is_beam_active('beam_two'):
+                beam_factorization_count += 1
+            n_loops = contribution_definition.n_loops
+            n_unresolved_particles = contribution_definition.n_unresolved_particles
+            # Beam factorization contributions are automatically of type RV because
+            # they must both generate local counterterms (for the form factors) and
+            # accept integrated ISR ones.
+            if beam_factorization_count > 0:
+                if contribution_definition.correlated_beam_convolution:
+                    target_type = 'BeamSoft'
+                else:
+                    target_type = 'RealVirtual'                
+            elif n_loops == 0 and n_unresolved_particles == 0:
+                target_type = 'Born'
+            elif n_loops == 1 and n_unresolved_particles == 0:
+                if correction_order < 1:
+                    target_type = 'LoopInduced_Born'
+                else:
                     target_type = 'Virtual'
-                elif contribution_definition.n_loops == 0 and \
-                     contribution_definition.n_unresolved_particles == 1:
-                    target_type = 'SingleReals'
-            elif contribution_definition.correction_order == 'NNLO':
-                if contribution_definition.n_loops == 0 and \
-                   contribution_definition.n_unresolved_particles == 2:
-                    target_type = 'DoubleReals'
-                else:
-                    raise MadGraph5Error("Some NNLO type of contributions are not implemented yet.")
-            elif contribution_definition.correction_order == 'NNNLO':
-                if contribution_definition.n_loops == 0 and \
-                   contribution_definition.n_unresolved_particles == 3:
-                    target_type = 'TripleReals'
-                else:
-                    raise MadGraph5Error("Some NNNLO type of contributions are not implemented yet.")
+            elif n_loops == 0 and n_unresolved_particles == 1:
+                target_type = 'SingleReals'
+            elif n_loops == 0 and n_unresolved_particles == 2:
+                target_type = 'DoubleReals'
+            elif n_loops == 0 and n_unresolved_particles == 3:
+                target_type = 'TripleReals'                
+            else:
+                raise MadGraph5Error("Some %s type of contributions are not implemented yet."%
+                                                  contribution_definition.correction_order)
 
-            else:
-                target_type = 'Unknown'
             target_class = Contribution_classes_map[target_type]
             if not target_class:
                 raise MadGraph5Error("Could not determine the class for contribution of type '%s' to be added for"%target_type+
@@ -147,15 +159,21 @@
         
         # Initialize an IR subtraction module if necessary
         self.IR_subtraction = None
-        if self.contribution_definition.n_unresolved_particles > 0:
+
+        if self.contribution_definition.n_unresolved_particles > 0 or \
+                                     self.contribution_definition.has_beam_factorization():
             self.IR_subtraction = subtraction.IRSubtraction(
                 self.model,
-                coupling_types=self.contribution_definition.correction_couplings,
-                n_unresolved=self.contribution_definition.n_unresolved_particles )
+                coupling_types  = self.contribution_definition.correction_couplings,
+                n_unresolved    = self.contribution_definition.n_unresolved_particles,
+                beam_types      = self.contribution_definition.get_beam_types(),
+                currents_scheme = self.options['subtraction_currents_scheme'],
+                mappings_scheme = self.options['subtraction_mappings_scheme'] )
         
-        # The following two attributes dicate the type of Exporter which will be assigned to this contribution
+        # The following two attributes dictate the type of Exporter which will be assigned to this contribution
         self.output_type             = 'default'
         self.format                  = 'standalone'
+        
         self.export_dir              = 'None'
         if self.options['group_subprocesses'] == 'Auto':
             self.collect_mirror_procs = True
@@ -166,22 +184,98 @@
 
         # Options relevant only for LO diagram generation
         self.ignore_six_quark_processes = False
-        self.diagram_filter             = False 
+        if 'diagram_filter' in opts:
+            self.diagram_filter         = opts['diagram_filter']
+        else:
+            self.diagram_filter         = False
         self.optimize                   = False
         # Flag to indicate whether this contribution supports decay chains
         # Typically only LO contributions support decay chain amplitudes
         self.supports_decay_chain       = False
         
+        if 'loop_filter' in opts:
+            self.loop_filter            = opts['loop_filter']
+        else:
+            self.loop_filter            = None
+
         # Specifies if this contribution output requires its own MODEL in Source
         self.requires_its_own_model     = False
         
+        # The following two attributes dictate the type of Exporter which will be assigned to this contribution
+        self.output_type             = self.get_output_type()
+        self.format                  = self.get_output_format()
+        
         # Specify the MultiProcessClass to use to generate amplitudes
-        self.MultiProcessClass          = diagram_generation.MultiProcess
+        self.MultiProcessClass       = self.get_multi_process_class()
         
         # Add default phase-space topology specifiers
-        self.processes_to_topologies = None
-        self.topologies_to_processes = None
-        
+        self.processes_to_topologies = {}
+        self.topologies_to_processes = {}
+        
+        # Add an empty processes map, specifying that it is obtained before anything
+        # was generated yet
+        self.processes_map = ({
+                'had_amplitudes'        :   False,
+                'had_matrix_elements'   :   False
+            }, {})
+        
+    def get_multi_process_class(self):
+        """ Function returning the appropriate multiprocess class to instantiate during
+        the amplitude generation for this particular contribution."""
+        
+        if self.output_type == 'madloop':
+            # Make sure to adjust the MultiProcessClass to be used in the case of a 'madloop' output
+            return loop_diagram_generation.LoopMultiProcess
+        elif self.output_type == 'default':
+            return diagram_generation.MultiProcess
+        else:
+            raise MadGraph5Error(
+                "Output type '%s' not reckognized by the base Contribution class."%self.output_type+
+                " It must be implemented explicitely in the daughter class.")
+
+    def get_output_type(self):
+        """ Function returning the appropriate output type to be used when outputting
+        the amplitudes from this particular contribution."""
+        
+        # Decide whether a MadLoop output is required based on the process definition
+        # underlying this contribution.
+        if (self.contribution_definition.process_definition.get('NLO_mode') in ['virt','sqrvirt']):
+            return 'madloop'
+        else:
+            return 'default'
+
+    def get_output_format(self):
+        """ Function returning the appropriate output format to be used when outputting
+        the amplitudes from this particular contribution."""
+        return 'standalone'
+        
+    def add_content_to_global_ME7_resources(self, global_ME7_dir, **opts):
+        """ Possibly add content to the global ME7 resources pool, depending on the output 
+        format requested."""
+        
+        if self.format == 'standalone' and self.output_type == 'madloop':        
+            destination = pjoin(global_ME7_dir,'Cards','MadLoopParams.dat')
+            if not os.path.exists(destination):
+                shutil.copyfile(pjoin(self.export_dir,'Cards','MadLoopParams.dat'), destination)
+
+
+    def set_helas_model(self):
+        """ Instantiate properly the helas model according to the exporter and the
+        output type."""
+
+        if self.exporter.exporter == 'cpp':       
+            self.helas_model = helas_call_writers.CPPUFOHelasCallWriter(self.model)
+        else:
+            assert self.exporter.exporter == 'v4'
+            if self.output_type == 'madloop':
+                assert (not self.options['_model_v4_path'])
+                self.helas_model = helas_call_writers.FortranUFOHelasCallWriter(self.model)
+            else:
+                if self.options['_model_v4_path']:
+                    self.helas_model = helas_call_writers.FortranHelasCallWriter(self.model)
+                else:
+                    self.helas_model = helas_call_writers.FortranUFOHelasCallWriter(self.model)
+
     def add_integrated_counterterm(self, integrated_CT_properties):
         """By default, do not support adding integrated counterterms."""
         
@@ -189,6 +283,10 @@
             "The contribution of type %s cannot receive" % type(self) +
             " contributions from integrated counterterms." )
 
+    def combine_initial_state_counterterms(self, *args, **opts):
+        """ By default, does nothing. For inheritance purposes only."""
+        pass
+
     def set_export_dir(self, prefix):
         """Assigns an export directory name."""
 
@@ -237,7 +335,42 @@
         return self.exporter.pass_information_from_cmd(cmd_interface)
 
     def generate_matrix_elements(self, group_processes=True):
-        """Generate the Helas matrix elements before exporting. Uses the main function argument 
+        """ Generate the matrix element in the manner appropriate to the selected output type."""
+        
+        if self.output_type == 'madloop':
+            return self.generate_loop_matrix_elements(group_processes=group_processes)
+        else:
+            return self.generate_tree_matrix_elements(group_processes=group_processes)
+    
+    def generate_loop_matrix_elements(self, group_processes=True):
+        """Generate the *loop* Helas matrix elements before exporting. Uses the main function argument 
+        'group_processes' to decide whether to use group_subprocess or not."""
+
+        cpu_time1 = time.time()
+        ndiags = 0
+
+        generation_mode = {'optimized_output': self.options['loop_optimized_output']}
+
+        self.all_matrix_elements = loop_helas_objects.LoopHelasProcess(
+            self.amplitudes,
+            compute_loop_nc = True,
+            matrix_element_opts = generation_mode,
+            optimized_output = self.options['loop_optimized_output'],
+            combine_matrix_elements=group_processes
+        )
+        ndiags = sum([len(me.get('diagrams')) for me in self.all_matrix_elements.get_matrix_elements()])
+        
+        # assign a unique id number to all process
+        uid = 0 
+        for me in self.all_matrix_elements.get_matrix_elements():
+            uid += 1 # update the identification number
+            me.get('processes')[0].set('uid', uid)
+
+        cpu_time2 = time.time()
+        return ndiags, cpu_time2 - cpu_time1     
+
+    def generate_tree_matrix_elements(self, group_processes=True):
+        """Generate the *tree* Helas matrix elements before exporting. Uses the main function argument 
         'group_processes' to decide whether to use group_subprocess or not."""
 
         cpu_time1 = time.time()
@@ -300,39 +433,51 @@
 
         return ndiags, cpu_time2 - cpu_time1
 
-    def set_helas_model(self):
-        """ Instantiate properly the helas model """
-
-        if self.exporter.exporter == 'cpp':       
-            self.helas_model = helas_call_writers.CPPUFOHelasCallWriter(self.model)
-        elif self.options['_model_v4_path']:
-            assert self.exporter.exporter == 'v4'
-            self.helas_model = helas_call_writers.FortranHelasCallWriter(self.model)
-        else:
-            assert self.exporter.exporter == 'v4'
-            self.helas_model = helas_call_writers.FortranUFOHelasCallWriter(self.model)
-
     def generate_code(self):
         """ Assuming the Helas Matrix Elements are now generated, we can write out the corresponding code."""
 
         matrix_elements = self.all_matrix_elements.get_matrix_elements()
-        
         calls=0
         cpu_time_start = time.time()
-        if isinstance(self.all_matrix_elements, group_subprocs.SubProcessGroupList) and \
-                                                                                self.exporter.grouped_mode:
-            modified, self.all_matrix_elements = self.exporter.modify_grouping(self.all_matrix_elements)
-            for me_number, me in enumerate(self.all_matrix_elements):
-                calls = calls + self.exporter.generate_subprocess_directory(me, self.helas_model, me_number)
-
-        else: # Non-grouped mode
-            for nb, me in enumerate(matrix_elements[:]):
-                new_calls = self.exporter.generate_subprocess_directory(me, self.helas_model, nb)
-                if isinstance(new_calls, int):
-                    if new_calls ==0:
-                        matrix_elements.remove(me)
-                    else:
-                        calls = calls + new_calls
+        
+        if self.output_type == 'madloop':
+            for me in matrix_elements:
+                # Choose the group number to be the unique id so that the output prefix
+                # is nicely P<proc_id>_<unique_id>_
+                calls = calls + self.exporter.generate_loop_subprocess(me, 
+                    self.helas_model,
+                    group_number = me.get('processes')[0].get('uid'),
+                    proc_id = None,
+                    config_map=None,
+                    unique_id=me.get('processes')[0].get('uid'))
+    
+            # If all ME's do not share the same maximum loop vertex rank and the
+            # same loop maximum wavefunction size, we need to set the maximum
+            # in coef_specs.inc of the HELAS Source. The SubProcesses/P* directory
+            # all link this file, so it should be properly propagated
+            if self.options['loop_optimized_output'] and len(matrix_elements)>1:
+                max_lwfspins = [m.get_max_loop_particle_spin() for m in \
+                                                                matrix_elements]
+                max_loop_vert_ranks = [me.get_max_loop_vertex_rank() for me in \
+                                                                matrix_elements]
+                if len(set(max_lwfspins))>1 or len(set(max_loop_vert_ranks))>1:
+                    self.exporter.fix_coef_specs(max(max_lwfspins),\
+                                                       max(max_loop_vert_ranks))
+        else:
+            if isinstance(self.all_matrix_elements, group_subprocs.SubProcessGroupList) and \
+                                                                                    self.exporter.grouped_mode:
+                modified, self.all_matrix_elements = self.exporter.modify_grouping(self.all_matrix_elements)
+                for me_number, me in enumerate(self.all_matrix_elements):
+                    calls = calls + self.exporter.generate_subprocess_directory(me, self.helas_model, me_number)
+    
+            else: # Non-grouped mode
+                for nb, me in enumerate(matrix_elements[:]):
+                    new_calls = self.exporter.generate_subprocess_directory(me, self.helas_model, nb)
+                    if isinstance(new_calls, int):
+                        if new_calls ==0:
+                            matrix_elements.remove(me)
+                        else:
+                            calls = calls + new_calls
         
         return calls, time.time() - cpu_time_start
 
@@ -344,7 +489,18 @@
         for amplitude in self.amplitudes:
             max_order = max(max_order, amplitude.get('diagrams').get_max_order(order))
         return max_order
+    
+    def subtract(self, ignore_integrated_counterterms=False, **opts):
+        """ Generate and export all necessary subtraction counterterterms and currents
+        (including beam factorisation counterterms as well [bulk ones are already generated
+        at this stage] ones.). Nothing to do in this base class."""
+        
+        # Assign an empty list of countrerterms for now to this contribution
+        self.counterterms = {}
 
+        # No information need to be returned at this stage since nothing was generated
+        return {}
+        
     def export(self, nojpeg=False, group_processes=True, args=[], **opts):
         """ Perform the export duties, that include generation of the HelasMatrixElements and 
         the actual output of the matrix element code by the exporter."""
@@ -384,7 +540,7 @@
         # Now investigate phase-space topologies and identify the list of kinematic configurations present
         # and, for each process key in the process map, what configuration it includes and with which defining diagrams.
         self.set_phase_space_topologies()
-        
+
         # The printout below summarizes the topologies identified and the information gathered about the
         # matrix elements and their diagrams
 #        misc.sprint('='*50)
@@ -550,14 +706,16 @@
         """ Adds all MEAccessors for the matrix elements generated as part of this contribution."""
         
         MEAccessors = []
+        
+        # This is not needed for factorization contributions as they always use processes already
+        # in other contributions.
+        if self.has_beam_factorization():
+            return MEAccessors
+
         for process_key, (defining_process, mapped_processes) in self.get_processes_map().items():
             # The PDGs of the hashed representations correspond to entry [0][0]
-            mapped_process_pdgs = [
-                (proc.get_initial_ids(), proc.get_final_ids())
-                for proc in mapped_processes ]
-            proc_dir = pjoin(
-                self.export_dir, 'SubProcesses',
-                'P%s' % self.process_dir_name(defining_process) )
+            mapped_process_pdgs = [ (proc.get_initial_ids(), proc.get_final_ids()) for proc in mapped_processes ]
+            proc_dir = pjoin(self.export_dir, 'SubProcesses', 'P%s' % self.process_dir_name(defining_process) )
             if not os.path.isdir(proc_dir):
                 raise MadGraph5Error(
                     "Cannot find the process directory '%s' for process %s." %
@@ -588,6 +746,46 @@
         all_MEAccessors.add_MEAccessors(
             MEAccessors, allow_overwrite=(not self.group_subprocesses) )
 
+    def add_current_accessors(
+        self, model, all_MEAccessors, root_path, current_set, currents_to_consider):
+        """Generate and add all subtraction current accessors to the MEAccessorDict."""
+
+        # Generate the computer code and export it on disk for the remaining new currents
+        current_exporter = subtraction.SubtractionCurrentExporter(model, root_path, current_set)
+        mapped_currents = current_exporter.export(currents_to_consider)
+        # Print to the debug log which currents were exported
+        log_string = "The following subtraction current implementation are exported "+\
+                     "in contribution '%s':\n"%self.short_name()
+        for (module_path, class_name, _), current_properties in mapped_currents.items():
+            if class_name != 'DefaultCurrentImplementation':
+                quote_class_name = "'%s'" % class_name
+                defining_current_str = str(current_properties['defining_current'])
+                line_pars = (quote_class_name, defining_current_str)
+                log_string += " > %-35s for representative current '%s'\n" % line_pars
+            else:
+                quote_default_name = "'DefaultCurrentImplementation'"
+                number_of_currents = len(current_properties['mapped_process_keys'])
+                line_pars = (quote_default_name, number_of_currents)
+                log_string += " > %-35s for a total of %d currents.\n" % line_pars
+        logger.debug(log_string)
+        # Instantiate the CurrentAccessors corresponding
+        # to all current implementations identified and needed
+        all_current_accessors = []
+        for (module_path, class_name, _), current_properties in mapped_currents.items():
+            all_current_accessors.append(accessors.VirtualMEAccessor(
+                current_properties['defining_current'],
+                module_path,
+                class_name,
+                '%s.subtraction_current_implementations_utils'%current_exporter.main_module_name,
+                current_properties['instantiation_options'], 
+                mapped_process_keys=current_properties['mapped_process_keys'],
+                root_path=root_path,
+                model=model
+            ))
+        # Add MEAccessors
+        all_MEAccessors.add_MEAccessors(all_current_accessors)
+        return mapped_currents
+
     def can_processes_be_integrated_together(self, processA, processB):
         """ Investigates whether processA and processB can be integrated together for this contribution."""
 
@@ -636,10 +834,6 @@
             all_integrands.extend(self.get_integrands_for_process_map(integrand_process_map, *args))
         
         return all_integrands
-        
-    def add_content_to_global_ME7_resources(self, global_ME7_dir, **opts):
-        """ Possibly add content to the global ME7 resources pool."""
-        pass
 
     def remove_superfluous_content(self):
         """At the end of the export of this contributions,
@@ -669,6 +863,11 @@
         # This stores the wanted couplings which should be exported by the overall ME7 exporter.
         global_wanted_couplings = []
         
+        # For some contributions, like the beam factorization ones, there is no additional
+        # export to be performed.
+        if self.has_beam_factorization():
+            return global_wanted_couplings
+
         # Handling of the model
         if self.options['_model_v4_path']:
             logger.info(
@@ -700,6 +899,11 @@
 
     def link_global_ME7_resources(self, global_ME7_dir):
         """ Remove some of the local resources and instead link the global ones from the ME7 exporters."""
+
+        # For some contributions, like the beam factorization ones, there is no additional
+        # export to be performed.
+        if self.has_beam_factorization(): 
+            return
         
         # Link to the global model if necessary
         if not self.requires_its_own_model:
@@ -714,6 +918,12 @@
 
     def make_model_symbolic_link(self):
         """ Links some of the MODEL resources to the other folders exported within this contribution."""
+        
+        # For some contributions, like the beam factorization ones, there is no additional
+        # export to be performed.
+        if self.has_beam_factorization():
+            return        
+        
         self.exporter.make_model_symbolic_link()
     
     def get_inverse_processes_map(self, force=False):
@@ -725,11 +935,7 @@
         
         # Sort PDGs only when processes have been grouped
         sort_PDGs = self.group_subprocesses
-        
-        if not self.all_matrix_elements.get_matrix_elements():
-            raise MadGraph5Error("The processes reverse map can only be built *after*"+
-                                                 " the matrix element have been exported.")
-        
+
         if hasattr(self, 'inverse_processes_map') and not force:
             return self.inverse_processes_map
             
@@ -761,9 +967,7 @@
             where defining_process is an actual instance of Process and [mapped processes] is a list 
             of processes instances mapped to that defining process.
         """
-        if not self.amplitudes:
-            return {}
-        
+
         # Attempt reusing a cached version
         if hasattr(self, 'processes_map') and not force:
             if (self.processes_map[0]['had_amplitudes'] == bool(self.amplitudes)) and \
@@ -771,6 +975,9 @@
                                           self.all_matrix_elements.get_matrix_elements())):
                 return self.processes_map[1]
         
+        if not self.amplitudes:
+            return {}
+        
         all_defining_procs = [amp.get('process') for amp in self.amplitudes]
         
         # The values are (defining_process_instance, list_of_mapped_process_instances)
@@ -803,7 +1010,11 @@
                 all_defining_procs[proc][1].extend([all_procs_pdgs[p] for p in all_procs_pdgs
                                             if p!=proc and p not in all_defining_procs[proc][1]])
 
-
+        # Insure that all processes have leg numbers that are consecutive
+        for process_key, (defining_process, mapped_processes) in all_defining_procs.items():
+            for p in [defining_process,]+mapped_processes:
+                for i, leg in enumerate(p.get('legs')):
+                    leg.set('number', i+1)
 
         # Cache the process map
         self.processes_map = ({
@@ -815,12 +1026,21 @@
 
     def process_dir_name(self, process):
         """ Given a specified process, return the directory name in which it is output."""
-        # This function is essentially useful in order to be shadowed by daughter classes     
-        return process.shell_string()       
-    
+        # This function is essentially useful in order to be shadowed by daughter classes
+        if self.output_type == 'madloop':
+            return "%d_%s_%s"%(process.get('id'), process.get('uid'),
+                                                      process.shell_string(print_id=False))
+        else:
+            return process.shell_string()
+
     def compile(self):
         """ Compiles the f2py shared library to provide easy access to this contribution."""
         
+        # This is not needed for factorization contributions as they always use processes already
+        # in other contributions.
+        if self.has_beam_factorization():
+            return
+        
         processes_map = self.get_processes_map()
         i_process=0
         for process_key, (defining_process, mapped_processes) in processes_map.items():
@@ -837,15 +1057,16 @@
             if not os.path.isfile(pjoin(Pdir, 'matrix_%s_py.so'%name)):
                 raise InvalidCmd("The f2py compilation of subprocess '%s' failed.\n"%Pdir+
                                  "Try running 'make matrix2py.so' by hand in this directory.")
-            ln(pjoin(Pdir, 'matrix_%s_py.so'%name), starting_dir=self.export_dir)
-    
+            ln(pjoin(Pdir, 'matrix_%s_py.so'%name), starting_dir=self.export_dir)           
+
     def get_nice_string_process_line(self, process_key, defining_process, format=0):
         """ Return a nicely formated process line for the function nice_string of this 
         contribution."""
         GREEN = '\033[92m'
-        ENDC = '\033[0m'        
-        return GREEN+'  %s'%defining_process.nice_string(print_weighted=False).\
-                                                               replace('Process: ','')+ENDC
+        ENDC = '\033[0m'
+
+        return GREEN+'  %s'%(
+            defining_process.nice_string(print_weighted=False).replace('Process: ',''))+ENDC
     
     def get_additional_nice_string_printout_lines(self, format=0):
         """ Return additional information lines for the function nice_string of this contribution."""
@@ -860,7 +1081,8 @@
         BLUE = '\033[94m'
         GREEN = '\033[92m'
         ENDC = '\033[0m'
-        res = ['%-30s:   %s'%('contribution_type',type(self))]
+        res = ['< %s%s%s >'%(BLUE,self.short_name(),ENDC)]
+        res.append('%-30s:   %s'%('contribution_type',type(self)))
         res.extend([self.contribution_definition.nice_string()])
         if not self.topologies_to_processes is None:
             res.append('%-30s:   %d'%('Number of topologies', 
@@ -875,7 +1097,8 @@
                 for amp in self.amplitudes:
                     res.append(GREEN+'  %s'%amp.get('process').nice_string(print_weighted=False).\
                                                                             replace('Process: ','')+ENDC)
-        elif self.amplitudes and self.all_matrix_elements.get_matrix_elements():
+        elif self.amplitudes and self.all_matrix_elements.get_matrix_elements() or \
+                                                             self.has_beam_factorization():
             processes_map = self.get_processes_map()
             if format < 1:
                 res.append('Generated and mapped processes for this contribution: %d (+%d mapped)'%
@@ -892,9 +1115,157 @@
             res.append(BLUE+'No amplitudes generated yet.'+ENDC)
             
         return '\n'.join(res).encode('utf-8')
+    
+    def has_beam_factorization(self):
+        """ Checks whether this contribution has as any beam factorization active."""
+
+        return self.contribution_definition.has_beam_factorization()
+    
+    def import_topologies_from_contrib(self, process_key, contrib):
+        """ Add to the current topologies of this contribution, the ones defined for the
+        process key of the contribution passed in argument."""
+        
+        # TODO a proper combination needs to be performed here. For now, I simply add them
+        # to the attribute processes_to_topologies and topologies_to_processes, but that
+        # doesn't work because the topology_key used is based on the 
+        #    (subprocess_group_number, config_number)
+        # which will clash between the various contributions added. So matching configurations
+        # must be combined and the different ones must be added with a topology key upgraded
+        # so as to make them unique (for example using a combination of:
+        #    contrib.contribution_definiton.process_definition.get('id')
+        # and
+        #    contrib.contribution_definiton.n_loops
+        #.which is best done when generating the topologies in the first place anyway.
+        # But they will need to be combined anyway.
+        #
+        # For now we implement a minimal behaviour below which is supposed to work for
+        # NLO with a simple 'add process' command, where we can basically copy the topologies
+        # directly from the Born
+        self.processes_to_topologies[process_key] = contrib.processes_to_topologies[process_key]
+        
+        for (topology_key, processes) in contrib.topologies_to_processes.items():
+                self.topologies_to_processes[topology_key] = processes
+
+    def add_beam_factorization_processes_from_contribution(self, contrib, beam_factorization_order):
+        """ 'Import' the processes_map and all generated attributes of the contribution
+        passed in argument to assign them to this beam factorization contributions, and
+        generate all necessary "bulk" BeamCurrents as we go along.
+        """
+        
+        beam_factorization_configuration = (
+           (not self.contribution_definition.beam_factorization['beam_one'] is None) and \
+            self.contribution_definition.beam_factorization['beam_one']['active'],
+           (not self.contribution_definition.beam_factorization['beam_two'] is None) and \
+            self.contribution_definition.beam_factorization['beam_two']['active'] )
+
+        # If both beams must be assigned a factorisation term, then make sure that the
+        # beam_factorization_order is at least two
+        if beam_factorization_configuration == (True,True) and beam_factorization_order < 2:
+            # We have nothing to add from this contribution.
+            return 
+        
+        # Make sure the source contribution does have a factorizing beam for the active
+        # ones of this contribution.
+        if ( (beam_factorization_configuration[0] and 
+              (contrib.contribution_definition.beam_factorization['beam_one'] is None)) or
+             (beam_factorization_configuration[1] and 
+               (contrib.contribution_definition.beam_factorization['beam_two'] is None)) ):
+            return
+        
+        # Now fetch the process map to load from
+        source_processes_map = contrib.get_processes_map()
+        processes_map = self.processes_map[1]
+        # To the new contribution instances, also assign references to the generated attributes
+        # of the original contributions (no deep copy necessary here)
+        for process_key, (defining_process, mapped_processes) in source_processes_map.items():
+            
+            # First add this process to the current processes map
+            if process_key in processes_map:
+                raise MadGraph5Error('MadNkLO attempts to add beam factorization terms for'+
+                    ' the %s a second time. This should never happen.'%(
+                             defining_process.nice_string().replace('Process','process')))
+            processes_map[process_key] = (
+                  defining_process.get_copy(), [mp.get_copy() for mp in mapped_processes] )
+            
+            # Create the corresponding BeamCurrents
+            all_beam_currents = []
+                        
+            beam_one_type = contrib.contribution_definition.beam_factorization['beam_one']['beam_type']
+            beam_two_type = contrib.contribution_definition.beam_factorization['beam_two']['beam_type']
+            beam_one_PDGs = contrib.contribution_definition.beam_factorization['beam_one']['beam_PDGs']
+            beam_two_PDGs = contrib.contribution_definition.beam_factorization['beam_two']['beam_PDGs']
+            # The double for-loop below creates all compatible assignments of n_loop to the beam
+            # factorisation currents applied to the first and second beam. So (2,1),(1,2) at N^3LO for instance
+            # In principle it would also be possible here to further differentiate P^(0) x P^(0) from P^(1), but
+            # we do not do this for now.
+            for factorization_order_beam_one in range(beam_factorization_order+1):
+                if factorization_order_beam_one > 0 and not beam_factorization_configuration[0]:
+                    break
+                for factorization_order_beam_two in range(beam_factorization_order-factorization_order_beam_one+1):
+                    if factorization_order_beam_one == factorization_order_beam_two == 0:
+                        continue
+                    if factorization_order_beam_two > 0 and not beam_factorization_configuration[1]:
+                        break
+                    initial_legs = { l.get('number'): l for l in defining_process.get_initial_legs() }
+                    all_beam_currents.append({
+                        'beam_one' : None if factorization_order_beam_one == 0 else
+                          subtraction.BeamCurrent({
+                            'beam_type'         : beam_one_type,
+                            'beam_PDGs'         : beam_one_PDGs,
+                            'distribution_type' : 'bulk',
+                            'n_loops'           : factorization_order_beam_one-1,
+                            'perturbation_couplings' : self.contribution_definition.correction_couplings,
+                            'squared_orders'    : { order : factorization_order_beam_one*2 for
+                                   order in self.contribution_definition.correction_couplings },
+                            'sqorders_types'    : { order : '<=' for order in 
+                                            self.contribution_definition.correction_couplings },
+                            'singular_structure': subtraction.SingularStructure(
+                                 substructures = [ subtraction.BeamStructure(
+                                                                 legs = [initial_legs[1],]), ]
+                            ),
+                          }),
+                        'beam_two' : None if factorization_order_beam_two == 0 else
+                          subtraction.BeamCurrent({
+                            'beam_type'         : beam_two_type,
+                            'beam_PDGs'         : beam_two_PDGs,
+                            'distribution_type' : 'bulk',
+                            'n_loops'           : factorization_order_beam_two-1,
+                            'perturbation_couplings' : self.contribution_definition.correction_couplings,
+                            'squared_orders'    : { order : factorization_order_beam_two*2 for
+                                   order in self.contribution_definition.correction_couplings },
+                            'sqorders_types'    : { order : '<=' for order in 
+                                            self.contribution_definition.correction_couplings },
+                            'singular_structure': subtraction.SingularStructure(
+                                 substructures = [ subtraction.BeamStructure(
+                                                                 legs = [initial_legs[2],]), ]
+                            ),
+                          })
+                    })
+
+            # Add these beam currents to the processes just added to this contribution
+            processes_map[process_key][0].set('beam_factorization', all_beam_currents)
+            for mp in processes_map[process_key][1]:
+                mp.set('beam_factorization', all_beam_currents)
+
+            # Also add the corresponding topologies
+            self.import_topologies_from_contrib(process_key, contrib)
+
+    def generate_all_counterterms(self, 
+                        ignore_integrated_counterterms=False, group_processes=True, **opts):
+        """ Generate all counterterms associated to the processes in this contribution.
+        Since we are in the mother class Contribution here, we only consider generating the
+        beam factorisation counterterm of type 'local_CT'."""
+        
+        # Initialise local and integrated counterterms to an empty list
+        self.counterterms = {}
+        all_integrated_counterterms = []
+        
+        # Nothing done by default
+        return all_integrated_counterterms
         
     def generate_amplitudes(self, force=False):
-        """ Generates the relevant amplitudes for this contribution."""
+        """ Generates the relevant amplitudes for this contribution and the construction
+        of the instances of currents for the beam factorisation terms."""
         
         # First check if the amplitude was not already generated
         if self.amplitudes and not force:
@@ -903,7 +1274,8 @@
         myproc = self.MultiProcessClass(self.contribution_definition.process_definition,
                     collect_mirror_procs = self.collect_mirror_procs,
                     ignore_six_quark_processes = self.ignore_six_quark_processes,
-                    optimize=self.optimize, diagram_filter=self.diagram_filter)
+                    optimize=self.optimize, diagram_filter=self.diagram_filter,
+                    loop_filter=self.loop_filter)
     
         for amp in myproc.get('amplitudes'):
             if amp not in self.amplitudes:
@@ -919,16 +1291,10 @@
         super(Contribution_B,self).__init__(contribution_definition, cmd_interface, **opts)
         self.ignore_six_quark_processes = self.options['ignore_six_quark_processes'] if \
             "ignore_six_quark_processes" in self.options else []
-        if 'diagram_filter' in opts:
-            self.diagram_filter = opts['diagram_filter']
         if 'optimize' in opts:
             self.optimize = opts['optimize']
         self.supports_decay_chain       = True
 
-class Contribution_LIB(Contribution_B):
-    """ Implements the handling of loop-induced Born-type of contribution."""
-    pass
-
 class Contribution_R(Contribution):
     """ Implements the handling of a single real-emission type of contribution."""
 
@@ -953,8 +1319,8 @@
         
         GREEN = '\033[92m'
         ENDC = '\033[0m'        
-        res =  GREEN+'  %s'%defining_process.nice_string(print_weighted=False).\
-                                                               replace('Process: ','')+ENDC
+        res =  GREEN+'  %s'%(defining_process.nice_string(print_weighted=False).\
+                                                              replace('Process: ',''))+ENDC
 
         if not self.counterterms:
             return res                                                               
@@ -983,10 +1349,11 @@
 
     def generate_all_counterterms(self, ignore_integrated_counterterms=False, group_processes=True):
         """ Generate all counterterms associated to the processes in this contribution."""
-        
-        self.counterterms = {}
+
         all_integrated_counterterms = []
+
         for process_key, (defining_process, mapped_processes) in self.get_processes_map().items():
+
             local_counterterms, integrated_counterterms =  self.IR_subtraction.get_all_counterterms(
                     defining_process, ignore_integrated_counterterms=ignore_integrated_counterterms)
             
@@ -1001,13 +1368,19 @@
             # matrix elements themselves are not added here
             local_counterterms = [ct for ct in local_counterterms if ct.is_singular()]
 
+            # Set the reduced flavor map of all the counterterms
+            for CT in integrated_counterterms:
+                CT.set_reduced_flavors_map(defining_process, mapped_processes, self.IR_subtraction)
+            for CT in local_counterterms:
+                CT.set_reduced_flavors_map(defining_process, mapped_processes, self.IR_subtraction)
+
 #            misc.sprint('Local CTs for %s'%(defining_process.nice_string()))
 #            for lc in local_counterterms:
 #                misc.sprint(str(lc))
 #            misc.sprint('Integrated CTs for %s'%(defining_process.nice_string()))
 #            for ic in integrated_counterterms:
 #                misc.sprint(str(ic))
-                
+
             self.counterterms[process_key] = local_counterterms
 
             # Now Add a additional information about the integrated_counterterms that will
@@ -1026,20 +1399,22 @@
                 initial_pdgs = process.get_initial_ids()
                 final_pdgs = tuple(process.get_final_ids_after_decay())
                 flavors_combinations.append( ( ( tuple(initial_pdgs),final_pdgs ), leg_numbers ) )
-#               It is not completely clear how to account for mirrored process at this stage,
-#               so we anticipate for now that this will only need to be processed at the
-#               integrand evaluation level.
-#               The problem with the line below is that for the process u d~ > u d~ a for instance
-#               the counterterm C(1,3) is only valid if leg #1 is u and leg #2 is d~. 
+#               Mirrored processes are not accounted for at this stage, as these are processed at the
+#               integrand evaluation level where mirrored copy of events are generated.
+#               Notice that this mirroring shouldn't be done simply as in the line below.
+##
+##                if process.get('has_mirror_process'):
+##                    flavors_combinations.append( 
+##                           ( ( tuple(reversed(initial_pdgs)), final_pdgs ), leg_numbers ) )
+##
+#               because if one takes the following process for instance:
+#                     u d~ > u d~ a
+#               then the counterterm C(1,3) is only valid if leg #1 is u and leg #2 is d~. 
 #               So swapping u(1) d~(2) for u(2) d~(1) would require propagating the indices
 #               substitution in the entire counterterm (most importantly, its singular structure)
 #               and create a copy of it. This pretty heavy (and the idea of the 'use_mirror_process'
-#               option is precisely to avoid this work), so we can hopefully bypass this
+#               option is precisely to avoid this work), so this is entirely bypassed
 #               by implementing this symmetrisation during the integrand evaluation.
-##                if process.get('has_mirror_process'):
-##                    flavors_combinations.append( 
-##                           ( ( tuple(reversed(initial_pdgs)), final_pdgs ), leg_numbers ) )
-
             
 #            misc.sprint('doing process:',defining_process.nice_string())
 #            misc.sprint('flavors_combinations',flavors_combinations)
@@ -1072,14 +1447,13 @@
                 except KeyError:
                     # Register the new integrated current topology
                     integrated_current_topologies[integrated_current_topology] = [integrated_counterterm,]
-            
-            
+
             # For each topology, count the number of counterterms that are associated to
             # various list of external leg numbers. Within a given topology, one expects
             # that there is the same number of counterterm repetition for each combination of external legs
             # Example: resolved process is
             #   e+(1) e-(2) > g(3) g(4) g(5) d(6) d~(7)
-            # All the following counterterm belong to the same topology, which we group here according
+            # All the following counterterms belong to the same topology, which we group here according
             # to which list of external leg numbers they include:
             #   (3,4,6) : C(C(S(3),4),6), C(C(S(4),3),6) 
             #   (3,5,6) : C(C(S(3),5),6), C(C(S(5),3),6) 
@@ -1125,6 +1499,32 @@
                 # flavors combination and the value is the multiplication factor.
                 reduced_flavors_combinations    = {}
                 resolved_flavors_combinations   = {}
+                # One subtlety occurs because of the backward flavor evolution in the initial
+                # states which will be registered in the event since it must be convoluted with
+                # the correct PDFs.
+                # Consider the following process:
+                #             g(1) s(2) > s(3) c(4) c~(5) u(6) u~(7)
+                #   + mapped  g(1) q(2) > q(3) qprime(4) qprime~(5) u(6) u~(7)
+                # and the counterterm C(2,3)C(4,5) which yields the following reduced mapped flavor:
+                #             g g > g u u~             
+                # with a multiplicity of naively n_f**2. However, the actual events that will be generated
+                # will be the following:
+                #             g s > g u u~
+                #             g d > g u u~
+                #             ...
+                # It is then clear that each of this flavor configuration in the ME7 event should not be
+                # multiplied by n_f**2 but instead just n_f. For this reason, the function 
+                # generate_all_counterterms of the contribution class also produces the dictionary 
+                # 'reduced_flavors_with_resolved_initial_states_combinations' which, in the example above, 
+                # would be (nf=5):
+                #     {  (21,3),(21,2,-2) : 5,
+                #        (21,1),(21,2,-2) : 5,
+                #        ...
+                #     }
+                # which will allow to divide the weight of each of these flavor configurations at the very
+                # by the corresponding multiplicity factor n_f.
+                reduced_flavors_with_resolved_initial_states_combinations = {}
+                
                 # Each integrated counterterm contributes proportionally to the ratio
                 # of the symmetry factors of the resolved process and reduced one,
                 # further divided by the symmetry factor of each group external leg
@@ -1147,14 +1547,15 @@
                 # There should always be a unique value of S_t for each reduced_flavor
                 symmetry_factors                = {}
                 for flavor_combination, leg_numbers in flavors_combinations:
+                    # Access the reduced flavors
+                    reduced_flavors = integrated_counterterm.get_reduced_flavors(
+                                                    defining_flavors = flavor_combination )
                     # Create a map from leg number to flavor
                     number_to_flavor_map = dict( 
                         [ ( number, flavor_combination[0][i] ) for i, number in enumerate(leg_numbers[0]) ] +
                         [ ( number, flavor_combination[1][i] ) for i, number in enumerate(leg_numbers[1]) ]
                     )
-                    reduced_flavors = integrated_counterterm.get_reduced_flavors(
-                      defining_flavors=number_to_flavor_map, IR_subtraction=self.IR_subtraction)
-#                    misc.sprint('From resolved flavor:',flavor_combination)
+#                    misc.sprint('From resolved flavor:',number_to_flavor_map)
 #                    misc.sprint('I get:',reduced_flavors)
 
                     resolved_symmetry_factor = misc.symmetry_factor(list(flavor_combination[1]))
@@ -1169,15 +1570,15 @@
                             l in singular_legs if l.state==subtraction.SubtractionLeg.FINAL]
                         singular_legs_symmetry_factor *= misc.symmetry_factor(singular_external_leg_pdgs)
 
-   #                 misc.sprint('Counterterm: %s'%str(integrated_counterterm))
-   #                 misc.sprint('    -> multiplicity                   = %d'%len(integrated_counterterms))
-   #                 misc.sprint('    -> CT with same topologies        = %s'%(' | '.join(str(CT) for CT in integrated_counterterms)) )                    
-   #                 misc.sprint('    -> resolved_flavors               = %s'%str(flavor_combination))
-   #                 misc.sprint('    -> reduced                        = %s'%str(reduced_flavors))
-   #                 misc.sprint('    -> resolved_symmetry_factor       = %d'%resolved_symmetry_factor)
-   #                 misc.sprint('    -> reduced_symmetry_factor        = %d'%reduced_symmetry_factor)
-   #                 misc.sprint('    -> singular_legs_symmetry_factor  = %f'%singular_legs_symmetry_factor)
-   #                 misc.sprint('    -> counterterm_repetition_factor  = %f'%repetitions_for_topology[integrated_current_topology])
+#                    misc.sprint('Counterterm: %s'%str(integrated_counterterm))
+#                    misc.sprint('    -> multiplicity                   = %d'%len(integrated_counterterms))
+#                    misc.sprint('    -> CT with same topologies        = %s'%(' | '.join(str(CT) for CT in integrated_counterterms)) )                    
+#                    misc.sprint('    -> resolved_flavors               = %s'%str(flavor_combination))
+#                    misc.sprint('    -> reduced                        = %s'%str(reduced_flavors))
+#                    misc.sprint('    -> resolved_symmetry_factor       = %d'%resolved_symmetry_factor)
+#                    misc.sprint('    -> reduced_symmetry_factor        = %d'%reduced_symmetry_factor)
+#                    misc.sprint('    -> singular_legs_symmetry_factor  = %f'%singular_legs_symmetry_factor)
+#                    misc.sprint('    -> counterterm_repetition_factor  = %f'%repetitions_for_topology[integrated_current_topology])
                         
                     # Now correct for the repetition number in that topology
                     # Notice it should always be an integer, the float here is just so that
@@ -1192,6 +1593,16 @@
                         reduced_flavors_combinations[reduced_flavors] = 1
                     resolved_flavors_combinations[flavor_combination] = reduced_flavors
                     
+                    # Now also build a multiplicity factor for combination of resolved initial-state
+                    # flavors and reduced final-state ones.
+                    reduced_flavors_with_resolved_initial_states = (flavor_combination[0], reduced_flavors[1])
+                    try:
+                        reduced_flavors_with_resolved_initial_states_combinations[
+                                            reduced_flavors_with_resolved_initial_states] += 1
+                    except KeyError:
+                        reduced_flavors_with_resolved_initial_states_combinations[
+                                            reduced_flavors_with_resolved_initial_states] = 1
+                    
                     if reduced_flavors in symmetry_factors:
                         # Sanity check, all reduced flavors should share the same S_t
                         if overall_symmetry_factor != symmetry_factors[reduced_flavors]:
@@ -1217,6 +1628,8 @@
                     'integrated_counterterm'             :   integrated_counterterm,
                     'resolved_flavors_combinations'      :   resolved_flavors_combinations,
                     'reduced_flavors_combinations'       :   reduced_flavors_combinations,
+                    'reduced_flavors_with_resolved_initial_states_combinations' :
+                                        reduced_flavors_with_resolved_initial_states_combinations,
                     'symmetry_factor'                    :   symmetry_factors.values()[0],
                     'multiplicity'                       :   len(integrated_counterterms)
                 })
@@ -1228,25 +1641,46 @@
         """Given the list of available reduced processes encoded in the MEAccessorDict 'all_MEAccessors'
         given in argument, remove all the counterterms whose underlying reduced process does not exist."""
 
-        for counterterm in list(counterterms):
+        all_removed = True
+        for counterterm_info in list(counterterms):
+            # Integrated counterterm come as dictionaries with additional metadata
+            if isinstance(counterterm_info, dict):
+                counterterm = counterterm_info['integrated_counterterm']
+            else:
+                counterterm = counterterm_info
             # Of course don't remove any counterterm that would be the real-emission itself.
             if not counterterm.is_singular():
                 continue
             try:
                 all_MEAccessors.get_MEAccessor(counterterm.process)
+                all_removed = False
             except MadGraph5Error:
                 # This means that the reduced process could not be found and
                 # consequently, the corresponding counterterm must be removed.
                 # Example: C(5,6) in e+ e- > g g d d~
-                counterterms.remove(counterterm)
-
-    def export(self, ignore_integrated_counterterms=False, **opts):
-        """ Overloads export so as to export subtraction currents as well."""
-        ret_value = super(Contribution_R, self).export(**opts)
+                # Useful debug lines below
+#                misc.sprint('-'*50)
+#                import pprint
+#                misc.sprint('Process key:\n%s'%(
+#                    pprint.pformat(dict(ProcessKey(counterterm.process).get_canonical_key()))))
+#                misc.sprint('-'*50)
+#                misc.sprint('All process keys in the accessor:\n')
+#                for key in all_MEAccessors:
+#                    misc.sprint(pprint.pformat(key))
+#                misc.sprint('-'*50)
+                counterterms.remove(counterterm_info)
+        
+        return all_removed
+
+    def subtract(self, ignore_integrated_counterterms=False, **opts):
+        """ Generate and export all necessary subtraction counterterterms and currents
+        (including beam factorization integrated counterterms [bulk ones are already generated
+        at this stage])."""
+
+        ret_value = super(Contribution_R,self).subtract(ignore_integrated_counterterms=False, **opts)
 
         # Fish out the group_processes option as it could be used when attempting to
         # generate all currents.
-        
         if 'group_processes' in opts:
             group_processes = opts['group_processes']
         else:
@@ -1258,7 +1692,7 @@
       
         # Add the integrated counterterms to be passed to the exporter
         ret_value.update({'integrated_counterterms': integrated_counterterms})
-
+        
         return ret_value
 
     def get_all_necessary_local_currents(self, all_MEAccessors):
@@ -1276,55 +1710,50 @@
                 if copied_current not in all_currents:
                     all_currents.append(copied_current)
 
+        # Also add the beam_factorization currents from process itself if it has any
+        for process_key, (defining_process, mapped_processes) in self.get_processes_map().items():
+            for beam_factorization_currents in defining_process['beam_factorization']:
+                for bfc in beam_factorization_currents.values():
+                    if bfc is not None:
+                        all_currents.append(bfc)
+
         # Now further remove currents that are already in all_MEAccessors
-        all_currents = [current
-                        for current in all_currents
-                        if current.get_key().get_canonical_key() not in all_MEAccessors]
-        return all_currents
-
-    @classmethod
-    def add_current_accessors(
-        cls, model, all_MEAccessors, root_path, current_set, currents_to_consider):
-        """Generate and add all subtraction current accessors to the MEAccessorDict."""
-
-        # Generate the computer code and export it on disk for the remaining new currents
-        current_exporter = subtraction.SubtractionCurrentExporter(
-            model, root_path, current_set)
-        mapped_currents = current_exporter.export(currents_to_consider)
-        # Print to the debug log which currents were exported
-        log_string = "The following subtraction current implementation are exported:\n"
-        for (module_path, class_name, _), current_properties in mapped_currents.items():
-            if class_name != 'DefaultCurrentImplementation':
-                quote_class_name = "'%s'" % class_name
-                defining_current_str = str(current_properties['defining_current'])
-                line_pars = (quote_class_name, defining_current_str)
-                log_string += " > %-35s for representative current '%s'\n" % line_pars
-            else:
-                quote_default_name = "'DefaultCurrentImplementation'"
-                number_of_currents = len(current_properties['mapped_process_keys'])
-                line_pars = (quote_default_name, number_of_currents)
-                log_string += " > %-35s for a total of %d currents." % line_pars
-        logger.debug(log_string)
-        # Instantiate the CurrentAccessors corresponding
-        # to all current implementations identified and needed
-        all_current_accessors = []
-        for (module_path, class_name, _), current_properties in mapped_currents.items():
-            all_current_accessors.append(accessors.VirtualMEAccessor(
-                current_properties['defining_current'],
-                module_path,
-                class_name,
-                '%s.subtraction_current_implementations_utils'%current_exporter.main_module_name,
-                current_properties['instantiation_options'], 
-                mapped_process_keys=current_properties['mapped_process_keys'],
-                root_path=root_path,
-                model=model
-            ))
-        # Add MEAccessors
-        all_MEAccessors.add_MEAccessors(all_current_accessors)
-        return mapped_currents
-
-    def remove_zero_counterterms(self, all_ME_accessors):
-
+        all_necessary_currents = {}
+        for current in all_currents:
+            current_key = current.get_key().get_canonical_key()
+            if current_key not in all_MEAccessors and current_key not in all_necessary_currents:
+                all_necessary_currents[current_key] = current
+
+        return all_necessary_currents.values()
+
+    def add_ME_accessors(self, all_MEAccessors, root_path):
+        """ Adds all MEAccessors for the matrix elements and currents generated as part of this contribution."""
+        
+        # Get the basic accessors for the matrix elements
+        super(Contribution_R, self).add_ME_accessors(all_MEAccessors, root_path)
+
+        # Only initialize all_counterterms_removed to True if there is CT to removed to begin with
+        all_counterterms_removed = any(len(CTs)>0 for CTs in self.counterterms.values())    
+        for process_key, counterterms in self.counterterms.items():
+            # Remove counterterms with non-existing underlying Born processes
+            if not self.remove_counterterms_with_no_reduced_process(all_MEAccessors, counterterms):
+                all_counterterms_removed = False
+        # Sanity check
+        if all_counterterms_removed:
+            logger.warning("All local counterterms from contribution '%s' were removed because"%self.short_name()+
+                " their corresponding reduced processes was not found in the list of processes exported.\n"+
+                " This is likely wrong and originates from a mismatch in the process definitions of the various"+
+                " contributions building this higher order computation.")
+
+        # Obtain all necessary currents
+        current_set = self.options['subtraction_currents_scheme']
+        currents_to_consider = self.get_all_necessary_local_currents(all_MEAccessors)
+        self.add_current_accessors(
+            self.model, all_MEAccessors, root_path, current_set, currents_to_consider )
+        
+    def remove_local_counterterms_set_to_zero(self, all_ME_accessors):
+        """ Remove all local counterterms involving currents whose implementation has set
+        them to zero."""
         for process_key, counterterms in self.counterterms.items():
             for counterterm in list(counterterms):
                 # misc.sprint("Considering CT %s" % str(counterterm))
@@ -1336,32 +1765,17 @@
                             counterterms.remove(counterterm)
                             break
 
-    def add_ME_accessors(self, all_MEAccessors, root_path):
-        """ Adds all MEAccessors for the matrix elements and currents generated as part of this contribution."""
-        
-        # Get the basic accessors for the matrix elements
-        super(Contribution_R, self).add_ME_accessors(all_MEAccessors, root_path)
-
-        for process_key, counterterms in self.counterterms.items():
-            # Remove counterterms with non-existing underlying Born processes
-            self.remove_counterterms_with_no_reduced_process(all_MEAccessors, counterterms)
-
-        # Obtain all necessary currents
-        current_set = self.options['subtraction_currents_scheme']
-        currents_to_consider = self.get_all_necessary_local_currents(all_MEAccessors)
-        self.add_current_accessors(
-            self.model, all_MEAccessors, root_path, current_set, currents_to_consider )
-        self.remove_zero_counterterms(all_MEAccessors)
-     
     def get_integrands_for_process_map(self, process_map, model, run_card, all_MEAccessors, ME7_configuration):
         """ Returns all the integrands implementing this contribution for the specified process_map.
         The instance of MEAccessorDict is necessary so as to be passed to the integrand instances.
         """
         
         relevant_counterterms = {}
-        self.remove_zero_counterterms(all_MEAccessors)
-        for process_key in process_map:
-            relevant_counterterms[process_key] = self.counterterms[process_key]
+        self.remove_local_counterterms_set_to_zero(all_MEAccessors)
+
+        if self.counterterms:
+            for process_key in process_map:
+                relevant_counterterms[process_key] = self.counterterms[process_key]
 
         return [
             ME7_integrands.ME7Integrand(
@@ -1387,10 +1801,7 @@
     def __init__(self, contribution_definition, cmd_interface, **opts):
         """ Bring in the couple of modifications necessary for this type of contributions."""
         super(Contribution_V,self).__init__(contribution_definition, cmd_interface, **opts)
-        # Make sure to adjust the MultiProcessClass to be used
-        self.MultiProcessClass          = loop_diagram_generation.LoopMultiProcess
-        self.output_type                = 'madloop'
-        
+
         # Store values being integration counterterms and their properties (as a dictionary), 
         # with the ProcessKey they are attached to as keys of this attribute's dictionary.
         self.integrated_counterterms    = {}
@@ -1408,8 +1819,8 @@
         contribution."""
         GREEN = '\033[92m'
         ENDC = '\033[0m'        
-        res = GREEN+'  %s'%defining_process.nice_string(print_weighted=False).\
-                                                               replace('Process: ','')+ENDC
+        res = GREEN+'  %s'%(defining_process.nice_string(print_weighted=False)
+                                                             .replace('Process: ',''))+ENDC
 
         if not self.integrated_counterterms:
             return res
@@ -1440,66 +1851,6 @@
             res += '\n'.join(long_res)
 
         return res
-    
-    @classmethod
-    def remove_counterterms_with_no_reduced_process(cls, all_MEAccessors, CT_properties_list):
-        """ Given the list of available reduced processes encoded in the MEAccessorDict 'all_MEAccessors'
-        given in argument, remove all the counterterms whose underlying reduced process does not exist."""
-
-        for CT_properties in list(CT_properties_list):
-            counterterm = CT_properties['integrated_counterterm']
-            # Of course don't remove any counterterm that would be the real-emission itself.
-            if not counterterm.is_singular():
-                continue
-            try:
-                all_MEAccessors.get_MEAccessor(counterterm.process)
-            except MadGraph5Error:
-                # This means that the reduced process could not be found and
-                # consequently, the corresponding counterterm must be removed.
-                # Example: C(5,6) in e+ e- > g g d d~
-                CT_properties_list.remove(CT_properties)
-
-    def generate_matrix_elements(self, group_processes=True):
-        """Generate the Helas matrix elements before exporting. Uses the main function argument 
-        'group_processes' to decide whether to use group_subprocess or not."""
-
-        cpu_time1 = time.time()
-        ndiags = 0
-
-        generation_mode = {'optimized_output': self.options['loop_optimized_output']}
-
-        self.all_matrix_elements = loop_helas_objects.LoopHelasProcess(
-            self.amplitudes,
-            compute_loop_nc = True,
-            matrix_element_opts = generation_mode,
-            optimized_output = self.options['loop_optimized_output'],
-            combine_matrix_elements=group_processes
-        )
-        ndiags = sum([len(me.get('diagrams')) for me in self.all_matrix_elements.get_matrix_elements()])
-        
-        # assign a unique id number to all process
-        uid = 0 
-        for me in self.all_matrix_elements.get_matrix_elements():
-            uid += 1 # update the identification number
-            me.get('processes')[0].set('uid', uid)
-
-        cpu_time2 = time.time()
-        return ndiags, cpu_time2 - cpu_time1        
-
-    def add_content_to_global_ME7_resources(self, global_ME7_dir, **opts):
-        """ Pass the MadLoopParams.dat card to the global ME7 resources."""
-        super(Contribution_V, self).add_content_to_global_ME7_resources(global_ME7_dir, **opts)
-        
-        destination = pjoin(global_ME7_dir,'Cards','MadLoopParams.dat')
-        if not os.path.exists(destination):
-            shutil.copyfile(pjoin(self.export_dir,'Cards','MadLoopParams.dat'), destination)
-
-    def set_helas_model(self):
-        """ Instantiate properly the helas model """
-
-        assert self.exporter.exporter == 'v4'
-        assert (not self.options['_model_v4_path'])
-        self.helas_model = helas_call_writers.FortranUFOHelasCallWriter(self.model)
 
     def get_all_necessary_integrated_currents(self, all_MEAccessors):
         """ Given the list of currents already available encoded in the all_MEAccessors,
@@ -1525,14 +1876,29 @@
         return all_currents
 
     @classmethod
-    def add_current_accessors(
-        cls, model, all_MEAccessors, root_path, current_set, currents_to_consider ):
-        """Generate and add all integrated current accessors to the MEAccessorDict.
+    def remove_counterterms_with_no_reduced_process(cls, all_MEAccessors, counterterms):
+        """Given the list of available reduced processes encoded in the MEAccessorDict 'all_MEAccessors'
+        given in argument, remove all the counterterms whose underlying reduced process does not exist.
         For now we can recycle the implementation of the Contribution_R class.
         """
-
-        return Contribution_R.add_current_accessors(
-            model, all_MEAccessors, root_path, current_set, currents_to_consider)
+        return Contribution_R.remove_counterterms_with_no_reduced_process(all_MEAccessors, counterterms)
+
+    def remove_integrated_counterterms_set_to_zero(self, all_ME_accessors):
+        """ Remove all integrated counterterms involving currents whose implementation has set
+        them to zero."""
+
+        for process_key, counterterms in self.integrated_counterterms.items():
+            for integrated_counterterm_properties in list(counterterms):
+                integrated_counterterm = integrated_counterterm_properties['integrated_counterterm']
+                # misc.sprint("Considering integrated CT %s" % str(integrated_counterterm))
+                if integrated_counterterm.is_singular():
+                    for current in integrated_counterterm.get_all_currents():                            
+                        accessor, _ = all_ME_accessors[current]
+                        if accessor.subtraction_current_instance.is_zero:
+                            # misc.sprint("Zero current found in integrated CT %s" % str(integrated_counterterm))
+                            logger.debug("Zero current found in integrated CT %s, removing it now."%str(integrated_counterterm))
+                            counterterms.remove(integrated_counterterm_properties)
+                            break
 
     def get_integrands_for_process_map(self, process_map, model, run_card, all_MEAccessors, ME7_configuration):
         """ Returns all the integrands implementing this contribution for the specified process_map.
@@ -1540,6 +1906,7 @@
         """
         
         relevant_counterterms = {}
+        self.remove_integrated_counterterms_set_to_zero(all_MEAccessors)
         if self.integrated_counterterms:
             for process_key in process_map:
                 relevant_counterterms[process_key] = self.integrated_counterterms[process_key]
@@ -1561,10 +1928,20 @@
         # Get the basic accessors for the matrix elements
         super(Contribution_V, self).add_ME_accessors(all_MEAccessors, root_path)
 
+        # Only initialize all_counterterms_removed to True if there is CT to removed to begin with
+        all_counterterms_removed = any(len(CTs)>0 for CTs in self.integrated_counterterms.values())
         for process_key, CT_properties in self.integrated_counterterms.items():
             # Remove integrated counterterms with non-existing underlying Born processes
-            self.remove_counterterms_with_no_reduced_process(all_MEAccessors, CT_properties)
-
+            if not self.remove_counterterms_with_no_reduced_process(all_MEAccessors, CT_properties):
+                all_counterterms_removed = False
+ 
+        # Sanity check
+        if all_counterterms_removed:
+            logger.warning("All integrated counterterms from contribution '%s' were removed because"%self.short_name()+
+                " their corresponding reduced processes was not found in the list of processes exported.\n"+
+                " This is likely wrong and originates from a mismatch in the process definitions of the various"+
+                " contributions building this higher order computation.")
+        
         # Obtain all necessary currents
         current_set = self.options['subtraction_currents_scheme']
         currents_to_consider = self.get_all_necessary_integrated_currents(all_MEAccessors)
@@ -1576,7 +1953,6 @@
         """ Figures out the permutation to apply to map the origin pdg orders
         (a 2-tuple of initial and final state pdgs) to the destination pdg orders (same format)."""
 
-
         # Create a look_up list from the user-provided list of PDGs, whose elements will progressively be set to zero
         # as they are being mapped
         look_up_list = list(list(a_tuple) for a_tuple in origin_pdg_orders)
@@ -1721,6 +2097,77 @@
         
         return all_mappings
 
+    def combine_initial_state_counterterms(self):
+        """ The initial state beam factorization current, including the integrated collinear
+        ones, return a flavor matrix that represents all possible backward evolution of the
+        initial state. For instance, if the BF1 process is:
+              d d~ > z
+        then there will be two [C(1,3)] coming from the real processes d d~ > z g and g d~ > z
+        which correspond to the two possible backward evolutions to a d-quark and a gluon
+        respectively. Given that those are accounted for at once in the beam-factorization
+        current via the flavor matrix, we must make sure to combine all such integrated
+        initial-state counterterms so as to avoid double-counting.
+        Alternatively, one can disable this combination and always backward evolve only to 
+        the particular initial-state flavor specified in this particular integrated counterterm.
+        
+        WARNING: Eventually this entire construction of the 'allowed_backward_evolved_flavors'
+        should be dropped in favour of further differentiation of the integrated collinear ISR
+        current that should apply only to particular current (for instance the one 
+        *only* backward-evolving to a gluon). When this will be done, we will be able to 
+        drop this construction altogether.
+        """
+
+        for process_key, (process, mapped_processes) in self.get_processes_map().items():
+            # It may be that a particular process does not have any counterterms (typically when
+            # considering forcing a beam-type by hand at generation time for debugging purposes, for 
+            # example:
+            #    generate e+ e- > j j j / a QED=2 --NLO=QCD --beam_types=proton@(1,-1,21)
+            # In this case, self.integrated_counterterms would be empty and this combination 
+            # should be skipped.
+            if process_key not in self.integrated_counterterms:
+                continue
+            for counterterm_characteristics in self.integrated_counterterms[process_key]:
+                # One should perform here the combination of the values of each key in
+                # self.integrated_counterterms.
+            
+                # By default include all possible bbackward evolutions.
+                counterterm_characteristics['allowed_backward_evolved_flavors'] = \
+                                                        {'beam_one': 'ALL', 'beam_two': 'ALL',}
+                if _COMBINE_INITIAL_STATE_INTEGRATED_COLLINEAR_COUNTERTERMS:
+                    raise MadGraph5Error('Combination of integrated initial-state collinear '+
+                                                                    'counterterms not implemented yet')
+                else:
+                    # Since in this case we don't combine different initial state integrated
+                    # collinear counterterms that differ only by the backward evolved flavor,
+                    # we must restrict each one to the particular backward evolved flavor they
+                    # correspond to.
+                    integrated_counterterm = counterterm_characteristics['integrated_counterterm']
+                    if integrated_counterterm.does_require_correlated_beam_convolution():
+                        # In this case the beam currents has a diagonal flavor matrix anyway
+                        continue
+                    beam_number_map = {'beam_one': 0, 'beam_two': 1}
+                    beam_currents = integrated_counterterm.get_beam_currents()
+                    active_beams = set([])
+                    # Check if the counterterm contains any current that could have a 
+                    # flavor matrix attached to either beam_one of beam_two
+                    for bc in beam_currents:
+                        for beam_name in beam_number_map:
+                            # If it is a pure integrated mass collinear counterterm, then
+                            # no masking should be applied
+                            if (bc[beam_name] is None) or all(ss.name()=='F' for ss in 
+                                      bc[beam_name].get('singular_structure').decompose()):
+                                continue
+                            # Now we can add this beam as a one for which a masking of the
+                            # allowed initial state flavours must be specified.
+                            active_beams.add(beam_name)
+                    # Now for each beam identified as potentially having a flavour matrix
+                    # then identify which flavours must be allowed as allowed ones after
+                    # the application of the flavor matrix
+                    for beam_name in active_beams:
+                        counterterm_characteristics['allowed_backward_evolved_flavors'][beam_name] = \
+                            tuple(set(fl[0][beam_number_map[beam_name]] for fl in 
+                                counterterm_characteristics['resolved_flavors_combinations']))
+
     def add_integrated_counterterm(self, integrated_CT_properties):
         """ Virtual contributions can receive integrated counterterms and they will
         be stored in the attribute list self.integrated_counterterms."""
@@ -1728,6 +2175,7 @@
         # Extract quantities from integrated_counterterm_properties
         integrated_counterterm = integrated_CT_properties['integrated_counterterm']
         # flavors_combinations = integrated_CT_properties['flavors_combinations']
+
         
         # Sort PDGs only when processes have been grouped
         sort_PDGs = self.group_subprocesses
@@ -1748,23 +2196,42 @@
         # Use ordered PDGs since this is what is used in the inverse_processes_map
         # Also overwrite some of the process properties to make it match the 
         # loop process definitions in this contribution
+        n_loops_of_reduced_process = self.contribution_definition.n_loops - \
+            integrated_counterterm.get_n_loops_in_beam_factorization_of_host_contribution()
+        
         counterterm_reduced_process_key = ProcessKey(integrated_counterterm.process,
-            sort_PDGs=sort_PDGs, n_loops=1, NLO_mode='virt', 
-            perturbation_couplings=self.contribution_definition.correction_couplings)\
-                                                                       .get_canonical_key()
+            sort_PDGs   = sort_PDGs,
+            n_loops     = n_loops_of_reduced_process).get_canonical_key()
         
         # misc.sprint(inverse_processes_map.keys(), len(inverse_processes_map.keys()))
         # misc.sprint(counterterm_reduced_process_key)
-        if counterterm_reduced_process_key not in inverse_processes_map:
+
+        # Now look for a matching process key using only a limited set of criteria
+        potential_matches = []
+        criteria_considered = ['PDGs','n_loops','id']
+        for process_key in inverse_processes_map:
+            if all(pk[1] == crpk[1] for (pk, crpk) in zip(process_key,
+                           counterterm_reduced_process_key) if pk[0] in criteria_considered):
+                # Now we can simply add this integrated counterterm in the group of the
+                # defining process to which the reduced process is mapped
+                potential_matches.append(inverse_processes_map[process_key])
+
+        if len(potential_matches)==0:
             # The reduced process of the integrated counterterm is not included in this
             # contribution so we return here False, letting the ME7 exporter know that
-            # we cannot host this integrated counterterm.
+            # we cannot host this integrated counterterm here.
             return False
         
-        # Now we can simply add this integrated counterterm in the group of the
-        # defining process to which the reduced process is mapped
-        defining_key, virtual_process_instance = \
-                                     inverse_processes_map[counterterm_reduced_process_key]
+        if len(potential_matches)>1:
+            logger.warning(
+                ('MadNkLO found the following processes:\n%s\nto be potential matches'+
+                 'to this integrated counterterm:\n%s\nThe first matching process will be'+
+                 ' selected, but this may indicate that the current list of criteria considered'+
+                 ' (%s) should be extended.')%(
+                     '\n'.join(p[1].nice_string() for p in potential_matches),
+                     integrated_counterterm.nice_string(), ','.join(criteria_considered)) )
+
+        defining_key, virtual_process_instance = potential_matches[0]
         
         integrated_counterterm_properties = dict(integrated_CT_properties)
        
@@ -1795,6 +2262,23 @@
         # Obtain the corresponding permutation
         basic_permutation = self.get_basic_permutation(origin_pdgs_order,destination_pdgs_order)
 
+        # If the integrated counterterm requires only a single active beam factorization then
+        # this beam-factorization contribution should indeed host this counterterm only if:
+        #  a) The initial state legs *does* require swapping and the factorized beam number of
+        #     this integrated CT does *not* match the factorized beam number of thhis contribution
+        #  b) The initial state legs does *not* require swapping and the factorized beam number 
+        #     of this integrated CT *does* match the factorized beam number of this contribution
+        active_beams = {k:self.contribution_definition.is_beam_active(k) for k in ['beam_one','beam_two'] }
+        integrated_CT_active_beams = integrated_counterterm.get_necessary_beam_convolutions()
+        if active_beams.values().count(True)==1 and len(integrated_CT_active_beams)==1:
+            require_initial_state_swapping = ( basic_permutation[0] == 1 and 
+                                               basic_permutation[1] == 0 )
+            integrated_CT_active_beam = list(integrated_CT_active_beams)[0]
+            if (require_initial_state_swapping and active_beams[integrated_CT_active_beam]==True) or \
+                ((not require_initial_state_swapping) and active_beams[integrated_CT_active_beam]==False):
+                # Then this integrated CT should not be hosted in this contribution
+                return False
+        
         # Now we need to distribute all flavors that are parents of the integrated current
         # in all possible ways. In the example below, the parent of C(3,7) is a d-quark and
         # we must assign it as both leg #3 and leg #5 of the following reduced integrated
@@ -1822,11 +2306,13 @@
 
         # Now obtain all possible ways of distributing these parent flavors over the
         # flavors of the underlying reduced process.
+        # Initial-state are always distinguishible (and they were treated as so when 
+        # generating the list of integrated counterterms to dispatch) so that we must not
+        # distribute over identical initial_state_parent_flavors:
+        initial_state_parent_flavors = {}
         flavor_permutations = self.distribute_parent_flavors(
-            (initial_state_parent_flavors, final_state_parent_flavors),
-            destination_pdgs_order)
-#        misc.sprint('For integrated counterterm %s, flavor permutations are: %s'%(str(integrated_counterterm),str(flavor_permutations)))
-        
+            (initial_state_parent_flavors, final_state_parent_flavors), destination_pdgs_order)
+
         # Now combine these permutations with the basic permutation to obtain the final
         # list of permutations to consider
         for flavor_permutation in flavor_permutations:
@@ -1835,6 +2321,17 @@
                 combined_permuation[index] = basic_permutation[flavor_permutation.get(index, index)]
             permutations.append(combined_permuation)
 
+#       Example of how to print out information of a specific contribution.
+#        if isinstance(self, Contribution_RV) and False:
+#            misc.sprint('-'*50)
+#            misc.sprint(momenta_dict)
+#            misc.sprint(initial_state_parent_flavors)
+#            misc.sprint(final_state_parent_flavors)
+#            misc.sprint(destination_pdgs_order)            
+#            misc.sprint(integrated_counterterm.process.nice_string())
+#            misc.sprint('For integrated counterterm %s, flavor permutations are: %s'%(str(integrated_counterterm),str(permutations)))
+#            misc.sprint('-'*50)
+
         # Store the mapping to apply to the virtual ME inputs
         integrated_counterterm_properties['input_mappings'] = permutations
         
@@ -1858,69 +2355,177 @@
         # Let the ME7 exporter that we could successfully host this integrated counterterm
         return True
 
-    def process_dir_name(self, process):
-        """ Given a specified process, return the directory name in which it is output."""
-        
-        # For MadLoop ME's there is extra prefixes to avoid symbol and resource file clashes.
-        return "%d_%s_%s"%(process.get('id'), process.get('uid'),
-                                                      process.shell_string(print_id=False))
-
-
-    def generate_code(self):
-        """ Assuming the Helas Matrix Elements are now generated, we can write out the corresponding code."""
-
-        matrix_elements = self.all_matrix_elements.get_matrix_elements()
-            
-        cpu_time_start = time.time()
-
-        # Pick out the matrix elements in a list
-        matrix_elements = self.all_matrix_elements.get_matrix_elements()
-        # MadLoop standalone fortran output
-        calls = 0
-        for me in matrix_elements:
-            # Choose the group number to be the unique id so that the output prefix
-            # is nicely P<proc_id>_<unique_id>_
-            calls = calls + self.exporter.generate_loop_subprocess(me, 
-                self.helas_model,
-                group_number = me.get('processes')[0].get('uid'),
-                proc_id = None,
-                config_map=None,
-                unique_id=me.get('processes')[0].get('uid'))
-
-        # If all ME's do not share the same maximum loop vertex rank and the
-        # same loop maximum wavefunction size, we need to set the maximum
-        # in coef_specs.inc of the HELAS Source. The SubProcesses/P* directory
-        # all link this file, so it should be properly propagated
-        if self.options['loop_optimized_output'] and len(matrix_elements)>1:
-            max_lwfspins = [m.get_max_loop_particle_spin() for m in \
-                                                            matrix_elements]
-            max_loop_vert_ranks = [me.get_max_loop_vertex_rank() for me in \
-                                                            matrix_elements]
-            if len(set(max_lwfspins))>1 or len(set(max_loop_vert_ranks))>1:
-                self.exporter.fix_coef_specs(max(max_lwfspins),\
-                                                   max(max_loop_vert_ranks))
-        
-        return calls, time.time() - cpu_time_start
-
-
 class Contribution_RV(Contribution_R, Contribution_V):
-    """ Implements the handling of a real-virtual type of contribution."""
-    pass
+    """ Implements the general handling of contribution with arbitrary number of reals, virtuals
+    and beam factorization contributions."""
+    
+    def get_integrands_for_process_map(self, process_map, model, run_card, all_MEAccessors, ME7_configuration):
+        """ Returns all the integrands implementing this contribution for the specified process_map.
+        The instance of MEAccessorDict is necessary so as to be passed to the integrand instances.
+        """
+        
+        relevant_counterterms = {}
+        self.remove_local_counterterms_set_to_zero(all_MEAccessors)
+        if self.counterterms:
+            for process_key in process_map:
+                relevant_counterterms[process_key] = self.counterterms[process_key]
+        
+        relevant_integrated_counterterms = {}
+        self.remove_integrated_counterterms_set_to_zero(all_MEAccessors)
+        if self.integrated_counterterms:
+            for process_key in process_map:
+                relevant_integrated_counterterms[process_key] = self.integrated_counterterms[process_key]
+
+        return [ ME7_integrands.ME7Integrand(model, run_card,
+                   self.contribution_definition,
+                   process_map,
+                   self.topologies_to_processes,
+                   self.processes_to_topologies,
+                   all_MEAccessors,
+                   ME7_configuration,
+                   subtraction_mappings_scheme=self.options['subtraction_mappings_scheme'],
+                   counterterms=relevant_counterterms,
+                   integrated_counterterms=relevant_integrated_counterterms
+                ) ]
+
+    def get_additional_nice_string_printout_lines(self, format=0):
+        """ Return additional information lines for the function nice_string of this contribution."""
+        res = Contribution_R.get_additional_nice_string_printout_lines(self, format=format)
+        res.extend(Contribution_V.get_additional_nice_string_printout_lines(self, format=format))
+        return res
+
+    def get_nice_string_process_line(self, process_key, defining_process, format=0):
+        """ Return a nicely formated process line for the function nice_string of this 
+        contribution."""
+        GREEN = '\033[92m'
+        ENDC = '\033[0m'        
+        res = GREEN+'  %s'%(defining_process.nice_string(print_weighted=False)
+                                                             .replace('Process: ',''))+ENDC
+
+        if self.integrated_counterterms is None or self.counterterms is None:
+            return res
+
+        if format<2:
+            if process_key in self.counterterms:
+                res += ' | %d local counterterms'%len([ 1
+                    for CT in self.counterterms[process_key] if CT.is_singular() ])
+            else:
+                res += ' | 0 local counterterm'                
+            if process_key in self.integrated_counterterms:
+                res += ' and %d integrated counterterms'%len(self.integrated_counterterms[process_key])
+            else:
+                res += ' and 0 integrated counterterm'
+        else:
+            long_res = [' | with the following local and integrated counterterms:']
+            for CT in self.counterterms[process_key]:
+                if CT.is_singular():
+                    if format==2:
+                        long_res.append( '   | %s' % str(CT))
+                    elif format==3:
+                        long_res.append( '   | %s' % CT.__str__(
+                            print_n=True, print_pdg=True, print_state=True ) )
+                    elif format>3:
+                        long_res.append(CT.nice_string("   | "))
+            for CT_properties in self.integrated_counterterms[process_key]:
+                CT = CT_properties['integrated_counterterm']
+                if format==2:
+                    long_res.append( '   | %s' % str(CT))
+                elif format==3:
+                    long_res.append( '   | %s' % CT.__str__(
+                        print_n=True, print_pdg=True, print_state=True ))
+                elif format==4:
+                    long_res.append(CT.nice_string("   | "))
+                elif format>4:
+                    long_res.append(CT.nice_string("   | "))
+                    for key, value in CT_properties.items():
+                        if not key in ['integrated_counterterm', 'matching_process_key']:
+                            long_res.append( '     + %s : %s'%(key, str(value)))
+            res += '\n'.join(long_res)
+
+        return res
+
+
+class Contribution_BS(Contribution_RV):
+    """ A class implementing the Beam-Soft factorization contributions (BS) originating
+    from soft counterterms in the colorful currents scheme which can recoil against the
+    initial states, as when using the SoftBeamsRecoilNLOWalker mappings."""
+
+    def __init__(self, contribution_definition, cmd_interface, **opts):
+        """ Make sure that this contribution is instantiated by a contribution definition
+        that enforces a correlated two-beams convolution."""
+        
+        if not contribution_definition.correlated_beam_convolution:
+            raise MadGraph5Error("A soft beam convolution contribution must be instantiated"+
+                " from a contribution definition that specifies a *correlated* convolution "+
+                "of the two beams.")
+        super(Contribution_RV,self).__init__(contribution_definition, cmd_interface, **opts)
+
+    def add_beam_factorization_processes_from_contribution(self, contrib, beam_factorization_order):
+        """ 'Import' the processes_map and all generated attributes of the contribution
+        passed in argument to assign them to this beam factorization contributions. Contrary
+        to the implementation of this function in the base class, here we do not instantiate
+        'bulk' beam factorization currents because they are not needed. The BS contribution is
+        only meant to be a receptacle for the integrated soft counterterm in the 'colorful'
+        currents scheme when the 'ppToOneWalker' mappings scheme is used (where the soft
+        recoils equally against both initial states).
+        """
+
+        # Make sure the source contribution does have a factorizing beam for the active
+        # ones of this contribution.
+        if ( contrib.contribution_definition.beam_factorization['beam_one'] is None ) or \
+           ( contrib.contribution_definition.beam_factorization['beam_two'] is None ):
+            return
+    
+        # Now fetch the process map to load from
+        source_processes_map = contrib.get_processes_map()
+        processes_map = self.processes_map[1]
+        # To the new contribution instances, also assign references to the generated attributes
+        # of the original contributions (no deep copy necessary here)
+        for process_key, (defining_process, mapped_processes) in source_processes_map.items():
+            
+            # First add this process to the current processes map
+            if process_key in processes_map:
+                raise MadGraph5Error('MadNkLO attempts to add beam factorization terms for'+
+                    ' the %s a second time. This should never happen.'%(
+                             defining_process.nice_string().replace('Process','process')))
+            processes_map[process_key] = (
+                  defining_process.get_copy(), [mp.get_copy() for mp in mapped_processes] )
+            
+            # Contrary to what is done in the base class, no 'bulk' beam factorization 
+            # current need to be added here
+            
+            # Also add the corresponding topologies
+            self.import_topologies_from_contrib(process_key, contrib)
+
+    def generate_all_counterterms(self, ignore_integrated_counterterms=False, group_processes=True, *args, **opts):
+        """ Contrary to the other contributions the BS one has no physical contribution that require subtraction
+        and is only a place-holder to receive the integrated soft counterterms. Therefore nothing needs to be
+        done a this stage here. """
+    
+        for process_key, (defining_process, mapped_processes) in self.get_processes_map().items():
+            self.counterterms[process_key] = []
+        all_integrated_counterterms = []
+        return all_integrated_counterterms
 
 class ContributionList(base_objects.PhysicsObjectList):
     """ A container for storing a list of contributions."""
     
     contributions_natural_order = [
-        ('LO',    (Contribution_B, Contribution_LIB) ),
-        ('NLO',   (Contribution_R, Contribution_V) ),
-        ('NNLO',  (Contribution_RR, ) ),
-        ('NNNLO', (Contribution_RRR, ) )
+        ('LO',    (Contribution_B,) ),
+        ('NLO',   (Contribution_R, Contribution_V, Contribution_RV) ),
+        ('NNLO',  (Contribution_RR, Contribution_RV) ),
+        ('NNNLO', (Contribution_RRR, Contribution_RV) )
     ]
     
     def is_valid_element(self, obj):
         """Test if object obj is a valid instance of Contribution."""
         return isinstance(obj, Contribution)
     
+    def get_loop_induced_contributions(self):
+        """ Return a list of loop-induced contributions."""
+        return ContributionList([contrib for contrib in self if
+                                        contrib.contribution_definition.is_loop_induced()])
+    
     def get_contributions_of_order(self, correction_order):
         """ Returns a list of all contributions of a certain correction_order in argument."""
         return ContributionList([contrib for contrib in self if
@@ -1933,13 +2538,19 @@
                 correction_classes = tuple(correction_classes)
             else:
                 correction_classes = (correction_classes,)                
-        return ContributionList([contrib for contrib in self if isinstance(contrib, correction_classes)])
+        return ContributionList([contrib for contrib in self if type(contrib) in correction_classes])
 
     def nice_string(self, format=0):
         """ A nice representation of a list of contributions. 
         We can reuse the function from ContributionDefinitions."""
         return base_objects.ContributionDefinitionList.contrib_list_string(self, format=format)
 
+    def get_max_correction_order(self):
+        """ Returns the maximum correction order found  among these contributions."""
+        
+        return 'N'*(max(contrib.contribution_definition.correction_order.count('N')
+                                                                 for contrib in self))+'LO'
+
     def sort_contributions(self):
         """ Sort contributions according to the order dictated by the class_attribute contributions_natural_order"""
         
@@ -1948,7 +2559,14 @@
             for contrib_type in contribution_types:
                 selected_contribs = self.get_contributions_of_order(correction_order).\
                                         get_contributions_of_type(contrib_type)
-                new_order.extend(selected_contribs)
+                # Now sort contributions of this type using the ordering F1, F2, F1F2
+                factorization_order = [(False,False),(True, False),(False, True), (True, True)]
+                new_order.extend(sorted(selected_contribs,
+                    key = lambda contrib: factorization_order.index((
+                        contrib.contribution_definition.is_beam_active('beam_one'),
+                        contrib.contribution_definition.is_beam_active('beam_two')
+                    ))
+                ))
                 for contrib in selected_contribs:
                     self.pop(self.index(contrib))
     
@@ -1961,7 +2579,6 @@
         
         # Keep track of the return values
         return_values = []
-        
         remaining_contribs = list(self)
         for correction_order, contribution_types in self.contributions_natural_order:
             for contrib_type in contribution_types:
@@ -1992,13 +2609,19 @@
             if log: logger.info('%s the %d contribution%s of unknown type...'%(log,
                         len(remaining_contribs), 's' if len(remaining_contribs)>1 else ''))                    
         for i, contrib in enumerate(remaining_contribs):
-            if log: logger.info('%s %d/%d'%(contrib_type.__name__, i+1, len(remaining_contribs)))
+            if log: logger.info('%s %d/%d'%(type(contrib).__name__, i+1, len(remaining_contribs)))
             try:
                 contrib_function = getattr(contrib, method)
             except AttributeError:
                 raise MadGraph5Error("The contribution\n%s\n does not have function '%s' defined."%(
                                                                             contrib.nice_string(), method))
-            return_values.append((contrib, contrib_function(*method_args, **method_opts)))
+            # If method opts depend on the contribution, it can be passed as a callable
+            # that can be evaluated here
+            if not isinstance(method_opts, dict) and callable(method_opts):
+                method_options = method_opts(contrib)
+            else:
+                method_options = dict(method_opts)
+            return_values.append((contrib, contrib_function(*method_args, **method_options)))
         
         return return_values
 
@@ -2006,9 +2629,10 @@
 # by the interface when using a PLUGIN system where the user can define his own class of Contribution.
 # Notice that this Contribution must be placed after all the Contribution daughter classes have been declared.
 Contribution_classes_map = {'Born': Contribution_B,
-                            'LoopInduced_Born': Contribution_LIB,
                             'Virtual': Contribution_V,
                             'SingleReals': Contribution_R,
+                            'RealVirtual': Contribution_RV,
                             'DoubleReals': Contribution_RR,
                             'TripleReals': Contribution_RRR,
+                            'BeamSoft': Contribution_BS,
                             'Unknown': None}

=== modified file 'madgraph/core/subtraction.py'
--- madgraph/core/subtraction.py	2018-06-19 17:36:59 +0000
+++ madgraph/core/subtraction.py	2019-01-14 20:53:09 +0000
@@ -25,6 +25,8 @@
 import os
 import importlib
 import shutil
+import copy
+import traceback
 
 if __name__ == '__main__':
     sys.path.append(os.path.join(os.path.dirname(
@@ -43,6 +45,12 @@
 logger = logging.getLogger('madgraph')
 pjoin = os.path.join
 
+# Specify here the list of subtraction currents schemes and subtraction_mappings_scheme
+# that uses soft counterterms that recoils against the initial states. This is necessary in
+# order to decide when a contribution with a correlated convolution of both beams must be setup.
+_currents_schemes_requiring_soft_beam_factorization = ['colorful', 'colorful_pp']
+_mappings_schemes_requiring_soft_beam_factorization = ['SoftBeamsRecoilNLO', ]
+
 #=========================================================================================
 # Multinomial function
 #=========================================================================================
@@ -347,6 +355,7 @@
 class SingularStructure(object):
     """Object that represents a hierarchical structure of IR singularities."""
 
+
     def __init__(self, *args, **opts):
         """Initialize a hierarchical singular structure."""
 
@@ -399,6 +408,26 @@
         tmp_str += ")"
         return tmp_str
 
+    def does_require_correlated_beam_convolution(self):
+        """ Checks whether this integrated counterterm requires a host contribution featuring
+        correlated convolution of the beams (i.e. 'BS', 'VS', etc... contributions).
+        This is the case only for integrated counterterms with disjoint structures having each
+        at least one soft *BeamCurrents* that originated from a colorful subtraction currents 
+        scheme with mappings such as ppToOneWalker that recoils the soft momentum equally 
+        against both initial-state beams."""
+
+        # One probably needs to do something more refined for higher orders, since the
+        # structure C(3,4,S(5)) probably requires a collinear mapping and not a soft one.
+        # However, maybe for such cases the mappings yielding correlated beam convolutions
+        # can be avoided altogether since one would be guaranteed to have final-state recoilers
+        # for the singly-unresolved limit S(5).
+        # Implement a caching too since this is called by the ME7Integrand at runtime.
+        if not hasattr(self,'does_require_correlated_beam_convolution_cached'):
+            self.does_require_correlated_beam_convolution_cached = \
+                all(any(c.name()=='S' for c in sub_ss.decompose()) for sub_ss in self.substructures)
+        
+        return self.does_require_correlated_beam_convolution_cached
+
     def get_subtraction_prefactor(self):
         """Determine the prefactor due to the nested subtraction technique."""
         
@@ -420,6 +449,14 @@
         
         return pref
 
+    def get_beam_factorization_legs(self):
+        """ Return the set of all the beam factorization legs involved in this singular
+        structure. Note that beam factorization structures should always be at the top level
+        and never nested."""
+        return set(sum([
+            [ l.n for l in struct.get_all_legs() ] 
+                for struct in self.substructures if isinstance(struct, BeamStructure)],[]))
+
     def annihilate(self):
         """When an operator cannot act on this structure,
         remove this structure altogether
@@ -454,6 +491,14 @@
             total_unresolved += substructure.count_unresolved()
         return total_unresolved
 
+    def count_couplings(self):
+        """Count the number of couplings covered by this structure."""
+
+        total_couplings = 0
+        for substructure in self.substructures:
+            total_couplings += substructure.count_couplings()
+        return total_couplings
+
     def get_canonical_representation(self, track_leg_numbers=True):
         """Creates a canonical hashable representation of self."""
         
@@ -508,12 +553,6 @@
         inner.append(foo)
         return inner
 
-    name_dictionary = {
-        '' : 'SingularStructure',
-        'S': 'SoftStructure',
-        'C': 'CollStructure'
-    }
-
     @staticmethod
     def from_string(string, process, log=False):
         """Reconstruct a singular structure from a string representation, e.g.
@@ -522,13 +561,20 @@
         :return: The reconstructed singular structure if successful, else None.
         """
 
+        singular_structures_name_dictionary = {
+            '' : SingularStructure,
+            'S': SoftStructure,
+            'C': CollStructure,
+            'F': BeamStructure
+        }
+
         msg = "SingularStructure.from_string:"
         openpos = string.find('(')
         closepos = string.rfind(')')
         name = string[:openpos]
-        if not SingularStructure.name_dictionary.has_key(name):
+        if not singular_structures_name_dictionary.has_key(name):
             if log:
-                print msg, "no singular structure with name", name
+                logger.warning('%s %s: %s'%(msg, "no singular structure with name", name))
             return None
         legs = []
         substructures = []
@@ -558,7 +604,7 @@
                     leg_n = int(before_comma)
                 except:
                     if log:
-                        print msg, "could not convert", before_comma, "to integer"
+                        logger.warning('%s %s %s'%("could not convert", before_comma, "to integer"))
                     return None
                 for leg in process['legs']:
                     if leg['number'] == leg_n:
@@ -567,33 +613,45 @@
                         break
                 if before_comma:
                     if log:
-                        print msg, "no leg number", leg_n, "found"
+                        logger.warning('%s %s %d %s'%(msg, "no leg number", leg_n, "found"))
                     return None
         if before_comma:
             if log:
-                print msg, "bracket matching failed for", before_comma[:-1]
+                logger.warning("%s %s '%s'"%(msg, "bracket matching failed for:", before_comma[:-1]))
             return None
-        return eval(SingularStructure.name_dictionary[name])(
-            substructures=substructures, legs=legs)
+        return singular_structures_name_dictionary[name](substructures=substructures, legs=legs)
+
+class BeamStructure(SingularStructure):
+
+    def count_couplings(self):
+        return 1
+
+    def count_unresolved(self):
+        return 0
+
+    def name(self):
+        return "F"
 
 class SoftStructure(SingularStructure):
 
+    def count_couplings(self):
+        return self.count_unresolved()
+
     def count_unresolved(self):
-
         return len(self.get_all_legs())
 
     def name(self):
-
         return "S"
     
 class CollStructure(SingularStructure):
 
+    def count_couplings(self):
+        return self.count_unresolved()
+
     def count_unresolved(self):
-
         return len(self.get_all_legs())-1
 
     def name(self):
-
         return "C"
     
 #=========================================================================================
@@ -762,6 +820,30 @@
         # WARNING: this statement has to be checked
         # By default, not skipping overlaps ensures correct subtraction
         return True
+    
+class BeamOperator(SingularOperator):
+    """Object that represents a beam factorization elementary singular operator."""
+
+    def name(self):
+
+        return "F"
+
+    def get_structure(self):
+
+        return BeamStructure(legs=self)
+
+    def __init__(self, *args, **opts):
+
+        super(BeamOperator, self).__init__(*args, **opts)
+
+    def act_here_needed(self, structure):
+
+        assert isinstance(structure, SingularStructure)
+
+        # Beam factorization limits should never be needed within any other singular structure
+        # WARNING: this statement has to be checked
+        
+        return False
 
 #=========================================================================================
 # SingularOperatorList
@@ -796,7 +878,6 @@
 #=========================================================================================
 # Current
 #=========================================================================================
-
 class Current(base_objects.Process):
 
     def default_setup(self):
@@ -824,6 +905,11 @@
 
         return self['singular_structure'].count_unresolved()
 
+    def count_couplings(self):
+        """Count the number of couplings expected in this current."""
+
+        return self['singular_structure'].count_couplings()
+
     def discard_leg_numbers(self,discard_initial_leg_numbers=True):
         """Discard all leg numbers in the singular_structure
         and in the parent_subtraction_leg,
@@ -852,31 +938,63 @@
             readable_string += " " + str(self.get('orders'))
         return readable_string
 
+    def does_require_correlated_beam_convolution(self):
+        """ Checks whether this integrated counterterm requires a host contribution featuring
+        correlated convolution of the beams (i.e. 'BS', 'VS', etc... contributions).
+        This is the case only for integrated counterterms with disjoint structures having each
+        at least one soft *BeamCurrents* that originated from a colorful subtraction currents 
+        scheme with mappings such as ppToOneWalker that recoils the soft momentum equally 
+        against both initial-state beams."""
+
+        # For basic currents, this is never the case
+        return False
+
     def get_reduced_flavors(self, defining_flavors, IR_subtraction_module, routing_dict):
         """Given the defining flavors dictionary specifying the flavors of some leg numbers,
         add entries corresponding to the subtraction legs in this node and using the function
         get_parent_PDGs from the subtraction module."""
         
-        def get_parent(PDGs):
+        get_particle = IR_subtraction_module.model.get_particle
+        
+        def get_parent(PDGs, is_initial):
             res = IR_subtraction_module.parent_PDGs_from_PDGs(PDGs)
             if len(res)!=1:
-                raise MadGraph5Error("Multiple mother PDGs assignment not supported"+
-                                            " yet by the function get_reduced_flavors")
-            return res[0]
+                str1 = "Multiple mother PDGs assignment "
+                str2 = " is not yet supported by the function get_reduced_flavors"
+                misc.sprint("get_parent called with PDGs = " + str(PDGs))
+                raise MadGraph5Error(str1 + str(res) + str2)
+            # Now flip back the identity of the parent PDG if in the initial state
+            if is_initial:
+                return get_particle(res[0]).get_anti_pdg_code()
+            else:
+                return res[0]
         
         all_legs = self['singular_structure'].get_all_legs()
+        
+        # If all legs are unresolved, as it is the case in the pure soft configurations,
+        # then nothing needs to be done
+        if self.count_unresolved()==len(all_legs):
+            return
+        
         leg_numbers = [leg.n for leg in all_legs]
-        get_particle = IR_subtraction_module.model.get_particle
-        # Swap the flavors of the initial states
-        leg_flavors = [ ( 
-                get_particle(defining_flavors[leg.n]).get_anti_pdg_code() if 
-                leg.state==SubtractionLeg.INITIAL else get_particle(defining_flavors[leg.n]).get_pdg_code() ) 
-            for leg in all_legs ]
 
-        if len(leg_numbers)>0:            
-            mother_number = Counterterm.get_ancestor(frozenset(leg_numbers), routing_dict)
-            if mother_number not in defining_flavors:
-                defining_flavors[mother_number] = get_parent(leg_flavors)
+        try:
+            # Swap the flavors of the initial states
+            leg_flavors = [ ( 
+                    get_particle(defining_flavors[leg.n]).get_anti_pdg_code() if 
+                    leg.state==SubtractionLeg.INITIAL else get_particle(defining_flavors[leg.n]).get_pdg_code() ) 
+                for leg in all_legs ]
+    
+            if len(leg_numbers)>0: 
+                mother_number = Counterterm.get_ancestor(frozenset(leg_numbers), routing_dict)
+                if mother_number not in defining_flavors:
+                    defining_flavors[mother_number] = get_parent(leg_flavors, 
+                                   any(leg.state==SubtractionLeg.INITIAL for leg in all_legs) )
+        except Exception as e:
+            raise MadGraph5Error('\nError when computing reduced quantities for current:\n%s\n'%str(self)+
+                                 'Defining flavors:\n%s\n'%str(defining_flavors)+
+                                 'routing dict:\n%s\n'%str(routing_dict)+
+                                 'Exception encountered:\n%s'%traceback.format_exc())
 
     def get_key(self, track_leg_numbers=False):
         """Return the ProcessKey associated to this current."""
@@ -888,9 +1006,48 @@
                     'singular_structure',
                     'n_loops',
                     'squared_orders',
+                    'perturbation_couplings',
                     'resolve_mother_spin_and_color' ],
                 track_leg_numbers = track_leg_numbers)
 
+    def get_initial_state_leg_map(self, leg_numbers_map):
+        """ Returns None if this current does not involve initial states, otherwise it returns 
+        pair of SubtractionLeg instances in the form of the 2-tuple 
+                             (mother_IS_leg, daughter_IS_leg)
+        Notice that beam factorization across more than one (i.e. two) beams must be 
+        factorized into two currents, so we should be guaranteed having to return at most
+        one pair of IS legs here."""
+        
+    
+        # First get the ancestors of all outermost leg numbers of this singular structures
+        ancestors = Counterterm.get_ancestor( 
+         [ leg.n for leg in self['singular_structure'].get_all_legs() ],  leg_numbers_map )
+        
+        # Filter initial state legs and make sure there is only one
+        ancestors = [leg for leg in ancestors if leg.state==SubtractionLeg.INITIAL]
+        
+        # This current has no initial-state
+        if len(ancestors) == 0:
+            return None
+        
+        if len(ancestors)!=1:
+            raise MadGraph5Error('Beam factorization current should involve'+
+                ' exactly one ancestor initial state leg (the two beam factorization must be factorized).')
+        initial_state_ancestor = ancestors[0]
+        
+        # Now find which external initial-state leg is connected to the initial
+        # state ancestor identified.
+        daughters = Counterterm.get_daughter_legs(
+                                                   initial_state_ancestor, leg_numbers_map)
+        # Filter initial state legs and make sure there is only one
+        daughters = [leg for leg in daughters if leg.state==SubtractionLeg.INITIAL]
+        if len(daughters)!=1:
+            raise CurrentImplementationError('Beam factorization current should involve'+
+                ' exactly one daughter initial state leg (the two beam factorization must be factorized).')        
+        initial_state_daughter = daughters[0]
+        
+        return (initial_state_ancestor, initial_state_daughter)
+
     def get_copy(self, copied_attributes=()):
         """Return a copy of this current with a deep-copy of its singular structure."""
 
@@ -901,22 +1058,92 @@
         
         return copied_current
 
+    def get_number_of_couplings(self):
+        """ Return a naive count of the number of couplings involved in this current."""
+        return self.get('n_loops') + self['singular_structure'].count_couplings()
+
     def __eq__(self, other):
         """Compare two currents using their ProcessKey's."""
 
         return self.get_key().key_dict == other.get_key().key_dict
 
 #=========================================================================================
+# BeamCurrent
+#=========================================================================================
+class BeamCurrent(Current):
+    """ A special class of current representing beam factorization terms 
+    (aka PDF counterterms).
+
+    If we have:
+        PDF_CT(xi) = F_+(xi) + [F] \delta(xi-1)
+    (which can also be explicitely written)
+        PDF_CT(xi) = F(xi) + {F(xi)} \delta(xi-1) + [F] \delta(xi-1)
+
+    Then this current represent the two pieces of the xi-dependent distribution F_+(xi),
+    the first piece (F(xi)) with 'distribution_type' set to 'bulk' and the latter 
+    ({F(xi)} \delta(xi-1)) with this attribute set to 'counterterm'.
+    """
+    
+    def default_setup(self):
+        """Default values for all properties specific to the BeamCurrent,
+        additional to Process.
+        """
+        super(BeamCurrent, self).default_setup()
+        self['n_loops'] = -1
+        self['squared_orders'] = {}
+        self['resolve_mother_spin_and_color'] = True
+        self['singular_structure'] = SingularStructure()
+
+        # Which piece of the plus distribution is intended to be returned
+        # This can be 'bulk' or 'counterterm'.
+        self['distribution_type'] = 'bulk'
+        
+        # List of the pdgs of the particles making up this beam
+        self['beam_PDGs'] = []
+        # Name of the beam factorization type that should be employed
+        self['beam_type'] = None
+
+    def does_require_correlated_beam_convolution(self):
+        """ Checks whether this integrated counterterm requires a host contribution featuring
+        correlated convolution of the beams (i.e. 'BS', 'VS', etc... contributions).
+        This is the case only for integrated counterterms with disjoint structures having each
+        at least one soft *BeamCurrents* that originated from a colorful subtraction currents 
+        scheme with mappings such as ppToOneWalker that recoils the soft momentum equally 
+        against both initial-state beams."""
+
+        # Return true only if the singular structures of this current are pure soft.
+        return self['singular_structure'].does_require_correlated_beam_convolution()
+
+    def get_key(self, track_leg_numbers=False):
+        """Return the ProcessKey associated to this current."""
+
+        from madgraph.core.accessors import ProcessKey
+        return ProcessKey(
+                self,
+                allowed_attributes = [
+                    'singular_structure',
+                    'n_loops',
+                    'squared_orders',
+                    'perturbation_couplings',
+                    'resolve_mother_spin_and_color',
+                    'beam_PDGs',
+                    'beam_type',
+                    'distribution_type' ],
+                track_leg_numbers = track_leg_numbers)
+    
+    def __str__(self, **opts):
+        base_str = Current.__str__(self, **opts)
+        # Denote "bulk" contributions embedded within colons ':%s:'
+        if self['distribution_type']=='bulk':
+            return ':%s:'%base_str
+        return base_str
+
+#=========================================================================================
 # Integrated current
 #=========================================================================================
-
 class IntegratedCurrent(Current):
     """Class for integrated currents."""
 
-    # TODO
-    # For now, it behaves exactly as a local 4D current,
-    # but it is conceptually different.
-
     def __str__(
         self,
         print_n=True, print_pdg=True, print_state=True,
@@ -928,7 +1155,50 @@
             print_n=print_n, print_pdg=print_pdg, print_state=print_state,
             print_loops=print_loops, print_orders=print_orders
         )
-        return '[integrated] %s' % res
+        # Denote integrated contributions embedded within '[%s]'
+        return '[%s]' % res
+
+class IntegratedBeamCurrent(IntegratedCurrent):
+    """ A special class of current representing the integrated beam factorization terms 
+
+    If we have:
+        PDF_CT(xi) = F_+(xi) + [F] \delta(xi-1)
+    (which can also be explicitely written)
+        PDF_CT(xi) = F(xi) + :F(xi): \delta(xi-1) + [F] \delta(xi-1)
+
+    Then this current returns the end-point of the distribution: [F] \delta(xi-1)
+    """
+
+    def default_setup(self):
+        """Default values for all properties specific to the BeamCurrent,
+        additional to Process.
+        """
+        super(IntegratedBeamCurrent, self).default_setup()
+        self['n_loops'] = 0
+        self['squared_orders'] = {}
+        self['resolve_mother_spin_and_color'] = True
+        self['singular_structure'] = SingularStructure()
+
+        # List of the pdgs of the particles making up this beam
+        self['beam_PDGs'] = []
+        # Name of the beam factorization type that should be employed
+        self['beam_type'] = None
+
+    def get_key(self, track_leg_numbers=False):
+        """Return the ProcessKey associated to this current."""
+
+        from madgraph.core.accessors import ProcessKey
+        return ProcessKey(
+                self,
+                allowed_attributes = [
+                    'singular_structure',
+                    'n_loops',
+                    'squared_orders',
+                    'perturbation_couplings',
+                    'resolve_mother_spin_and_color',
+                    'beam_PDGs',
+                    'beam_type' ],
+                track_leg_numbers = track_leg_numbers)
 
 #=========================================================================================
 # CountertermNode
@@ -999,6 +1269,14 @@
         for node in self.nodes:
             total_unresolved += node.count_unresolved()
         return total_unresolved
+    
+    def get_number_of_couplings(self):
+        """Count the number of couplings expected in the currents of this counterterm node."""
+        
+        total_couplings = self.current.get_number_of_couplings()
+        for node in self.nodes:
+            total_couplings += node.get_number_of_couplings()
+        return total_couplings
 
     def find_leg(self, number):
         """Find the SubtractionLeg with number specified."""
@@ -1027,13 +1305,12 @@
         add entries corresponding to the subtraction legs in this node and using the function
         get_parent_PDGs from the subtraction module."""
 
+        for node in self.nodes:
+            node.get_reduced_flavors(defining_flavors, IR_subtraction_module, routing_dict)  
+
         self.current.get_reduced_flavors(
                                     defining_flavors, IR_subtraction_module, routing_dict)
 
-        for node in self.nodes:
-            node.current.get_reduced_flavors(
-                                    defining_flavors, IR_subtraction_module, routing_dict)  
-
     def recursive_singular_structure(self, intermediate_leg_ns):
         """Reconstruct recursively the singular structure of this CountertermNode."""
 
@@ -1048,11 +1325,11 @@
         ])
         return type(structure)(legs=legs, substructures=substructures)
 
-    def split_loops(self, n_loops):
+    def distribute_loops(self, n_loops):
         """Split a counterterm in several ones according to the individual
         loop orders of its currents.
         """
-
+        
         assert isinstance(n_loops, int)
 
         # Generate all combinations of nodes with total number of loops
@@ -1067,7 +1344,7 @@
             # at any loop number between 0 and n_loops
             first_without_loops = self.nodes[0]
             for loop_n in range(n_loops + 1):
-                for first_with_loops in first_without_loops.split_loops(loop_n):
+                for first_with_loops in first_without_loops.distribute_loops(loop_n):
                     node_combinations += [[first_with_loops, ], ]
             # Add the other nodes one by one recursively,
             # taking care of never exceeding n_loops
@@ -1081,7 +1358,7 @@
                         # Distribute the order of this node
                         # in all possible ways within the current itself
                         # (recursion step)
-                        for ith_node in self.nodes[i_current].split_loops(new_loop_n):
+                        for ith_node in self.nodes[i_current].distribute_loops(new_loop_n):
                             new_node_combinations.append(
                                 combination + [ith_node, ] )
                 node_combinations = new_node_combinations
@@ -1100,7 +1377,7 @@
 
         return result
 
-    def split_orders(self, target_squared_orders):
+    def distribute_orders(self, target_squared_orders):
         """Split a counterterm in several ones according to the individual
         coupling orders of its currents.
         """
@@ -1110,17 +1387,12 @@
         # Eventually, it should create new instances, set orders there, and return the
         # new list of counterterm copies with orders
         # TODO actually create a new counterterm with orders set
-
         for node in self.nodes:
-            node.split_orders(target_squared_orders)
+            node.distribute_orders(target_squared_orders)
 
         if isinstance(self.current, Current):
             self.current.set(
-                'squared_orders',
-                { 'QCD':
-                        (self.current.get('n_loops') * 2 +
-                         self.current[
-                             'singular_structure'].count_unresolved() * 2) } )
+                'squared_orders', { 'QCD': self.current.get_number_of_couplings() * 2 } )
         else:
             # Here the squared order reduced process should be changed accordingly
             # and decreased for each squared order "sucked up" by the currents
@@ -1178,10 +1450,21 @@
         except KeyError:
             self.prefactor = self.get_prefactor(**opts)
 
+        # Both attributes below will be set via a call to `set_reduced_flavors_map`
+        # during the generation of the counterterms
+        # The dictionary below takes the form:
+        #   { ( (reolved_initial_state_PDGs,), (resolved_final_state_PDGs, ) ) :
+        #               ( (resolved_initial_state_PDGs,), (resolved_final_state_PDGs, ) ) }
+        # In the future, the values can be upgraded to list of mapped reduced flavors.
+        # But for now there can only be one at a time.
+        # Also, an additional 'None' key is added which points to the reduced flavor
+        # configuration of the defining process.
+        self.reduced_flavors_map     = None
+
     @classmethod
     def get_ancestor(cls, particles, momentum_dict):
         """Recursively explore the momentum dictionary to find an ancestor
-        of the given particle set.
+        of the given particle set (i.e. these are leg numbers here).
         """
     
         try:
@@ -1204,13 +1487,18 @@
                     "in this momentum routing dictionary:\n" + str(momentum_dict) )
 
     def __str__(self, print_n=True, print_pdg=False, print_state=False):
-
         return self.reconstruct_complete_singular_structure().__str__(
             print_n=print_n, print_pdg=print_pdg, print_state=print_state )
 
     def nice_string(self, lead="    ", tab="    "):
+        
+        CT_type = ''
+        if type(self) == Counterterm:
+            CT_type = '[local] '
+        elif type(self) == IntegratedCounterterm:
+            CT_type = '[integrated] '
 
-        tmp_str  = lead + self.process.nice_string(0, True, False)
+        tmp_str  = lead + CT_type + self.process.nice_string(0, True, False)
         tmp_str += " ("
         tmp_str += " ".join(
             str(leg['number'])
@@ -1273,23 +1561,32 @@
         
         return None
 
-    def get_daughter_pdgs(self, leg_number, state):
-        """Walk down the tree of currents to find the pdgs 'attached'
-        to a given leg number of the reduced process.
-        """
+    @classmethod
+    def get_daughter_legs(cls, leg_number, momentum_dict):
+        """ Walk down the momentum dictionary tree to find the outter-most daughter legs
+        coming from the leg_number specified."""
         
         external_leg_numbers = []
         intermediate_leg_number = [leg_number, ]
         while len(intermediate_leg_number) > 0:
             next_leg = intermediate_leg_number.pop(0)
-            # Check if this leg is final
-            daughters = self.momenta_dict[next_leg]
+            # Check if this leg is external
+            daughters = momentum_dict[next_leg]
             if  daughters == frozenset([next_leg,]):
                 external_leg_numbers.append(next_leg)
             else:
                 intermediate_leg_number = list(daughters) + intermediate_leg_number
         
-        # Now that we have the external leg numbers, find the corresponding pdg
+        return external_leg_numbers
+
+    def get_daughter_pdgs(self, leg_number, state):
+        """Walk down the tree of currents to find the pdgs 'attached'
+        to a given leg number of the reduced process.
+        """
+        
+        external_leg_numbers = self.get_daughter_legs(leg_number, self.momenta_dict)
+        
+        # Now that we have the external leg numbers, find the corresponding PDG
         pdgs = []
         for n in external_leg_numbers:
             leg = self.find_leg(n)
@@ -1394,6 +1691,14 @@
             total_unresolved += node.count_unresolved()
         return total_unresolved
 
+    def get_number_of_couplings(self):
+        """Count the number of couplings expected in the currents of this counterterm."""
+        
+        total_couplings = 0
+        for node in self.nodes:
+            total_couplings += node.get_number_of_couplings()
+        return total_couplings
+
     def get_all_currents(self):
         """Return a list of all currents involved in this counterterm."""
 
@@ -1402,36 +1707,90 @@
             currents += node.get_all_currents()
         return currents
 
-    def get_reduced_flavors(self, defining_flavors=None, IR_subtraction=None):
-        """Given the defining flavors corresponding to the resolved process (as a dictionary),
-         return a *list* of flavors corresponding to the flavor assignment of the reduced process 
-         given the defining flavors.
-        """
+    def get_beam_currents(self):
+        """ Returns a list of dictionaries of the form
+              {   'beam_one' : beam_current_1 # None for no convolution
+                  'beam_two' : beam_current_2 # None for no convolution
+              } """
+        
+        # Copy the beam currents from the reduced process to avoid border effects
+        # when further post-processing it below.
+        beam_currents = [dict(bf) for bf in self.current['beam_factorization']]
+
+        # First fetch all currents making up this counterterm
+        for current in self.get_all_currents():        
+            if type(current) not in [BeamCurrent]:
+                continue
+
+            # Now get all legs of its singular structure
+            all_legs = current['singular_structure'].get_all_legs()
+    
+            if any(l.n==1 for l in all_legs):
+                for bcs in beam_currents:
+                    if bcs['beam_one'] is not None:
+                        raise MadGraph5Error('The beam factorization currents from the reduced'+
+                ' process and the ones in the counterterms must never apply to the same beam (#1).')
+                    bcs['beam_one'] = current
+            if any(l.n==2 for l in all_legs):
+                for bcs in beam_currents:
+                    if bcs['beam_two'] is not None:
+                        raise MadGraph5Error('The beam factorization currents from the reduced'+
+                ' process and the ones in the counterterms must never apply to the same beam (#2).')
+                    bcs['beam_two'] = current
+    
+        return beam_currents
+
+    def set_reduced_flavors_map(self, defining_process, mapped_processes, IR_subtraction):
+        """Given the defining and mapped processes subject to this counterterm, 
+        this function sets the instance attribute 'reduced_flavor_map' which can be 
+        used later to easily access the corresponding reduced flavor."""
+        
         if self.count_unresolved() > 1 and any(
             any(substruct.name()=='S' for substruct in current['singular_structure'].substructures)
-                                                   for current in self.get_all_currents()):
+                for current in self.get_all_currents() if isinstance(current,(
+                                    IntegratedBeamCurrent,IntegratedCurrent,BeamCurrent))):
             raise MadGraph5Error("The function 'get_reduced_flavors' of the class Counterterm"+
-                " does not yet support NNLO soft structures such as: %s"%str(self))
-
-        # Now construct the reduced flavors list
-        if defining_flavors is None or IR_subtraction is None:
-            # If no defining flavors are specified then simply use the flavors already
-            # filled in the reduced process.
-            reduced_flavors = self.process.get_cached_initial_final_pdgs()
-        else:            
-            number_to_flavors_map = dict(defining_flavors)
+                " does not yet support integrated NNLO soft structures such as: %s"%str(self))
+
+        self.reduced_flavors_map = {}
+        for i_proc, process in enumerate([defining_process,]+mapped_processes):
+            number_to_flavors_map = { l['number'] : l['id'] for l in process.get('legs') }
+
             # Update the number_to_flavors_map dictionary by walking through the nodes
             for node in self.nodes:
                 node.get_reduced_flavors(
                                  number_to_flavors_map, IR_subtraction, self.momenta_dict)
             reduced_process_leg_numbers = self.process.get_cached_initial_final_numbers()
-
             reduced_flavors = (
                 tuple( number_to_flavors_map[n] for n in reduced_process_leg_numbers[0] ),
                 tuple( number_to_flavors_map[n] for n in reduced_process_leg_numbers[1] )
             )
-            
-        return reduced_flavors
+            self.reduced_flavors_map[process.get_cached_initial_final_pdgs()] = reduced_flavors
+            # Also add a None key for the default reduced flavors corresponding to the 
+            # resolved ones of the defining process.
+            if i_proc==0:
+                self.reduced_flavors_map[None] = reduced_flavors                
+
+    def get_reduced_flavors(self, defining_flavors=None):
+        """Given the defining flavors as a tuple (initial_state_PDGs, final_state_PDGs),
+         returns in the same format the flavors corresponding to the flavor assignment of
+         the reduced process  given the defining flavors.
+         If the defining defining_flavors is set to None, then the default reduced_flavors
+         will be returned, corresponding to the defining process that yielded this counterterm.
+        """
+
+        if self.reduced_flavors_map is None:
+            raise MadGraph5Error("The function 'set_reduced_flavors_map' of the counterterm"+
+                            " should be called before the function 'get_reduced_flavors'.")
+
+        try:
+            return self.reduced_flavors_map[defining_flavors]
+        except KeyError:
+            raise MadGraph5Error(
+                "The reduced flavor configuration for the resolved configuration "+
+                "'%s' was not found in the"%str(defining_flavors)+
+                " reduced flavors map:\n%s\nof this counterterm:\n%s"%(
+                                         str(self.reduced_flavors_map),self.nice_string()))
 
     def get_reduced_quantities(self, reduced_PS_dict, defining_flavors=None):
         """Given the PS *dictionary* providing the reduced kinematics returned
@@ -1462,25 +1821,134 @@
 class IntegratedCounterterm(Counterterm):
     """A class for the integrated counterterm."""
 
+    def __str__(self, print_n=True, print_pdg=False, print_state=False):
+        """ A nice short string for the structure of this integrated CT."""
+        
+        reconstructed_ss = self.reconstruct_complete_singular_structure()
+        
+        # We deconstruct the reconstructed ss so as to be able to render the information about
+        # wheter beam currents are local counterterms or bulk ones.
+        res = []
+        for i, node in enumerate(self.nodes):
+            if isinstance(node.current, BeamCurrent) and node.current['distribution_type']=='bulk':
+                template = ':%s:,'
+            else:
+                template = '%s,'
+            res.append(template%(reconstructed_ss.substructures[i].__str__(
+                                 print_n=print_n, print_pdg=print_pdg, print_state=print_state )))
+
+        return '[%s]'%(''.join(res))
+
     def get_integrated_current(self):
         """ This function only makes sense in the current context where we don't allow
         the factorization of the integrated counterterm into several integrated currents."""
         
         # For now we only support a basic integrated counterterm which is not broken
         # down into subcurrents but contains a single CountertermNode with a single
-        # current in it that contains the whole singular structure describing this
-        # integrated counterterm.
-        if len(self.nodes) != 1 or len(self.nodes[0].nodes) != 0:
+        # current *or* a BeamCurrent and no CountermNode
+        if len(self.nodes) == 1 and len(self.nodes[0].nodes) == 0:
+            assert isinstance(self.nodes[0].current, (IntegratedCurrent, BeamCurrent, IntegratedBeamCurrent))
+            return self.nodes[0].current
+        else:
             raise MadGraph5Error(
-                """For now, MadEvent7 only support simple integrated
+                """For now, MadNkLO only support simple integrated
                 counterterms that consists of single current encompassing the full
-                singular structure that must be analytically integrated over.""")
-        assert(isinstance(self.nodes[0].current, IntegratedCurrent))
-
-        return self.nodes[0].current
+                singular structure that must be analytically integrated over, not:\n%s"""%str(self))
+
+
+    def get_beam_currents(self):
+        """ Returns a list of dictionaries of the form
+              {   'beam_one' : beam_current_1 # None for no convolution
+                  'beam_two' : beam_current_2 # None for no convolution
+      <optional>  'correlated_convolution' : beam_current # None for not correlated convolution
+              }
+           Note that whenever 'correlated_convolution' is not set to None, the other beam_currents
+           will be None.
+        """
+        
+        beam_currents = [dict(bf) for bf in self.current['beam_factorization']]
+
+        # First fetch all integrated currents making up this integrated counterterm
+        # In the current implementation, there can only be one for now.
+        integrated_current = self.get_integrated_current()
+        
+        if type(integrated_current) not in [BeamCurrent, IntegratedBeamCurrent]:
+            return beam_currents
+        
+        if integrated_current['singular_structure'].does_require_correlated_beam_convolution():
+            for bcs in beam_currents:
+                if bcs['beam_one'] is not None or bcs['beam_two'] is not None:
+                    raise MadGraph5Error('The beam factorization currents from the reduced'+
+        ' process should all be None if the integrated counterterm necessitates a correlated'+
+        ' beam convolution.')
+                bcs['correlated_convolution'] = integrated_current
+            return beam_currents
+
+        # Now get all legs of its singular structure
+        all_legs = integrated_current['singular_structure'].get_all_legs()
+        
+        if any(l.n==1 for l in all_legs):
+            for bcs in beam_currents:
+                if bcs['beam_one'] is not None:
+                    raise MadGraph5Error('The beam factorization currents from the reduced'+
+            ' process and the ones in the integrated CT must never apply to the same beam (#1).')
+                bcs['beam_one'] = integrated_current
+        if any(l.n==2 for l in all_legs):
+            for bcs in beam_currents:
+                if bcs['beam_two'] is not None:
+                    raise MadGraph5Error('The beam factorization currents from the reduced'+
+            ' process and the ones in the integrated CT must never apply to the same beam (#2).')
+                bcs['beam_two'] = integrated_current
+
+        return beam_currents
+
+    def does_require_correlated_beam_convolution(self):
+        """ Checks whether this integrated counterterm requires a host contribution featuring
+        correlated convolution of the beams (i.e. 'BS', 'VS', etc... contributions).
+        This is the case only for integrated counterterms with pure-soft *BeamCurrent* that originated
+        from a colorful subtraction currents scheme with mappings such as ppToOneWalker that
+        recoils the soft momentum equally against both initial-state beams."""
+
+        # Fetch all integrated currents making up this integrated counterterm
+        # In the current implementation, there can only be one for now.
+        integrated_current = self.get_integrated_current()
+        return integrated_current.does_require_correlated_beam_convolution()
+
+    def get_necessary_beam_convolutions(self):
+        """ Returns a set of beam names ('beam_one' or 'beam_two') that must be active
+        in the contribution that will host this counterterm"""
+        
+        necessary_beam_convolutions = set([])
+        
+        # First analyze the reduced process
+        for bft in self.process['beam_factorization']:
+            for beam_name, beam_current in bft.items():
+                if not beam_current is None:
+                    necessary_beam_convolutions.add(beam_name)
+            
+        # Then fetch all integrated currents making up this integrated counterterm
+        # In the current implementation, there can only be one for now.
+        integrated_current = self.get_integrated_current()
+
+        if type(integrated_current) not in [BeamCurrent,]:
+            return necessary_beam_convolutions
+        
+        # Check if this integrated counterterm necessitate a soft-recoil that demands
+        # a correlated convolution of both beams
+        if integrated_current.does_require_correlated_beam_convolution():
+            return set(['beam_one','beam_two'])
+        
+        # Now get all legs of its singular structure
+        all_legs = integrated_current['singular_structure'].get_all_legs()
+        
+        if any(l.n==1 for l in all_legs):
+            necessary_beam_convolutions.add('beam_one')
+        if any(l.n==2 for l in all_legs):
+            necessary_beam_convolutions.add('beam_two')
+        return necessary_beam_convolutions
 
     def get_reduced_kinematics(self, input_reduced_PS):
-        """Given the PS *dictionary* providing the reduced kinematics corresponding to
+        """Given the PS specifying the reduced kinematics corresponding to
         the current PS point thrown at the virtual contribution, return a *dictionary* of 
         momenta corresponding to the reduced kinematics."""
 
@@ -1510,20 +1978,128 @@
 
         return reduced_PS
 
+    def n_loops_in_host_contribution(self):
+        """ Returns the number of loops that the contribution hosting this integrated counterterm
+        is supposed to have. Factorization contributions render this quantity more subtle than
+        the naive n_loops() call to self."""
+        
+        # First account for all the loops of the building blocks
+        host_contrib_n_loops  = self.n_loops()
+        host_contrib_n_loops += self.process.get_n_loops_in_beam_factorization()
+        
+        # We must then add one loop for each unresolved parton of this integrated CT.
+        host_contrib_n_loops += self.count_unresolved()
+        
+        
+        # We must now adjust the number of loops expected in the host contributions because
+        # some integrated counterterms don't yield loops in the host contributions, but instead
+        # beam factorisation convolutions.
+        beam_numbers_for_which_n_loops_is_already_decremented = []
+        for current in self.get_all_currents():
+            if not type(current)==BeamCurrent:
+                continue
+            # All the soft structures appearing as 'BeamCurrents' come from the colorful currents scheme
+            # when the a mapping such as ppToOneWalker is used, which recoils equally against both beams.
+            # In this case, the corresponding integrated soft CT is placed in dedicated contributions
+            # such as 'BS', 'VS' etc... which are defined without an increment on n_loops corresponding
+            # to the unresolved soft particle, so it must be decremented here
+            for c in current['singular_structure'].decompose():
+                if c.name()=='S':
+                    host_contrib_n_loops -= 1
+
+            if not current.does_require_correlated_beam_convolution():
+                # Finally the ISR beamCurrent contributions do have a non-zero number of unresolved contributions
+                # which increased the loop count in the line above. However, given that those must be placed
+                # together with the beam factorization currents, we should not account for these loops.
+                beam_numbers_in_currents = set(sum([[l.n for l in ss.get_all_legs() if l.n in [1,2]] 
+                        for ss in current['singular_structure'].substructures if ss.name()=='C'],[]))
+                for beam_number in beam_numbers_in_currents:
+                    if beam_number not in beam_numbers_for_which_n_loops_is_already_decremented:
+                        beam_numbers_for_which_n_loops_is_already_decremented.append(beam_number)
+                        host_contrib_n_loops -= 1
+
+        # Similarly, the *integrated* beamCurrent F terms have zero number of unresolved legs
+        # but they must be placed together with the virtuals which have a loop count. We must
+        # therefore increase the loop count by *at most* one for each of the two beams that has
+        # at least one integratedBeamCurrent F structure attached
+        beam_numbers_for_which_n_loops_is_already_incremented = []
+        for current in self.get_all_currents():
+            if not type(current)==IntegratedBeamCurrent:
+                continue
+            
+            # We don't need to look recursively inside the singular structures since the beam
+            # factorization ones are supposed to be at the top level since they factorize
+            beam_numbers_in_currents = set(sum([[l.n for l in ss.legs] for ss in 
+                current['singular_structure'].substructures if ss.name()=='F'],[]))
+            for beam_number in beam_numbers_in_currents:
+                assert (beam_number in [1,2]), "Inconsistent beam number found."
+                if beam_number not in beam_numbers_for_which_n_loops_is_already_incremented:
+                    beam_numbers_for_which_n_loops_is_already_incremented.append(beam_number)
+                    host_contrib_n_loops += 1
+                
+        return host_contrib_n_loops
+
+    def get_n_loops_in_beam_factorization_of_host_contribution(self):
+        """ Return the effective number of loops which will be part of the beam factorization
+        of the host contribution."""
+        
+        # Start by counting the number of loops in the beam factorization of the reduced
+        # process factorized
+        n_loops = self.process.get_n_loops_in_beam_factorization()
+
+        # Then add those of the xi-dependent integrated ISR counterterm and local F ones.
+        beam_numbers_in_currents = set([])
+        soft_beam_factorisation_n_loop = 0
+        for current in self.get_all_currents():
+            if not type(current)==BeamCurrent:
+                continue
+            # For the integrated soft counterterms that go in the BS, VS etc... contributions,
+            # we must only subtract exactly one as this correspond to exactly one convolution,
+            # irrespectively of the number of unresolved particles of the counterterm.
+            if self.does_require_correlated_beam_convolution():
+                soft_beam_factorisation_n_loop = 1
+            else:   
+                # We don't need to look recursively inside the singular structures since the beam
+                # factorization ones are supposed to be at the top level since they factorize
+                beam_numbers_in_current = set(sum([[l.n for l in ss.legs if l.n in [1,2]] for ss in 
+                                  current['singular_structure'].substructures if ss.name()=='C'],[]))
+                beam_numbers_in_currents |= beam_numbers_in_current
+            n_loops += current['n_loops'] + current.count_unresolved()
+        # We must remove one loop per beam number encountered in integrated ISR collinear CT
+        # since the leading beam factorization counterterm has zero loop but count_unresolved()
+        # counts for one there.
+        n_loops -= (len(beam_numbers_in_currents) + soft_beam_factorisation_n_loop)
+        
+        return n_loops
+
 #=========================================================================================
 # IRSubtraction
 #=========================================================================================
 
 class IRSubtraction(object):
 
-    def __init__(self, model, n_unresolved, coupling_types=('QCD', )):
+    _allowed_model_names = ['sm', 'loop_sm', 'loopsm', 
+                            'simple_qcd','loop_qcd_qed_sm','hc_nlo_x0_ufo']
+
+    def __init__(self, model, n_unresolved, coupling_types=('QCD', ), beam_types=(None,None),
+                 currents_scheme = 'colorful', mappings_scheme = 'LorentzNLO'):
         """Initialize a IR subtractions for a given model,
         correction order and type.
         """
-        
-        self.model = model
-        self.n_unresolved = n_unresolved
+        self.model          = model
+        self.n_unresolved   = n_unresolved
         self.coupling_types = coupling_types
+        self.beam_types     = beam_types
+        
+        self.currents_scheme = currents_scheme
+        self.mappings_scheme = mappings_scheme
+        # Decide is soft recoil against initial states
+        if self.currents_scheme in _currents_schemes_requiring_soft_beam_factorization and \
+           self.mappings_scheme in _mappings_schemes_requiring_soft_beam_factorization:
+            self.soft_do_recoil_against_initial_states = True
+        else:
+            self.soft_do_recoil_against_initial_states = False
+            
         # Map perturbed coupling orders to the corresponding interactions and particles.
         # The entries of the dictionary are
         # 'interactions', 'pert_particles' and 'soft_particles'
@@ -1561,7 +2137,8 @@
 
         if not any(
             self.model.get('name').lower().startswith(name)
-            for name in ['sm', 'loop_sm', 'loopsm', 'simple_qcd'] ):
+            for name in self._allowed_model_names
+        ):
             raise InvalidCmd(
                 "parent_PDGs_from_PDGs is implemented for SM only, "
                 "not in model %s." % self.model.get('name') )
@@ -1607,8 +2184,10 @@
                 pdgs.append(cur)
             else:
                 pdgs.append(self.model.get_particle(cur).get_anti_pdg_code())
+
         # Get the parent PDG of this leg set
         parent_PDGs = self.parent_PDGs_from_PDGs(pdgs)
+
         # Cross back if the leg set referred to initial-state
         final_PDGs = []
         for parent_PDG in parent_PDGs:
@@ -1618,6 +2197,7 @@
                 )
             else:
                 final_PDGs.append(parent_PDG)
+
         return final_PDGs
         
     def can_become_soft(self, legs):
@@ -1683,6 +2263,18 @@
                     coll_set = (coll_initial_leg, ) + coll_initial_set
                     if self.can_become_collinear(coll_set):
                         elementary_operator_list.append(CollOperator(*coll_set))
+                        
+        # Finally add an elementary BeamOperator if the process acted on has factorization
+        # Beam currents attached to it.
+        beam_names_added = []
+        for pbft in process.get('beam_factorization'):
+            for beam_name in ['beam_one','beam_two']:
+                if beam_name not in beam_names_added and (not pbft[beam_name] is None):
+                    beam_names_added.append(beam_name)
+                    elementary_operator_list.append(
+                             BeamOperator(*pbft[beam_name]['singular_structure'].get_all_legs()))
+            
+        
         return SingularOperatorList(elementary_operator_list)
 
 
@@ -1767,8 +2359,10 @@
 
             # The squared orders of the reduced process will be set correctly later
             reduced_process = reduced_process.get_copy(
-                ['legs', 'n_loops', 'legs_with_decays']
+                ['legs', 'n_loops', 'legs_with_decays','beam_factorization']
             )
+            # We need a deep copy of the beam_factorization here.
+            reduced_process['beam_factorization'] = copy.deepcopy(reduced_process['beam_factorization'])
             # Empty legs_with_decays: it will be regenerated automatically when asked for
             reduced_process['legs_with_decays'][:] = []
             # TODO If resetting n_loops, what about orders?
@@ -1778,7 +2372,6 @@
         nodes = []
 
         # 2. Recursively look into substructures
-
         current_args = set(structure.legs)
         for substructure in structure.substructures:
             node = self.get_counterterm(
@@ -1808,12 +2401,24 @@
             # Replace soft structures with their flattened versions
             elif current_structure.name() == "S":
                 current_args.add(current_structure)
+            # Remove all non-identity beam factorization currents for that leg of 
+            # type 'bulk' to be the identity operator ('None').
+            elif current_structure.name() == "F":
+                if len(current_structure.legs)!=1 or current_structure.legs[0].n not in [1,2]:
+                    raise MadGraph5Error("Only beam factorizations attached to leg number 1 or"+
+                                      " 2 are supported for now; not %s."%str(current_structure))
+                current_args.add(current_structure)
+                beam_names = {1:'beam_one', 2:'beam_two'}
+                # Change the distribution type of all beam_factorization terms of that beam
+                # in the reduced process to be the identity (i.e. None).
+                for bft in reduced_process.get('beam_factorization'):
+                    beam_name = beam_names[current_structure.legs[0].n]
+                    if not bft[beam_name] is None:
+                        bft[beam_name]=None
+
             # Other structures need to be implemented
             else:
-                raise MadGraph5Error(
-                    "Unrecognized current of type %s" %
-                    str(type(current_structure))
-                )
+                raise MadGraph5Error("Unrecognized current of type %s" % str(type(current_structure)))
             # Add this node
             nodes.append(node)
 
@@ -1833,12 +2438,30 @@
         # 3. Else build the current and update
         #    the reduced process as well as the dictionary
         current_type = type(structure)(*current_args)
-        current = Current({
-            'singular_structure': current_type })
+        if current_type.name()=='F':
+            # The beam factorization singular structure always have a single leg with a number
+            # matching the beam number.
+            beam_number = current_type.legs[0].n
+            assert (beam_number in [1,2]), "Beam factorization structure are expected to "+\
+                                                        "have a single leg withh number in [1,2]."
+            current = BeamCurrent({
+                'beam_type'         : self.beam_types[beam_number-1][0],
+                'beam_PDGs'         : self.beam_types[beam_number-1][1],
+                'distribution_type' : 'counterterm',
+                'singular_structure': current_type })
+        else:
+            current = Current({'singular_structure': current_type })
         structure_legs = current_type.get_all_legs()
         structure_leg_ns = frozenset(leg.n for leg in structure_legs)
+        if current_type.name()=='F':
+            # We must not remove legs from the reduced process for beam factorization 
+            # contributions
+            structure_leg_ns = frozenset() 
+        else:
+            structure_leg_ns = frozenset(leg.n for leg in structure_legs)
+           
         parent = None
-        if structure.name() == "C":
+        if structure.name() in ["C"]:
             # Add entry to dictionary
             parent_index = len(momenta_dict_so_far) + 1
             momenta_dict_so_far[parent_index] = structure_leg_ns
@@ -1853,10 +2476,11 @@
         elif structure.name() == "S":
             # No parent
             pass
+        elif structure.name() == "F":
+            # No need to propagate parents for F structures
+            pass
         else:
-            raise MadGraph5Error(
-                "Building unrecognized current of type %s" %
-                str(type(structure)) )
+            raise MadGraph5Error("Building unrecognized current of type %s" %str(type(structure)) )
         # Remove legs of this structure
         legs_to_remove = []
         for leg in reduced_process['legs']:
@@ -1888,13 +2512,15 @@
         # First check that the local counterterm is singular, because if not then we
         # should of course not return any integrated counterterm.
         if not local_counterterm.is_singular():
-            return None
-
-        complete_singular_structure = local_counterterm. \
-            reconstruct_complete_singular_structure()
+            return []
+
+        complete_singular_structure = local_counterterm.reconstruct_complete_singular_structure()
+
+        # Useful to also have access to the flatten complete singular structure
+        flatten_singular_structure = complete_singular_structure.decompose()
 
         reduced_process = local_counterterm.process.get_copy(
-            ['legs', 'n_loops', 'legs_with_decays', 'squared_orders'] )
+            ['legs', 'n_loops', 'legs_with_decays', 'squared_orders', 'beam_factorization'] )
 
         # The following sums all the loop numbers in all nodes and reduced process,
         # the latter of which must then be removed
@@ -1923,17 +2549,117 @@
         #
         ########
         
-        integrated_current = IntegratedCurrent({
-            'n_loops': n_loops,
-            'squared_orders': squared_orders,
-            'resolve_mother_spin_and_color': True,
-            'singular_structure': complete_singular_structure })
-
-        return IntegratedCounterterm(
-            process=reduced_process,
-            nodes=[CountertermNode(current=integrated_current), ],
-            momenta_dict=bidict(local_counterterm.momenta_dict),
-            prefactor=-1. * local_counterterm.prefactor )
+        ########
+        #
+        # TODO handle the mixed final-initial local counterterm case. This issue is more general, however and pertains
+        # to the definition of the integrated countertem(S?) of any set of disjoint local counterterms. 
+        # For now only the case of a single initial state collinear local CT is supported.
+        #
+        ########
+        
+        # First tackle the special case of single initial-state collinear structure
+        initial_state_legs = [leg for leg in complete_singular_structure.substructures[0].legs if
+                                                       leg.state == SubtractionLeg.INITIAL]
+ 
+        # A list to hold all integrated CT generated here.
+        integrated_CTs = []
+        
+        # Check if this local counterterm has a soft component that recoils symmetrically
+        # against both beams.
+        has_soft_symmetric_ISR_recoil = (self.soft_do_recoil_against_initial_states and 
+                complete_singular_structure.does_require_correlated_beam_convolution())
+
+        # Handle the specific case of single initial-state pure collinear counterterm or 
+        # soft recoil against IS.
+        if (len(initial_state_legs) == 1 and not any(ss.name()=='S' 
+                   for ss in flatten_singular_structure)) or has_soft_symmetric_ISR_recoil:
+            
+            if has_soft_symmetric_ISR_recoil:
+                # In the case of soft recoiling against initial states, the beam_type should not matter.
+                beam_type = None
+                beam_PDGs = None
+            else:  
+                beam_number = initial_state_legs[0].n
+                beam_names = {1:'beam_one', 2:'beam_two'}
+                        
+                assert (beam_number in [1,2]), "MadNkLO only supports initial state legs with number in [1,2]."
+                beam_type = self.beam_types[beam_number-1][0] if self.beam_types[beam_number-1] else None
+                beam_PDGs = self.beam_types[beam_number-1][1] if self.beam_types[beam_number-1] else None
+
+            beam_current_options = {
+                    'beam_type'     :   beam_type,
+                    'beam_PDGs'     :   beam_PDGs,
+                    'n_loops'       :   n_loops,
+                    'squared_orders':   squared_orders,
+                    'singular_structure': complete_singular_structure
+                }
+
+            # For now there will be no additional CT node, but for NNLO there may be because
+            # the integrated ISR has to be combined with local FSR counterterms
+            additional_CT_nodes = []
+            
+            # Handle the specific case of NLO single beam factorization counterterm
+            if complete_singular_structure.substructures[0].name() == 'F':
+
+                integrated_CTs.append(IntegratedCounterterm(
+                    process     = reduced_process,
+                    nodes       = [
+                            CountertermNode(current=IntegratedBeamCurrent(beam_current_options)), 
+                        ]+additional_CT_nodes,
+                    momenta_dict= bidict(local_counterterm.momenta_dict),
+                    prefactor   = -1. * local_counterterm.prefactor ))
+
+            # Handle the specific case of ISR recoils, either because of ISR collinear or
+            # soft counterterms in colorful and ppToOneWalker mapping.
+            elif has_soft_symmetric_ISR_recoil or complete_singular_structure.substructures[0].name() == 'C':
+
+                # Now add the three pieces of the integrated ISR. Ultimately we may choose to instead
+                # generate only the first piece ('bulk') and the other two via an iterative application
+                # of the subtraction procedure. For now however, we will generate all three at once.
+                
+                # Starting with the integrated endpoint:
+                integrated_CTs.append(IntegratedCounterterm(
+                    process     = reduced_process,
+                    nodes       = [
+                            CountertermNode(current=IntegratedBeamCurrent(beam_current_options)), 
+                        ]+additional_CT_nodes,
+                    momenta_dict= bidict(local_counterterm.momenta_dict),
+                    prefactor   = -1. * local_counterterm.prefactor ))
+    
+                # Then the bulk and counterterm integrated ISR pieces, which are both xi dependent.
+                beam_current_options['distribution_type'] = 'bulk'
+                integrated_CTs.append(IntegratedCounterterm(
+                    process     = reduced_process,
+                    nodes       = [
+                            CountertermNode(current=BeamCurrent(beam_current_options)), 
+                        ]+additional_CT_nodes,
+                    momenta_dict= bidict(local_counterterm.momenta_dict),
+                    prefactor   = -1. * local_counterterm.prefactor ))
+                
+                beam_current_options['distribution_type'] = 'counterterm'
+                integrated_CTs.append(IntegratedCounterterm(
+                    process     = reduced_process,
+                    nodes       = [
+                            CountertermNode(current=BeamCurrent(beam_current_options)), 
+                        ]+additional_CT_nodes,
+                    momenta_dict= bidict(local_counterterm.momenta_dict),
+                    prefactor   = +1. * local_counterterm.prefactor ))
+            
+        else:
+            # Here is the general solution chosen for arbitrary final state local CT
+            integrated_current = IntegratedCurrent({
+                'n_loops': n_loops,
+                'squared_orders': squared_orders,
+                'resolve_mother_spin_and_color': True,
+                'singular_structure': complete_singular_structure })
+
+            integrated_CTs.append(IntegratedCounterterm(
+                process     = reduced_process,
+                nodes       = [CountertermNode(current=integrated_current), ],
+                momenta_dict= bidict(local_counterterm.momenta_dict),
+                prefactor   = -1. * local_counterterm.prefactor ))
+        
+        return integrated_CTs
 
     @staticmethod
     def get_all_currents(counterterms):
@@ -1953,7 +2679,7 @@
         ignore_integrated_counterterms=False):
         """Generate all counterterms for the corrections specified in this module
         and the process given in argument."""
-
+        
         if max_unresolved_in_elementary is None:
             max_unresolved_in_elementary = self.n_unresolved
         if max_unresolved_in_combination is None:
@@ -1962,33 +2688,34 @@
             process, max_unresolved_in_elementary)
         combinations = self.get_all_combinations(
             elementary_operators, max_unresolved_in_combination)
+
         all_counterterms = []
         all_integrated_counterterms = []
         for combination in combinations:
             template_counterterm = self.get_counterterm(combination, process)
             if not ignore_integrated_counterterms:
-                template_integrated_counterterm = \
+                template_integrated_counterterms = \
                                        self.get_integrated_counterterm(template_counterterm)
             else:
-                template_integrated_counterterm = None
+                template_integrated_counterterms = []
 
-            counterterms_with_loops = template_counterterm.split_loops(process['n_loops'])
+            counterterms_with_loops = template_counterterm.distribute_loops(process['n_loops'])
             # TODO
-            # For the time being, split_loops is given None instead of the squared orders
+            # For the time being, distribute_loops is given None instead of the squared orders
             # because they should be retrieved from the process by looking at individual
             # matrix elements
             # That is, every process has a list of possible coupling orders assignations
             # so we should loop over them
             for counterterm_with_loops in counterterms_with_loops:
-                all_counterterms.extend( counterterm_with_loops.split_orders(None) )
-            # Now also distribute the template integrated counterterm if it is not None
-            if not template_integrated_counterterm is None:
-                integrated_counterterms_with_loops = template_integrated_counterterm.split_loops(
+                all_counterterms.extend( counterterm_with_loops.distribute_orders(None) )
+
+            # Now also distribute the template integrated counterterms
+            for template_integrated_counterterm in template_integrated_counterterms:
+                integrated_counterterms_with_loops = template_integrated_counterterm.distribute_loops(
                     process['n_loops'] )
                 for integrated_counterterm_with_loops in integrated_counterterms_with_loops:
                     all_integrated_counterterms.extend(
-                        integrated_counterterm_with_loops.split_orders(None) )
-
+                        integrated_counterterm_with_loops.distribute_orders(None) )
         return all_counterterms, all_integrated_counterterms
 
 #=========================================================================================
@@ -2118,8 +2845,10 @@
         for current in currents:
             found_current_class = None
             for (dir_name, module_path, class_name, implementation_class) in all_classes:
+#                misc.sprint('Testing class %s for current: %s'%(class_name,str(current)))
                 instantiation_options = implementation_class.does_implement_this_current(
                     current, self.model )
+#                misc.sprint('Result: %s'%str(instantiation_options))
                 if instantiation_options is None:
                     continue
                 if found_current_class is not None:
@@ -2165,7 +2894,7 @@
         # (it should never be used in production)
         if currents_with_default_implementation:
             currents_str = '\n'.join(
-                ' > %s' % str(crt)
+                ' > %s (type: %s)' % (str(crt), type(crt) )
                 for crt in currents_with_default_implementation )
             msg = """No implementation was found for the following subtraction currents:
 %s

=== modified file 'madgraph/integrator/ME7_integrands.py' (properties changed: +x to -x)
--- madgraph/integrator/ME7_integrands.py	2018-08-13 20:45:37 +0000
+++ madgraph/integrator/ME7_integrands.py	2019-01-14 20:53:09 +0000
@@ -73,6 +73,7 @@
 import madgraph.integrator.phase_space_generators as phase_space_generators
 import madgraph.integrator.pyCubaIntegrator as pyCubaIntegrator
 import madgraph.integrator.vegas3_integrator as vegas3_integrator
+import madgraph.integrator.vectors as vectors
 #    import madgraph.various.histograms as histograms  # imported later to not slow down the loading of the code
 import models.check_param_card as check_param_card
 import models.model_reader as model_reader
@@ -81,8 +82,6 @@
 from madgraph.iolibs.files import ln    
 from madgraph import InvalidCmd, MadGraph5Error, MG5DIR, ReadWrite
 
-
-
 class MadEvent7Error(Exception):
     pass
 
@@ -90,6 +89,551 @@
     pass
 
 #===============================================================================
+# ME7Events
+#===============================================================================
+class ME7Event(object):
+    """ A class to store a particular configuration produced by the integrands
+    when evaluated and which has undefinite flavor configuration to begin with."""
+    
+    def __init__(self, PS_point, weights_per_flavor_configurations, 
+                 requires_mirroring           = False,
+                 host_contribution_definition = 'N/A',
+                 # This entry is mostly for debugging purposes, but it is useful to be able
+                 # to know what is the list of defining PDGs of the resolved process from 
+                 # which this counterterm originates. For this reason the counterterm_structure
+                 # is (possible a list of) tuple(s) of the form:
+                 #   (
+                 #      integrated_CT_instance, 
+                 #      list_of_all combination_of_PDGs_of_the_resolved_process_this_integrated_CT_can_comes_from
+                 #   )
+                 counterterm_structure        = None,
+                 Bjorken_xs                   = (1.0,1.0),
+                 Bjorken_x_rescalings         = (1.0,1.0),
+                 is_a_mirrored_event          = False,
+                 # The user is free to add below any additional information that 
+                 # plays no structural role in in the MadNkLO construction. Example:
+                 # test_IR_poles will store here information such as the input_mapping used
+                 # and the convolutional mask employed for nicer printouts.
+                 additional_information       = {}
+                ):
+        """Initialize this particular event with a unique PS point
+        and a dictionary of the form:
+             ( (initial_state_flavors,), (final_state_flavors,) ) : weight
+        to specify all possible flavor configurations applying to this event.
+        For PDF counterterms and integrated collinear CT counterterms, we must
+        provide the possibility of applying a rescaling to the Bjorken x that will be
+        used to compute each of the two beams. By default this rescaling is of course one."""
+        
+        self.PS_point = PS_point.to_list()
+        self.weights_per_flavor_configurations = weights_per_flavor_configurations
+        self.requires_mirroring = requires_mirroring
+        self.is_a_mirrored_event = is_a_mirrored_event
+        self.host_contribution_definition = host_contribution_definition
+        self.counterterm_structure = counterterm_structure
+        # This can either be None or a dictionary of the format:
+        # {  'beam_one' : 'ALL' or (PDG_allowed_1, PDG_allowed_2, PDG_allowed_3, ...)
+        #    'beam_two' : 'ALL' or (PDG_allowed_1, PDG_allowed_2, PDG_allowed_3, ...)
+        # }
+        self.Bjorken_x_rescalings = Bjorken_x_rescalings
+        self.Bjorken_xs = Bjorken_xs
+        self.additional_information = additional_information
+
+        # This entry specifies which flavor has been selected for this event.
+        # None implies that it has not been selected yet.
+        self.selected_flavors = None
+        
+    def get_copy(self):
+        """ Returns a shallow copy of self, except for the 
+        'weights_per_flavor_configurations' attribute. """
+        
+        return ME7Event(
+            self.PS_point, dict(self.weights_per_flavor_configurations), 
+            requires_mirroring           = self.requires_mirroring,
+            host_contribution_definition = self.host_contribution_definition,
+            counterterm_structure        = self.counterterm_structure,
+            Bjorken_x_rescalings         = self.Bjorken_x_rescalings,
+            Bjorken_xs                   = self.Bjorken_xs,
+            is_a_mirrored_event          = self.is_a_mirrored_event,
+            additional_information       = copy.deepcopy(self.additional_information)
+        )
+        
+    def set_Bjorken_rescalings(self, *args):
+        """ Assigns the specified Bjorken x rescaling to this event."""
+        self.Bjorken_x_rescalings = tuple(args)
+
+    def set_Bjorken_xs(self, *args):
+        """ Assigns the specified Bjorken x's."""
+        self.Bjorken_xs = tuple(args)
+
+    def apply_PDF_convolution(self, pdf_accessor, all_pdf, all_mu_f_squared):
+        """ Convolute the weights_per_flavor_configurations."""
+        
+        # In order for the + distributions of the PDF counterterms and integrated
+        # collinear ISR counterterms to hit the PDF only (and not the matrix elements or
+        # observables functions), a change of variable is necessary: xb_1' = xb_1 * xi1
+        # In this case, the xi1 to factor to apply is given by the Bjorken_x_rescalings
+        # factor, which must be used both for rescaling the argument of the PDF as well
+        # as bringing a factor 1/xi for the jacobian of the change of variable.
+        for flavors in self.weights_per_flavor_configurations:
+            PDFs = 1.
+            for i, flavor in enumerate(flavors[0]):
+                PDFs *= pdf_accessor(all_pdf[i], flavor, 
+                     self.Bjorken_xs[i]*self.Bjorken_x_rescalings[i], all_mu_f_squared[i])
+                PDFs *= self.Bjorken_x_rescalings[i]
+            self.weights_per_flavor_configurations[flavors] *= PDFs
+
+    def store_information(self, key, value):
+        """ This function registers a new additional information that will be attached to this
+        event."""
+        self.additional_information[key] = value
+
+    def get_information(self, key):
+        """ Tries to fetch a particular additional information registered in this event.
+        If not present, it returns None."""
+        return self.additional_information.get(key, None)
+
+    def get_total_weight(self):
+        """ Return the total weight of this event for all flavor considerations to consider."""
+
+        summed_weight = None
+        for wgt in self.weights_per_flavor_configurations.values():
+            if summed_weight is None:
+                summed_weight = copy.copy(wgt)
+            else:
+                summed_weight += wgt
+        if summed_weight is None:
+            return 0.
+        
+        return summed_weight
+
+    def generate_mirrored_event(self):
+        """ Creates a mirror event of self if necessary."""
+        if not self.requires_mirroring or self.is_empty():
+            return None
+        
+
+        # Now that the mirroring was applied to this event, it no longer needs any
+        self.requires_mirroring = False
+        # Neither does the mirrored event we instantiate below
+        mirrored_event = self.get_copy()
+        # Set the mirrored event flag "is_a_mirrored_event" appropriately
+        mirrored_event.is_a_mirrored_event = True
+
+        # First swap initial state flavors
+        mirrored_event.weights_per_flavor_configurations = {}
+        for flavor_config, wgt in self.weights_per_flavor_configurations.items():
+            swapped_flavor_config = ( (flavor_config[0][1],flavor_config[0][0]),
+                                           flavor_config[1] )        
+            mirrored_event.weights_per_flavor_configurations[swapped_flavor_config] = wgt
+        
+        # Then swap the kinematic configurations: 
+        # a) Initial-state momenta must be swapped
+        # b) Final-state momenta must have their z-component swapped (for this to be correct
+        # however, first make sure that initial momenta are indeed on the z-axis in the c.o.m frame)
+        # Normally all events generated should be in the c.o.m frame already, but it may be
+        # that in test_IR_limits, for debugging purposes only, the option 'boost_back_to_com' was
+        # set to False, in which case we must here first *temporarily boost* it back to the c.o.m.
+        if __debug__:
+            # Two initial states
+            assert(len(self.weights_per_flavor_configurations.keys()[0][0])==2)
+            sqrts = math.sqrt((self.PS_point[0]+self.PS_point[1]).square())
+            # Assert initial states along the z axis
+            assert(abs(self.PS_point[0][1]/sqrts)<1.0e-9)
+            assert(abs(self.PS_point[1][1]/sqrts)<1.0e-9)
+            assert(abs(self.PS_point[0][2]/sqrts)<1.0e-9)
+            assert(abs(self.PS_point[1][2]/sqrts)<1.0e-9)
+        # If initial states are back to back we can directly proceed with a simple swap of the
+        # z-axis, otherwise we must first boost to the c.o.m
+        PS_point_to_swap = self.PS_point.get_copy()
+        boost_vector = None
+        if abs((self.PS_point[0]+self.PS_point[1])[3]/sqrts)>1.0e-9:
+            boost_vector = (PS_point_to_swap[0]+PS_point_to_swap[1]).boostVector()
+            for vector in PS_point_to_swap:
+                vector.boost(-boost_vector)
+        # debug: now make sure the event is back to back
+        if __debug__:
+            sqrts = math.sqrt((self.PS_point[0]+self.PS_point[1]).square())
+            assert(abs((PS_point_to_swap[0]+PS_point_to_swap[1])[3]/sqrts)<1.0e-9)
+        mirrored_event.PS_point = vectors.LorentzVectorList([ PS_point_to_swap[1],PS_point_to_swap[0] ])
+        for vector in PS_point_to_swap[2:]:
+            mirrored_event.PS_point.append(vectors.LorentzVector(
+                                               [vector[0],vector[1],vector[2],-vector[3]]))
+        # And now if we had boosted the event we must now boost it back
+        if boost_vector is not None:
+            for vector in PS_point_to_swap:
+                vector.boost(boost_vector)
+        
+        # Then swap Bjorken x's and rescaling.
+        mirrored_event.Bjorken_x_rescalings = (self.Bjorken_x_rescalings[1], self.Bjorken_x_rescalings[0])
+        mirrored_event.Bjorken_xs = (self.Bjorken_xs[1], self.Bjorken_xs[0])
+        
+        return mirrored_event
+
+    def filter_flavor_configurations(self, flavor_cut, **opts):
+        """ Apply the flavor cut on this event for each flavor configuration, removing
+        all those failings. If there is None left, this returns False."""
+        
+        new_weights_per_flavor_configurations = {}
+        for flavors, weight in self.weights_per_flavor_configurations.items():
+            if flavor_cut(self.PS_point, flavors, xb_1=self.Bjorken_xs[0], 
+                                                        xb_2=self.Bjorken_xs[1], **opts):
+                new_weights_per_flavor_configurations[flavors] = weight
+        
+        self.weights_per_flavor_configurations = new_weights_per_flavor_configurations
+        
+        return (not self.is_empty())
+
+    def apply_flavor_blind_cuts(self, flavor_blind_cut, *args, **opts):
+        """ Apply the flavor-blind cut to this event, returning False if it failed."""
+        
+        return flavor_blind_cut(self.PS_point, *args, xb_1=self.Bjorken_xs[0], 
+                                                        xb_2=self.Bjorken_xs[1], **opts)
+
+    def select_a_flavor_configuration(self):
+        """ Select a particular flavor configuration for this event which may be used
+        by the flavor-sensitive observables. The selection is performed randomly with 
+        a probability proportional to the weight of each flavor configuration."""
+
+        weights_by_flavor_configurations = self.weights_per_flavor_configurations.items()
+        index_selected=0
+        abs_flavor_wgts_sum = sum(abs(value) for value in 
+                                         self.weights_per_flavor_configurations.values())
+        running_sum = abs(weights_by_flavor_configurations[index_selected][1])
+        rv_flavor = random.random()*abs_flavor_wgts_sum
+        while rv_flavor > running_sum:
+            index_selected +=1
+            running_sum += abs(weights_by_flavor_configurations[index_selected][1])
+        self.selected_flavors = weights_by_flavor_configurations[index_selected][0]
+
+    def select_epsilon_expansion_term(self, term_name):
+        """ Select a particular EpsilonExpansion term in the weights of the flavor 
+        matrix of this event. If the weights are float already, then they will be
+        set to zero unless the term_name is 'finite'."""
+
+        if not isinstance(term_name, int):
+            is_finite_specified = (term_name.lower()=='finite')
+        for flavors, weight in self.weights_per_flavor_configurations.items():
+            if isinstance(weight, float):
+                if not is_finite_specified:
+                    del self.weights_per_flavor_configurations[flavors]
+            else:
+                self.weights_per_flavor_configurations[flavors] = weight.get_term(term_name)
+
+    def convolve_flavors(self, flavor_matrix, leg_index, initial_state=True):
+        """ Convolves the flavor matrix in argument (of the form: 
+            {   start_flavor_a :  {  end_flavors_1 : wgt_x,
+                                     end_flavors_2 : wgt_y,
+                                     [...]
+                                 },
+                start_flavor_b :  {  end_flavors_1 : wgt_z,
+                                     end_flavors_2 : wgt_w,
+                                     [...]
+                                 },
+                [...]
+            }
+            where 'wgt_x' can be EpsilonExpansion instances.
+            and modifies self.weights_per_flavor_configurations with the convolved result.
+            leg index specifies which leg must be convolved.
+               (e.g. beam #1 correpsonds to leg_index=0 and initial_state = True.
+            Note that end_flavors is a tuple of flavor with identical matrix elements (such as in q -> q g, which is flavor independent)
+        """
+            
+        def substitute_flavor(flavor_config, new_flavor):
+            """ Substitute in the flavor config the convoluted flavor with the one in 
+            argument and returns the corresponding new flavor configuration as a two tuples,
+            for the initial and final states respectively."""
+            return ( 
+                tuple( (pdg if i!=leg_index else new_flavor) 
+                       for i, pdg in enumerate(flavor_config[0]) ) if initial_state else flavor_config[0],
+                tuple( (pdg if i!=leg_index else new_flavor) 
+                       for i, pdg in enumerate(flavor_config[1]) ) if (not initial_state) else flavor_config[1]
+            )
+            
+        if len(flavor_matrix)==0:
+            # shortcut in case an empty convolutional flavor matrix is specified which can happen
+            # if the beam currents are zero but don't have the attribute is_zero = True
+            self.weights_per_flavor_configurations = {}
+            return
+
+        new_flavor_configurations = {}
+        for flavor_config, wgt in self.weights_per_flavor_configurations.items():
+            start_flavor = flavor_config[0][leg_index] if initial_state else flavor_config[1][leg_index]
+            if start_flavor in flavor_matrix:
+                new_flavor_configs = []
+                for end_flavors, matrix_wgt in flavor_matrix[start_flavor].items():
+                    new_flavor_configs.extend([
+                    ( substitute_flavor(flavor_config, end_flavor), 
+                                base_objects.EpsilonExpansion(matrix_wgt)*wgt )
+                                                           for end_flavor in end_flavors ])
+
+                for new_flavor_config, new_wgt in new_flavor_configs:
+                    try:
+                        new_flavor_configurations[new_flavor_config] += new_wgt
+                    except KeyError:
+                        new_flavor_configurations[new_flavor_config] = new_wgt
+        
+        # Now assign the newly create flavor configurations
+        self.weights_per_flavor_configurations = new_flavor_configurations
+
+    def is_empty(self):
+        """
+        Check if the weight vector (weights_per_flavor_configurations) is an empty dictionary.
+        For now does not check if all weights are zero
+        :return: bool
+        """
+        return len(self.weights_per_flavor_configurations) == 0
+
+    def __add__(self, other):
+        """ overload the '+' operator."""
+        new_event = self.get_copy()
+        return new_event.__iadd__(other)
+            
+    def __iadd__(self, other):
+        """ overload the '+=' operator."""
+        if __debug__: assert(self.requires_mirroring == other.requires_mirroring)
+        if __debug__: assert(self.PS_point == other.PS_point)
+        if __debug__: assert(self.Bjorken_x_rescalings  == other.Bjorken_x_rescalings)
+        if __debug__: assert(self.Bjorken_xs  == other.Bjorken_xs)
+        for flavor_configuration, wgt in other.weights_per_flavor_configurations.items():
+            try:
+                self.weights_per_flavor_configurations[flavor_configuration] += wgt
+            except KeyError:
+                self.weights_per_flavor_configurations[flavor_configuration] = wgt
+        return self
+
+    def __mul__(self, multiplier):
+        """ overload the '*' operator with a float or an EpsilonExpansion. """
+        new_event = self.get_copy()
+        return new_event.__imul__(multiplier)
+    
+    def __imul__(self, multiplier):
+        """ overload the '*=' operator with a float or an EpsilonExpansion. """
+        if __debug__: assert( isinstance(multiplier, (float, base_objects.EpsilonExpansion)) )
+        for flavor_configuration, wgt in self.weights_per_flavor_configurations.items():
+            try:
+                self.weights_per_flavor_configurations[flavor_configuration] = multiplier*wgt
+            except KeyError:
+                self.weights_per_flavor_configurations[flavor_configuration] = multiplier*wgt
+        return self
+
+    def __str__(self):
+        return self.nice_string()
+
+    def nice_string(self):
+        """ Returns a nice and short string representation of this event."""
+        BLUE = misc.bcolors.BLUE
+        GREEN = misc.bcolors.GREEN
+        ENDC = misc.bcolors.ENDC
+        res = []
+        if self.counterterm_structure is None:
+            res.append("%s%svent from %s contribution '%s'%s%s:"%(
+                GREEN,
+                'E' if not self.is_a_mirrored_event else 'Mirrored e',
+                self.host_contribution_definition.correction_order,
+                self.host_contribution_definition.short_name(),
+                ' (requires mirroring)' if self.requires_mirroring else '',
+                ENDC))
+        else:
+            if isinstance(self.counterterm_structure, list):
+                all_ct_strucs = self.counterterm_structure
+            else:
+                all_ct_strucs = [self.counterterm_structure,]
+            
+            counterterm_structure_str_elems = []
+            for ct_struct, all_resolved_PDGs in all_ct_strucs:
+                # We show only one representative resolved_PDGs structure: the first one.
+                if isinstance(all_resolved_PDGs, dict):
+                    representative_resolved_PDGs = all_resolved_PDGs.keys()[0]
+                else:
+                    representative_resolved_PDGs = all_resolved_PDGs[0]                    
+                counterterm_structure_str_elems.append('%s@%s'%(str(ct_struct), 
+                                        str(representative_resolved_PDGs).replace(' ','')))
+            
+            if len(all_ct_strucs)>1:
+                counterterm_structure_str = '( %s )'%(' + '.join(counterterm_structure_str_elems))
+            else:
+                counterterm_structure_str = counterterm_structure_str_elems[0]
+
+            res.append('%s%sounterterm event from limit%s %s of %s contribution %s%s'%(
+                GREEN,
+                'C' if not self.is_a_mirrored_event else 'Mirrored c',
+                's' if isinstance(self.counterterm_structure, list) and len(self.counterterm_structure)>1 else '',
+                counterterm_structure_str,
+                self.host_contribution_definition.correction_order,
+                self.host_contribution_definition.short_name(),
+                ENDC))
+            
+        res.append('%s  Kinematic configuration:%s'%(BLUE, ENDC))
+        res.extend('     %s'%line for line in str(self.PS_point).split('\n'))
+        if not any(bs is None for bs in self.Bjorken_x_rescalings):
+            res.append("     Bjorken x's scalings: %s"%(', '.join('%.16e'%bs for bs in self.Bjorken_x_rescalings)))
+        if not any(bs is None for bs in self.Bjorken_xs):
+            res.append("     Bjorken x's         : %s"%(', '.join('%.16e'%bs for bs in self.Bjorken_xs)))
+        res.append('%s  Flavors configurations:%s'%(BLUE, ENDC))
+        for flavors in sorted(self.weights_per_flavor_configurations.keys()):
+            wgt = self.weights_per_flavor_configurations[flavors]
+            res.append('%s  |  %s'%('%-25s'%('     %s -> %s'%(
+                ' '.join('%s'%pdg for pdg in flavors[0]),
+                ' '.join('%s'%pdg for pdg in flavors[1]))), 
+                wgt.__str__(format='.16e') if isinstance(wgt,base_objects.EpsilonExpansion) 
+                                                                         else '%.16e'%wgt ))
+        if not self.selected_flavors is None:
+            res.append('%s  Selected flavors configuration : %s%s -> %s'%(BLUE,  ENDC,
+             ' '.join('%s'%pdg for pdg in self.selected_flavors[0]), 
+             ' '.join('%s'%pdg for pdg in self.selected_flavors[1]) ))
+        
+        if len(self.additional_information)>0:
+            res.append('%s  Additional information:%s'%(BLUE,  ENDC))            
+            for key, value in self.additional_information.items():
+                res.append('    %s : %s'%(str(key), str(value)))
+
+        return '\n'.join(res)
+
+    def counterterm_structure_short_string(self):
+        """ Return a short string specifying the counterterm structure(s) of this event."""
+
+        if self.counterterm_structure is None:
+            if not self.is_a_mirrored_event:
+                return 'Event'
+            else:
+                return '<->Event'
+
+        event_str = str(self.counterterm_structure)
+        if isinstance(self.counterterm_structure, list):
+            event_ct_structs = self.counterterm_structure
+        else:
+            event_ct_structs = [self.counterterm_structure,]
+        
+        event_str_elems = []
+        for ct_struct, all_resolved_PDGs in event_ct_structs:
+            # We show only one representative resolved_PDGs structure: the first one.
+            if isinstance(all_resolved_PDGs, dict):
+                representative_resolved_PDGs = all_resolved_PDGs.keys()[0]
+            else:
+                representative_resolved_PDGs = all_resolved_PDGs[0]
+            str_ct_struct = str(ct_struct)
+            # Clean-up of the embedding overall parenthesis for the short string
+            if str_ct_struct.startswith("("):
+                str_ct_struct = str_ct_struct[1:]
+            if str_ct_struct.endswith(",)"):
+                str_ct_struct = str_ct_struct[:-2]                
+            event_str_elems.append('%s@%s'%(str_ct_struct, 
+                                    str(representative_resolved_PDGs).replace(' ','')))
+        
+        if len(event_ct_structs)>1:
+            event_str = '( %s )'%(' + '.join(event_str_elems))
+        else:
+            event_str = event_str_elems[0]
+        input_mapping = self.get_information('input_mapping')
+        if input_mapping is not None:
+             event_str += 'x{%s}'%(','.join('%d:%d'%(k,input_mapping[k]) for k in
+                                                        sorted(input_mapping.keys()) ))                 
+        if self.is_a_mirrored_event:
+            event_str = '<->%s'%event_str
+        beam_convolution_masks = self.get_information('beam_convolution_masks') 
+        if beam_convolution_masks is not None:
+            for (beam_name, beam_short_name) in [('beam_one','F1'),('beam_two','F2')]:
+                if beam_convolution_masks[beam_name] != 'ALL':
+                    event_str += '@%s<-(%s)'%(beam_short_name,'|'.join('%d'%pdg for pdg in
+                                                    beam_convolution_masks[beam_name]))
+        return event_str
+
+class ME7EventList(list):
+    """ A class handling a collection of ME7Events, with helper function to 
+    collectively act on them. This customized list is mainly intended to
+    store a list of ME7Events that all originate from the same integrand 
+    evaluation call."""
+    
+    def __init__(self, *args, **opts):
+        super(ME7EventList,self).__init__(*args, **opts)
+
+    def apply_PDF_convolution(self, *args,**opts):
+        """ Convolute the weights_per_flavor_configurations of each event 
+        with the PDF luminosity."""
+        for event in self:
+            event.apply_PDF_convolution(*args, **opts)
+
+    def filter_flavor_configurations(self, flavor_cut, **opts):
+        """ Apply the flavor cut on each event in this list, removing those with
+        no valid flavor configuration left. This returns false if there is no events left."""
+        
+        new_event_list = []
+        for event in self:
+            if event.filter_flavor_configurations(flavor_cut, **opts):
+                new_event_list.append(event)
+        self[:] = new_event_list
+
+        return len(self) != 0
+
+    def filter_with_flavor_blind_cuts(self, cut_function, *args, **opts):
+        """ Apply the kinematic flavor-blind cut on each event in this list, removing those with
+        failing the cut. This returns false if there is no events left."""
+        
+        new_event_list = []
+        for event in self:
+            if event.apply_flavor_blind_cuts(cut_function, *args, **opts):
+                new_event_list.append(event)
+        self[:] = new_event_list
+
+        return len(self) != 0
+    
+    def get_total_weight(self):
+        """ Return the total weight over all events in this list."""
+        if len(self) > 0:
+            return sum(event.get_total_weight() for event in self)
+        else:
+            return 0.
+
+    def select_a_flavor_configuration(self):
+        """ Select a particular flavor configuration to be used in the flavor-sensitive
+        observables for each event. """
+        
+        for event in self:
+            event.select_a_flavor_configuration()
+    
+    def select_epsilon_expansion_term(self, term_name):
+        """ Select a particular EpsilonExpansion term in the weights of the flavor 
+        matrix of these events. If the weights are float already, then they will be
+        set to zero unless the term_name is 'finite'."""
+        
+        for event in self:
+            event.select_epsilon_expansion_term(term_name)
+    
+    def apply_observables(self, observable_list, integrator_weight=1.):
+        """ Applies the observable list to the entire set of events of this list."""
+
+        for event in self:
+            observable_list.apply_observables(event, integrator_weight)
+
+    def generate_mirrored_events(self):
+        """ Creates the mirror configurations if necessary."""
+        for event in list(self):
+            mirrored_event = event.generate_mirrored_event()
+            if mirrored_event:
+                self.append(mirrored_event)
+
+    def __str__(self):
+        return self.nice_string()
+
+    def nice_string(self):
+        """ Returns a nice string representation of this list of event."""
+        
+        BLUE = misc.bcolors.BLUE
+        GREEN = misc.bcolors.GREEN
+        ENDC = misc.bcolors.ENDC
+        
+        n_events = len(self)
+        res=['%s>> Start of a list of %d event%s:%s'%(
+                                        GREEN, n_events, 's' if n_events > 1 else '',ENDC)]
+        res.append('')
+        for i_event, event in enumerate(self):
+             res.append('#%d '%(i_event+1) + event.nice_string())
+             res.append('')
+        res.append('%s>> End of a list of %d event%s.%s'%(
+                                        GREEN, n_events, 's' if n_events > 1 else '',ENDC))
+        return '\n'.join(res)
+
+#===============================================================================
 # ME7Integrand
 #===============================================================================
 class ME7Integrand(integrands.VirtualIntegrand):
@@ -98,6 +642,8 @@
     # Maximum size of the cache for PDF calls
     PDF_cache_max_size = 1000
     
+    FROZEN_DIMENSIONS = {}#{'x_1': 0.13, 'x_2': 0.26, 'x_3': 0.35}
+
     def __new__(cls, model, 
                      run_card,
                      contribution_definition,
@@ -108,40 +654,45 @@
                      ME7_configuration, **opt):
         all_args = [model, run_card, contribution_definition, processes_map,
                     all_MEAccessors, ME7_configuration ]
+
         if cls is ME7Integrand:
             target_type = 'Unknown'
-            if contribution_definition.correction_order == 'LO':
-                if contribution_definition.n_loops == 0 and \
-                  contribution_definition.n_unresolved_particles == 0:
-                   target_type = 'Born'
-                elif contribution_definition.n_loops == 1 and \
-                  contribution_definition.n_unresolved_particles == 0:
-                   target_type = 'LoopInduced_Born'
-            elif contribution_definition.correction_order == 'NLO':
-                if contribution_definition.n_loops == 1 and \
-                   contribution_definition.n_unresolved_particles == 0:
+            correction_order = contribution_definition.correction_order.count('N')
+            beam_factorization_count = 0
+            if contribution_definition.is_beam_active('beam_one'):
+                beam_factorization_count += 1
+            if contribution_definition.is_beam_active('beam_two'):
+                beam_factorization_count += 1
+            n_loops = contribution_definition.n_loops
+            n_unresolved_particles = contribution_definition.n_unresolved_particles
+            # Beam factorization contributions are automatically of type RV because
+            # they must both generate local counterterms (for the form factors) and
+            # accept integrated ISR ones.
+            if beam_factorization_count > 0:
+                if contribution_definition.correlated_beam_convolution:
+                    target_type = 'BeamSoft'
+                else:
+                    target_type = 'RealVirtual'                          
+            elif n_loops == 0 and n_unresolved_particles == 0:
+                target_type = 'Born'
+            elif n_loops == 1 and n_unresolved_particles == 0:
+                if correction_order < 1:
+                    target_type = 'LoopInduced_Born'
+                else:
                     target_type = 'Virtual'
-                elif contribution_definition.n_loops == 0 and \
-                     contribution_definition.n_unresolved_particles == 1:
-                    target_type = 'SingleReals'
-            elif contribution_definition.correction_order == 'NNLO':
-                if contribution_definition.n_loops == 0 and \
-                   contribution_definition.n_unresolved_particles == 2:
-                    target_type = 'DoubleReals'            
-                else:
-                    raise MadGraph5Error("Some NNLO type of integrands are not implemented yet.")   
-            elif contribution_definition.correction_order == 'NNNLO':
-                if contribution_definition.n_loops == 0 and \
-                   contribution_definition.n_unresolved_particles == 3:
-                    target_type = 'TripleReals'            
-                else:
-                    raise MadGraph5Error("Some NNNLO type of integrands are not implemented yet.") 
+            elif n_loops == 0 and n_unresolved_particles == 1:
+                target_type = 'SingleReals'
+            elif n_loops == 0 and n_unresolved_particles == 2:
+                target_type = 'DoubleReals'
+            elif n_loops == 0 and n_unresolved_particles == 3:
+                target_type = 'TripleReals'                
             else:
-                target_type = 'Unknown'
+                raise MadGraph5Error("Some %s type of ME7Integrands are not implemented yet."%
+                                                  contribution_definition.correction_order)
             target_class = ME7Integrand_classes_map[target_type]
 
             if not target_class:
-                raise MadGraph5Error("Could not determine the class of integrand of type '%s' to be added for"%target_type+
+                raise MadGraph5Error("Could not determine the class of ME7Integrand of type '%s' to be added for"%target_type+
                                      " the contribution definiton:\n%s"%str(contribution_definition.nice_string()))
 
             return super(ME7Integrand, cls).__new__(target_class, *all_args, **opt)
@@ -176,6 +727,7 @@
         
         # The original ContributionDefinition instance at the origin this integrand 
         self.contribution_definition    = contribution_definition
+
         # The process map of the Contribution instance at the origin of this integrand.
         # The format is identical to the one generated from the function 'get_process_map' of a contribution.
         self.processes_map              = processes_map
@@ -197,13 +749,17 @@
         # This counter is incremented for each time self.__call__ is called and reinitialized in self.synchronize
         self.n_observable_calls = 0
 
-    def nice_string(self):
-        """ For now simply use the contribution_definition and class name for a nice readable representation."""
-        
-        res = []
-        res.append("Instance of class '%s', with the following contribution definition:"%(self.__class__.__name__))
-        res.append('\n'.join(' > %s'%line for line in self.contribution_definition.nice_string().split('\n')))
-        return '\n'.join(res)
+        # Store a flavor cut definition, which the ME7 interface can overwrite. This will be used
+        # by the function pass_flavor_sensitive cut and it is always None by default.
+        self.flavor_cut_function = None 
+
+    def has_integrated_counterterms(self):
+        """ A property of the class, indicating whether it can contain integrated counterterms."""
+        return False
+
+    def has_local_counterterms(self):
+        """ A property of the class, indicating whether it can contain local counterterms."""
+        return False
 
     def get_additional_nice_string_printout_lines(self):
         """ Additional printout information for nice_string. 
@@ -218,12 +774,18 @@
         return GREEN+'  %s'%defining_process.nice_string(print_weighted=False).\
                                                                replace('Process: ','')+ENDC
 
+    def get_short_name(self):
+        """ Returns the short-name for this integrand, typically extracted from the one
+        of the contribution definition."""
+        return self.contribution_definition.short_name()
+
     def nice_string(self, format=0):
         """ Nice string representation of self."""
         BLUE = '\033[94m'
         GREEN = '\033[92m'
         ENDC = '\033[0m'
-        res = ['%-30s:   %s'%('ME7Integrand_type',type(self))]
+        res = ['< %s%s%s >'%(BLUE,self.get_short_name(),ENDC)]
+        res.append('%-30s:   %s'%('ME7Integrand_type',type(self)))
         res.extend([self.contribution_definition.nice_string()])
         if not self.topologies_to_processes is None:
             res.append('%-30s:   %d'%('Number of topologies', 
@@ -244,6 +806,155 @@
             
         return '\n'.join(res).encode('utf-8')
 
+    def __str__(self):
+        return self.nice_string()
+
+    @staticmethod
+    def build_flavor_cut_function(flavor_cut_string):
+        """ Given a string defining a flavor cut function specification (taken from the run
+        card), build here the flavor cut function to apply at run-time. Example of a complicated
+        cut function:
+          
+          "{
+            'R':'* 21 > (21,2)x3 * | * 21 > (21,3)x3 *',
+            'V':'21 (-2,2) > 21x2 (range(1,7)+range(-1,-7,-1))x5
+           }"
+
+          Means that the real must have:
+            a) anything as first incoming particle and a gluon the second incoming particle
+            b) Three 'gluons or up quarks' in the final state  + anything
+               *or* Three 'gluons or strange quarks' in the final state  + anything
+          
+          Means that the virtual must have:
+            a) a gluon and a 'up quark or anti-up quark' as first and second particle
+            b) 2 gluons and 5 quarks or antiquarks *exactly* (all particles must fit in one category at least)
+            
+          Notice that one can also speficy just a string, in which case the cut will apply
+          to all contributions.
+    
+        """
+
+        # First check if the format is a dictionary that specifies a different function
+        # for different integrand or if it applies to all integrands.
+        try:
+            flavor_func_dic = eval(flavor_cut_string)
+            is_a_dict = isinstance(flavor_func_dic, dict)
+        except:
+            is_a_dict = False
+        if is_a_dict:
+            if self.get_short_name() in flavor_func_dic:
+                flavor_cut_string = flavor_func_dic[self.get_short_name()]
+            else:
+                return None
+        
+        if flavor_cut_string.lower() == 'none':
+            return None
+        if flavor_cut_string.lower() == 'hardcoded':
+            return 'hardcoded'
+        
+        patterns = [p.strip() for p in flavor_cut_string.split('|')]
+        processed_patterns = []
+        for pattern in patterns:
+            processed_pattern = {}
+            try:
+                initial_state_pattern, final_state_pattern = pattern.split('>')
+                initial_state_pattern = initial_state_pattern.strip()
+                final_state_pattern = final_state_pattern.strip()
+            except:
+                raise MadEvent7Error('Malformed flavor-cut pattern: %s.'%pattern)
+            initial_state_patterns = [isp.strip() for isp in re.split(r'\s+',initial_state_pattern)]
+            processed_pattern['initial_states'] = []
+            for isp in initial_state_patterns:
+                if isp=='*':
+                    processed_pattern['initial_states'].append(None)
+                else:
+                    try:
+                        specifier = eval(isp)
+                    except:
+                        raise MadEvent7Error('Malformed flavor-cut pattern: %s.'%isp)
+                    if isinstance(specifier,int):
+                        processed_pattern['initial_states'].append((specifier,))
+                    elif isinstance(specifier, (tuple, list)):
+                        processed_pattern['initial_states'].append(tuple(specifier))
+                    else:
+                        raise MadEvent7Error('Malformed flavor-cut pattern: %s.'%isp)
+            processed_pattern['final_states'] = []
+            processed_pattern['needs_exact_final_states'] = True
+            final_state_patterns = [fsp.strip() for fsp in re.split(r'\s+',final_state_pattern)]
+            for fsp in final_state_patterns:
+                try:
+                    matcher, multiplier = fsp.split('x')
+                except:
+                    multiplier = '1'
+                    matcher = fsp
+                try:
+                    multiplier = int(multiplier.strip())
+                except:
+                    raise MadEvent7Error('Malformed flavor-cut pattern: %s.'%multiplier)
+                matcher = matcher.strip()
+                if matcher=='*':
+                    processed_pattern['needs_exact_final_states'] = False
+                else:
+                    try:
+                        specifier = eval(matcher)
+                    except:
+                        raise MadEvent7Error('Malformed flavor-cut pattern: %s.'%matcher)
+                    if isinstance(specifier,int):
+                        processed_pattern['final_states'].append(((specifier,),multiplier))
+                    elif isinstance(specifier, (tuple, list)):
+                        processed_pattern['initial_states'].append((tuple(specifier),multiplier))
+                    else:
+                        raise MadEvent7Error('Malformed flavor-cut pattern: %s.'%matcher)
+            processed_patterns.append(processed_pattern)
+            
+        # Ok, now that we've processed the patterns to consider, we can create the function
+        # that checks it:
+        def pass_flavor_cut_func(flavor_config):
+            initial_pdgs, final_pdgs = flavor_config[0], flavor_config[1]
+            # Loop over all capturing patterns
+            for pattern in processed_patterns:
+                # First check the initial states
+                if len(initial_pdgs)!=len(pattern['initial_states']):
+                    continue
+                must_continue = False
+                for i, pdg in enumerate(initial_pdgs):
+                    if (pattern['initial_states'][i] is not None) and \
+                                    (pdg not in pattern['initial_states'][i]):
+                        must_continue = True
+                        break
+                if must_continue:
+                    continue
+                
+                # Then check the final states.
+                # First make sure all constraints specified by the user apply
+                must_continue = False
+                for (matcher, multiplier) in pattern['final_states']:
+                    if len([1 for pdg in final_pdgs if pdg in matcher])!=multiplier:
+                        must_continue = True
+                        break
+                if must_continue:
+                    continue
+                # And if 'needs_exact_final_states' then also make sure that all final
+                # states belong in at least one category
+                if pattern['needs_exact_final_states']:
+                    must_continue = False
+                    for pdg in final_pdgs:
+                        if not any(pdg in matcher for (matcher, multiplier) in pattern['final_states']):
+                            must_continue = True
+                            break
+                    if must_continue:
+                        continue
+                
+                # If we haven't "continued" by now, then it means that the pattern was
+                # capturing and we can return True
+                return True
+                
+            # No pattern matched upt to this point so we must return False
+            return False
+            
+        # We can now return the flavor cut function created
+        return pass_flavor_cut_func
+                
     def synchronize(self, model, run_card, ME7_configuration):
         """ Synchronize this integrand with the most recent run_card and model."""
 
@@ -258,6 +969,8 @@
         # A RunCardME7 instance, properly initialized with the values of the run_card.dat of this run
         self.run_card                   = run_card
         
+        self.flavor_cut_function = ME7Integrand.build_flavor_cut_function(self.run_card['flavor_cuts'])
+        
         # Set external masses
         all_processes = [p[0] for p in self.processes_map.values()]
         self.masses = all_processes[0].get_external_masses(self.model)
@@ -279,22 +992,34 @@
                                                             (self.run_card['lpp1'], self.run_card['lpp2']))
         
         # Always initialize the basic flat PS generator. It can be overwritten later if necessary.
+        simplified_beam_types = (
+            0 if self.contribution_definition.beam_factorization['beam_one'] is None else 1,
+            0 if self.contribution_definition.beam_factorization['beam_two'] is None else 1,
+        )
         self.phase_space_generator = phase_space_generators.FlatInvertiblePhasespace(
             self.masses[0], self.masses[1],
-            beam_Es    = (self.run_card['ebeam1'], self.run_card['ebeam2']),
-            beam_types = (self.run_card['lpp1'], self.run_card['lpp2']),
+            beam_Es             = (self.run_card['ebeam1'], self.run_card['ebeam2']),
+            beam_types          = simplified_beam_types,
+            is_beam_factorization_active = 
+                        ( self.contribution_definition.is_beam_active('beam_one'),
+                          self.contribution_definition.is_beam_active('beam_two') ),
+            correlated_beam_convolution = self.contribution_definition.correlated_beam_convolution
         )
 
         # Add a copy of the PS generator dimensions here.
         # Notice however that we could add more dimensions pertaining to this integrand only, and PS generation.
-        # This is in particular true for discrete integration dimension like sectors, helicities, etc... 
-        self.set_dimensions(integrands.DimensionList(self.phase_space_generator.dimensions))
+        # This is in particular true for discrete integration dimension like sectors, helicities, etc...
+        integrand_dimensions = integrands.DimensionList(self.phase_space_generator.dimensions)
+
+        # Now remove the frozen dimensions from the list of continuous ones.
+        integrand_dimensions = integrands.DimensionList(
+                            d for d in integrand_dimensions if d.name not in self.FROZEN_DIMENSIONS)
+        self.set_dimensions(integrand_dimensions)
         self.dim_ordered_names = [d.name for d in self.get_dimensions()]
         self.dim_name_to_position = dict((name,i) for i, name in enumerate(self.dim_ordered_names))
         self.position_to_dim_name = dict((v,k) for (k,v) in self.dim_name_to_position.items())
 
         self.collider_energy = self.run_card['ebeam1'] + self.run_card['ebeam2']
-        
         # Set the seed
         if self.run_card['iseed'] > 0:
             random.seed(self.run_card['iseed'])
@@ -462,8 +1187,8 @@
 
         return False
 
-    def find_counterterms_matching_regexp(
-        self, counterterms, limit_pattern=None):
+    def find_counterterms_matching_regexp(self,
+                                    counterterms, limit_pattern=None, integrated_CT=False):
         """ Find all mappings that match a particular limits given in argument
         (takes a random one if left to None). This function is placed here given that
         it can be useful for both the ME7Integrnd_V and ME7_integrand_R."""
@@ -471,7 +1196,11 @@
         # First select only the counterterms which are not pure matrix elements 
         # (i.e. they have singular structures) and also exclude here soft-collinear 
         # counterterms since they do not have an approach limit function.
-        selected_counterterms = [ ct for ct in counterterms if ct.is_singular() ]
+        if not integrated_CT:
+            selected_counterterms = [ ct for ct in counterterms if ct.is_singular() ]
+        else:
+            selected_counterterms = [ ct for ct in counterterms if ct['integrated_counterterm'].is_singular() ]
+
         if len(selected_counterterms)==0:
             return []
         
@@ -499,13 +1228,18 @@
                     list_limit_pattern = [limit_pattern]
                 new_list_limit_pattern = []
                 for limit_pattern in list_limit_pattern:
-                    if not limit_pattern.startswith('('):
+                    if not integrated_CT and not limit_pattern.startswith('('):
                         # If not specified as a raw string, we take the liberty of adding 
                         # the enclosing parenthesis.
                         limit_pattern = '(%s,)' % limit_pattern
+                    elif integrated_CT and not limit_pattern.startswith('['):
+                        # If not specified as a raw string, we take the liberty of adding 
+                        # the enclosing parenthesis.
+                        limit_pattern = '[%s,]' % limit_pattern                    
                     # We also take the liberty of escaping the parenthesis
                     # since this is presumably what the user expects.
                     limit_pattern = limit_pattern.replace('(', '\(').replace(')', '\)')
+                    limit_pattern = limit_pattern.replace('[', '\[').replace(']', '\]')
                     new_list_limit_pattern.append(limit_pattern)
                 limit_pattern_re = re.compile(r'^(%s)$'%(
                     '|'.join(limit_pattern for limit_pattern in new_list_limit_pattern) ))    
@@ -515,19 +1249,34 @@
             returned_counterterms.append(random.choice(selected_counterterms))
         else:
             for counterterm in selected_counterterms:
-                if re.match(limit_pattern_re, str(counterterm)):
+                ct = counterterm if not integrated_CT else counterterm['integrated_counterterm']
+                if re.match(limit_pattern_re, str(ct)):
                     returned_counterterms.append(counterterm)
 
         return returned_counterterms
 
-    def pass_flavor_blind_cuts(self, PS_point, process_pdgs, n_jets_allowed_to_be_clustered = None):
+    def pass_all_cuts(self, PS_point, flavors, 
+                          n_jets_allowed_to_be_clustered = None, xb_1 = None, xb_2 = None):
+        """ Calls both the flavor-blind and flavor-sensitive cuts."""
+        return self.pass_flavor_blind_cuts(PS_point, flavors, 
+                    n_jets_allowed_to_be_clustered = n_jets_allowed_to_be_clustered, 
+                    xb_1 = xb_1, xb_2 = xb_2 ) and \
+                self.pass_flavor_sensitive_cuts(PS_point, flavors, 
+                    xb_1 = xb_1, xb_2 = xb_2 )
+
+    def pass_flavor_blind_cuts(self, PS_point, process_pdgs, 
+                          n_jets_allowed_to_be_clustered = None, xb_1 = None, xb_2 = None):
         """ Implementation of a minimal set of isolation cuts. This can be made much nicer in the future and 
         will probably be taken outside of this class so as to ease user-specific cuts, fastjet, etc...
         We consider here a two-level cuts system, this first one of which is flavour blind.
         The 'n_jets_allowed_to_be_clustered' is an option that allows to overwrite the
         maximum number of jets that can be clustered and which is by default taken to be:
             self.contribution_definition.n_unresolved_particles 
-        This is useful when using this function to apply cuts to the reduced PS of the CTs."""
+        This is useful when using this function to apply cuts to the reduced PS of the CTs.
+        Finally, the options 'xb_1' and 'xb_2' allow to specify the boost bringing
+        the PS point back from the c.o.m. to the lab frame, necessary if some cuts are not
+        invariant under a boost along the beam axis."""
+
         debug_cuts = False
         # This is a temporary function anyway which should eventually be replaced by a full
         # fledged module for handling generation level cuts, which would also make use of fjcore.
@@ -535,6 +1284,12 @@
 
         from madgraph.integrator.vectors import LorentzVectorDict, LorentzVectorList, LorentzVector
         
+        # If the cuts depend on the boost to the lab frame in case of hadronic collision
+        # then the quantity below can be used:
+#        boost_vector_to_lab_frame = None
+#        if (xb_1 is not None) and (xb_2 is not None) and (xb_1, xb_2)!=(1.,1.):
+#            boost_vector_to_lab_frame = PS_point.get_boost_vector_to_lab_frame(xb_1, xb_2)
+        
         # These cuts are not allowed to resolve flavour, but only whether a particle is a jet or not
         def is_a_jet(pdg):
             return abs(pdg) in range(1,self.run_card['maxjetflavor']+1)+[21]
@@ -554,6 +1309,19 @@
         if n_jets_allowed_to_be_clustered is None: