LOEventGenerationBias: PY8_CKKW.f

File PY8_CKKW.f, 6.0 KB (added by Valentin Hirschi, 8 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