def run_L_D_simulation(self, L, D):
    # L = duplication length
    # D = number of DCJs in each branch.
    #
    param = self.sim_parameters
    # pre_dups (at root) and post_dups (at branches) to achieve 1.5 genes/family in average.
    pre_duplications = int(0.43 * param.num_genes / L)
    post_duplications = int(0.07 * param.num_genes / L)
    post_duplications = [int(0.5 * post_duplications), int(1.5 * post_duplications)]
    # post_duplications = [int(1 * post_duplications), int(1 * post_duplications)]
    param.pre_duplications = pre_duplications
    current_copy_number = None  # will init at root
    deletion_length_range = xrange(1, param.indel_length + 1)
    duplication_length_range = xrange(1, L + 1)
    idx = 1
    ev_tree = self.sim_tree
    for ev_node in ev_tree.preorder_node_iter():
        if ev_node.parent_node is None:
            # identity genome:
            ev_node.value = current_genome = model.Genome.identity(param.num_genes, param.num_chr)
            ev_node.events = {ev: 0 for ev in EventType.all}

            # add copy number information to track orthologous/paralogous, when duplications are present:
            for chromosome in current_genome.chromosomes:
                chromosome.copy_number = [1] * len(chromosome.gene_order)
            current_copy_number = current_genome.gene_count()
            # pre-duplications:
            for i in range(pre_duplications):
                Simulation.apply_random_segmental_duplication(current_genome,
                                                              range(1, param.duplication_length + 1),
                                                              current_copy_number)

            ev_node.events[EventType.DUPLICATION] = pre_duplications
            # ev_node.edge.length = pre_duplications

            if ev_node.label is None:
                ev_node.label = "Root"
        else:
            # evolve genome:
            if ev_node.is_internal():
                if ev_node.label is None:
                    ev_node.label = "M%02d" % idx
                    idx += 1
            else:  # complete labelling for leaves
                ev_node.label = ev_node.taxon.label

            current_genome = ev_node.parent_node.value.clone(ev_node.label)
            ev_node.value = current_genome
            pd = post_duplications.pop()
            ev_node.edge.length = D + pd

            # events
            events = [EventType.DUPLICATION] * pd + [EventType.REARRANGEMENT] * D

            ev_node.edge.events = {ev: 0 for ev in EventType.all}
            random.shuffle(events)
            for event in events:
                if event == EventType.DUPLICATION:
                    Simulation.apply_random_segmental_duplication(current_genome, duplication_length_range, current_copy_number)
                    ev_node.edge.events[event] += 1
                elif event == EventType.REARRANGEMENT:
                    # here, I can also have deletions:
                    ev = np.random.choice([RearrangementType.REVERSAL, EventType.DELETION], 1,
                                          p=[param.rearrangement_p, param.deletion_p])[0]
                    if ev == RearrangementType.REVERSAL:
                        Simulation.apply_random_reversal(current_genome)
                        ev_node.edge.events[event] += 1
                    else:
                        Simulation.apply_random_deletion(current_genome, deletion_length_range)
                        ev_node.edge.events[EventType.DELETION] += 1

            ev_node.events = {ev: ev_node.parent_node.events[ev] + count for ev, count in
                              ev_node.edge.events.iteritems()}