LOEventGenerationBias: PY8_CKKW.f

File PY8_CKKW.f, 6.0 KB (added by Valentin, 4 years ago)
Line 
1C ************************************************************
2C Source for the library implementing a bias function that
3C implements the Sudakov weight in CKKW directly from Pythia8
4C ************************************************************
5C
6C The following lines are read by MG5aMC to set what are the
7C relevant parameters for this bias module.
8C
9C  parameters = {'BeamA': 0.0,
10C                'BeamB': 0.0}
11C
12
13      subroutine bias_wgt(p, original_weight, bias_weight)
14          implicit none
15C
16C Parameters
17C
18          include '../../maxparticles.inc'         
19          include '../../nexternal.inc'
20C
21C Accessingt the details of the event
22C
23          include '../../run_config.inc'
24          include '../../lhe_event_infos.inc'
25C
26C Arguments
27C
28          double precision p(0:3,nexternal)
29          double precision original_weight, bias_weight
30C
31C local variables
32C
33          integer lpp_to_beam(-3:3)
34          data lpp_to_beam/-11,-2212,-2212,0,2212,2212,11/
35c
36c local variables defined in the run_card
37c
38c         Bias module arguments
39          double precision BeamA, BeamB
40
41c         truly local variables
42          integer i,j, n_initial
43          double precision OutputBiasWeight
44          double precision Pythia8eCM
45          integer Pythia8nParticles
46          double precision Pythia8p(5,npart)
47          integer Pythia8BeamA
48          double precision PythiaBeamEnergyA
49          integer Pythia8BeamB
50          double precision PythiaBeamEnergyB
51          character*(maxEventLength) Pythia8EvtRecord
52          integer Pythia8Helicities(npart)
53          integer Pythia8ColorOne(npart)
54          integer Pythia8ColorTwo(npart)         
55          integer Pythia8ID(npart)
56          integer Pythia8Status(npart)
57          integer Pythia8MotherOne(npart)
58          integer Pythia8MotherTwo(npart)         
59          integer Pythia8SubprocessGroup
60          integer Pythia8MurScale
61          integer Pythia8AlphaQCD
62          integer Pythia8AlphaQED
63
64C
65C Global variables
66C
67C
68C Mandatory common block to be defined in bias modules
69C
70          double precision stored_bias_weight
71          data stored_bias_weight/1.0d0/         
72          logical impact_xsec, requires_full_event_info
73C         We only want to bias distributions, but not impact the xsec.
74          data impact_xsec/.True./
75C         Pythia8 will need the full information for the event
76C          (color, resonances, helicities, etc..)
77          data requires_full_event_info/.True./
78          common/bias/stored_bias_weight,impact_xsec,
79     &                requires_full_event_info
80C
81C Access the value of the run parameters in run_card
82C
83          include '../../run.inc'
84          include '../../cuts.inc'
85C
86C Read the definition of the bias parameter from the run_card   
87C
88          include '../bias.inc'
89
90C --------------------
91C BEGIN IMPLEMENTATION
92C --------------------
93
94C        Let's initialize the PY8 variables describing the event
95         Pythia8eCM             = sqrt(4d0*ebeam(1)*ebeam(2))
96         Pythia8EvtRecord       = event_record
97         Pythia8SubprocessGroup = ngroup
98         Pythia8MurScale        = sscale
99         Pythia8AlphaQCD        = aaqcd
100         Pythia8AlphaQED        = aaqed
101         Pythia8nParticles      = npart
102         PythiaBeamEnergyA      = ebeam(1)
103         PythiaBeamEnergyB      = ebeam(2)
104         do i=1,npart
105           Pythia8ID(i)         = jpart(1,i)
106           Pythia8MotherOne(i)  = jpart(2,i)
107           Pythia8MotherTwo(i)  = jpart(3,i)
108           Pythia8ColorOne(i)   = jpart(4,i)
109           Pythia8ColorTwo(i)   = jpart(5,i)           
110           Pythia8Status(i)     = jpart(6,i)
111           Pythia8Helicities(i) = jpart(7,i)           
112           do j=1,4
113             Pythia8p(j,i)=pb(mod(j,4),i)
114           enddo
115           Pythia8p(5,npart)=pb(4,i)
116         enddo
117         Pythia8BeamA = lpp_to_beam(lpp(1))
118         Pythia8BeamB = lpp_to_beam(lpp(2))
119         n_initial = 0
120         do i=1,npart
121           if (Pythia8Status(i).eq.-1) then
122             n_initial = n_initial + 1
123             if (lpp(n_initial).eq.0) then
124               if (n_initial.eq.1) then
125                 Pythia8BeamA = Pythia8ID(i)
126               elseif (n_initial.eq.2) then
127                 Pythia8BeamB = Pythia8ID(i)
128                 exit
129               endif
130             endif
131             if (n_initial.eq.2) then
132               exit
133             endif
134           endif
135         enddo
136C        If this is a 1 > N decay event, then enforce beamIDs to match
137C        those specified in the event record.
138         if (n_initial.eq.1) then
139           Pythia8BeamB = 0
140           do i=1,npart         
141             if (Pythia8Status(i).eq.-1) then
142               Pythia8BeamA = Pythia8ID(i)
143               exit
144             endif
145           enddo
146         endif
147C        Make sure to enforce the user-choice of beam if specified.
148         if (idnint(BeamA).ne.0) then
149           Pythia8BeamA = idnint(BeamA)
150         endif
151         if (idnint(BeamB).ne.0) then
152           Pythia8BeamB = idnint(BeamB)
153         endif
154
155C        Call PY8 to derive the bias weight.
156         call py8_bias_weight( Pythia8eCM,
157     &                         Pythia8BeamA,
158     &                         PythiaBeamEnergyA,
159     &                         Pythia8BeamB,
160     &                         PythiaBeamEnergyB,
161     &                         Pythia8EvtRecord,
162     &                         Pythia8p,
163     &                         Pythia8nParticles,
164     &                         Pythia8MurScale,
165     &                         Pythia8AlphaQCD,
166     &                         Pythia8AlphaQED,
167     &                         Pythia8ID,
168     &                         Pythia8MotherOne,
169     &                         Pythia8MotherTwo,
170     &                         Pythia8ColorOne,
171     &                         Pythia8ColorTwo,
172     &                         Pythia8Status,
173     &                         Pythia8Helicities,
174     &                         OutputBiasWeight    )
175 
176          bias_weight = OutputBiasWeight
177           
178          return
179
180      end subroutine bias_wgt