EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pepsiMaineRHIC_radcorr.v1.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pepsiMaineRHIC_radcorr.v1.F
1  program pepsimainerhic
2 
3  include "leptou.inc"
4  include "lujets.inc"
5  include "ludat1.inc"
6  include "mcRadCor.inc"
7  include "phiout.inc"
8  include "mc_set.inc"
9 
10  integer nevgen, i, j, lastgenevent, genevent
11  integer iev, nev, nprt
12  integer idum1, idum2, initseed
13  integer ltype, tracknr, nrtrack, t1, t2
14  integer*4 today(3), now(3)
15  real rval, alpha, alpha2
16  real mcfixedweight, mcextraweight, dilut, gendilut
17  real truey, trueq2, truex, truenu, truew2
18  real epx, epy, epz, eprim, ept, eptheta, epphi
19  double precision fdilut, f1, f2, r, da1, da2
20  DOUBLE PRECISION radgame, radgamp, radgamenucl, etheta, temp
21  DOUBLE PRECISION pbeame, pbeta, pgamma, ebeame, epznucl
22  DOUBLE PRECISION depol, gamma2, d, eta, eps, chi
23  DOUBLE PRECISION gendepol, gengamma2, gend, geneta, geneps
24  double precision genf1, genf2, genr, gena1, gena2, genchi
25  DOUBLE PRECISION dbrack, dxsec
26 
27  logical uselut, genlut
28  CHARACTER param*20
29 
30 c---------------------------------------------------------------------
31 c ascii output file
32 c ---------------------------------------------------------------------
33  integer asciilun
34  parameter(asciilun=29)
35  CHARACTER*256 outputfilename
36  CHARACTER*256 outname
37 
38  alpha = 1./137.
39  alpha2 = (1./137.)**2
40 c initialize event generation counter
41  nevgen = 0
42  genevent = 0
43  lastgenevent = 0.
44  dilut= 1.
45  gendilut= 1.
46  r = 0.
47  depol =0.
48  da1 = 0.
49  da2 = 0.
50  f1 = 0.
51  f2 = 0.
52  mcfixedweight=1.
53  mcextraweight=1.
54  weight=1.
55 c init all the parameters and read them in
56  ltype = 11
57  pbeam = 100.
58  ebeam = 4.
59 
60 c...Read parameters for linit call
61 c...Read output file name
62  READ(*,*) outname
63 c...Read beam target particle type
64  READ(*,*) ltype
65 c...Read beam and target particle energy
66  READ(*,*) pbeam, ebeam
67 c...Read number of events to generate, and to print.
68  READ(*,*) nev,nprt
69 c...Read min/max x of radgen lookup table
70  READ(*,*) mcset_xmin, mcset_xmax
71 c...Read min/max y of generation range
72  READ(*,*) mcset_ymin, mcset_ymax
73 c...Read min/max q2 of generation range
74  READ(*,*) mcset_q2min, mcset_q2max
75 c...Read information for cross section used in radgen
76  READ(*,*) genset_fstruct, genset_r
77 c...Read switch for rad-corrections
78  READ(*,*) qedrad
79 c...Read min/max w2 of generation range
80  READ(*,*) mcset_w2min,mcset_w2max
81 c...Read min/max nu of generation range
82  READ(*,*) mcset_numin,mcset_numax
83 c...Read beam and Target polarisation direction
84  READ(*,*) mcset_pbvalue, mcset_ptvalue
85 c...Read beam and Target polarisation
86  READ(*,*) mcset_pbeam, mcset_ptarget
87 c...Read all the pepsi/lepto parameters
88  100 READ(*,*, end=200) param, j, rval
89  if (param.eq.'LST') then
90  lst(j)=rval
91  write(*,*) 'LST(',j,')=',rval
92  endif
93  if (param.eq.'PARL') then
94  parl(j)=rval
95  write(*,*) 'PARL(',j,')=',rval
96  endif
97  if (param.eq.'PARJ') then
98  parj(j)=rval
99  write(*,*) 'PARJ(',j,')=',rval
100  endif
101  if (param.eq.'MSTJ') then
102  mstj(j)=rval
103  write(*,*) 'MSTJ(',j,')=',rval
104  endif
105  if (param.eq.'MSTU') then
106  mstu(j)=rval
107  write(*,*) 'MSTU(',j,')=',rval
108  endif
109  if (param.eq.'PARU') then
110  paru(j)=rval
111  write(*,*) 'PARU(',j,')=',rval
112  endif
113  goto 100
114  200 write(*,*) '*********************************************'
115  write(*,*) 'NOW all parameters are read by PEPSI'
116  write(*,*) '*********************************************'
117 
118 c ... pepsi polarization
119  if (mcset_pbvalue*mcset_ptvalue.lt.0.) then
120  lst(40) = 1.
121  elseif (mcset_pbvalue*mcset_ptvalue.gt.0.) then
122  lst(40) = -1.
123  else
124  lst(40) = 0.
125  endif
126 
127 c ... polarization for helium 3 !
128 c use lst(39)=1 for ideal helium 3 (no s prime and d state)
129 c use lst(39)=3 or (02 should work too ?) for a real helium 3.
130 c lst(39) = 0
131  if (lst(40).ne.0.and.parl(1).eq.3.and.parl(2).eq.2) then
132  lst(39) = 3
133  lst(36) = 1
134  endif
135 
136 c ... xmin,xmax
137  cut(1)=mcset_xmin
138  cut(2)=mcset_xmax
139 c ... ymin,ymax
140  cut(3)=mcset_ymin
141  cut(4)=mcset_ymax
142 c ... q2min, q2max
143  cut(5)=mcset_q2min
144  cut(6)=mcset_q2max
145 c ... w2min, w2max
146  cut(7)=mcset_w2min
147 c cut(8)=mcset_w2max
148 c ... numin, numax
149 c cut(9)=mcset_numin
150 c cut(10)=mcset_numax
151 
152  mcset_tara=parl(1)
153  mcset_tarz=parl(2)
154 
155 c getting the date and time of the event generation
156  call idate(today) ! today(1)=day, (2)=month, (3)=year
157  call itime(now) ! now(1)=hour, (2)=minute, (3)=second
158 c
159 c take date as the seed for the random number generation
160  initseed = today(1) + 10*today(2) + today(3) + now(1) + 5*now(3)
161  write(6,*) 'SEED = ', initseed
162  call rndmq(idum1,idum2,initseed,' ')
163 
164  if (abs(ltype).eq.11) then
165  masse=0.000510998910
166  else
167  write(*,*) "Only positron and electron beam supported"
168  endif
169 
170  if ((parl(1).eq.1).and.(parl(2).eq.1)) then
171  massp=0.938272013
172  elseif ((parl(1).eq.1).and.(parl(2).eq.0)) then
173  massp=0.93956536
174  elseif ((parl(2).eq.2).and.(parl(2).eq.1)) then
175  massp=1.875612793
176  else
177  write(*,*) "Currently only p,n,d as target supported"
178  endif
179 
180  pbeame=sqrt(pbeam**2+massp**2)
181  pbeta=pbeam/pbeame
182  pgamma=pbeame/massp
183  ebeame=sqrt(ebeam**2+masse**2)
184  ebeamenucl=pgamma*ebeame-pgamma*pbeta*(-ebeam)
185  epznucl=-pgamma*pbeta*ebeame+pgamma*(-ebeam)
186 c write(*,*) ebeamenucl, ebeame, epznucl, -ebeam
187  mcset_enebeam=sngl(ebeamenucl)
188 
189 c initialize pepsi
190  CALL timex(t1)
191  call linit(0,ltype,(-ebeam),pbeam,lst(23))
192  CALL timex(t2)
193 
194  if (qedrad.eq.1) then
195  uselut=.true.
196  genlut=.false.
197  call radgen_init(uselut,genlut)
198  write(*,*) 'I have initialized radgen'
199  elseif (qedrad.eq.2) then
200  write(*,*) 'radgen lookup table will be generated'
201  uselut=.true.
202  genlut=.true.
203  call radgen_init(uselut,genlut)
204  goto 500
205  endif
206 
207 c ---------------------------------------------------------------------
208 c Open ascii output file
209 c ---------------------------------------------------------------------
210  outputfilename=outname
211  open(asciilun, file=outputfilename)
212  write(*,*) 'the outputfile will be named: ', outname
213 
214 c the total integrated xsec, in pbarns
215  if (qedrad.eq.1) then
216  mcfixedweight=(mcset_q2max-mcset_q2min)*(mcset_ymax-mcset_ymin)
217  else
218  mcfixedweight = 1.
219  endif
220 
221 ! ============== Generate until happy (past cuts) =================
222 
223 c this is what we write in the ascii-file
224 
225  write(29,*)' PEPSI EVENT FILE '
226  write(29,*)'============================================'
227  write(29,30)
228 30 format('I, ievent, genevent, process, subprocess, nucleon,
229  & struckparton, partontrck, trueY, trueQ2, trueX, trueW2,
230  & trueNu, FixedWeight, weight, Extraweight, Sigcor, dilut, F1,
231  & F2, A1, A2, R, Depol, d, eta, eps, chi, radgamEnucl, nrTracks')
232 
233  write(29,*)'============================================'
234  write(29,*)' I K(I,1) K(I,2) K(I,3) K(I,4) K(I,5)
235  & P(I,1) P(I,2) P(I,3) P(I,4) P(I,5) V(I,1) V(I,2) V(I,3)'
236  write(29,*)'============================================'
237 
238  DO 300 iev=1,nev
239 
240 c count generated events
241 10 nevgen = nevgen+1
242 
243  weight=1.
244  if (qedrad.eq.0) then
245  If (iev.eq.1) then
246  write(*,*) 'Generating full event with LEPTO/PEPSI.'
247  endif
248  DO i=1,mstu(4)
249  DO j=1,5
250  k(i,j)=0.
251  p(i,j)=0.
252  v(i,j)=0.
253  enddo
254  enddo
255  mcradcor_ebrems=0.
256  uselut=.false.
257  genlut=.false.
258  lst(7) = 1
259  call lepto
260  if (lst(21).ne.0) then
261  write(*,*) 'PEPSI failed, trying again,
262  & no rad-corrections'
263  goto 10
264  endif
265  da2=0.
266  da1=0.
267  f1=0.
268  f2=0.
269  dilut=1.
270  r=0.
271  depol=0.
272  gamma2=0.
273  d=0.
274  eta=0.
275  eps=0.
276  chi=0.
277  trueq2=q2
278  truenu=u
279  truex=x
280  truey=y
281  truew2=w2
282  genevent=nevgen-lastgenevent
283  endif
284 
285  if (qedrad.eq.1) then
286  If (iev.eq.1) then
287  write(*,*) 'Generating only hadrons with LEPTO/PEPSI.'
288  endif
289  mcradcor_ebrems = 0.
290  geny=mcset_ymin*(mcset_ymax/mcset_ymin)**rlu(1)
291  genq2=mcset_q2min*(mcset_q2max/mcset_q2min)**rlu(2)
292 c geny = (mcset_ymax-mcset_ymin)*rlu(1) + mcset_ymin
293 c genq2 = (mcset_q2max-mcset_q2min)*rlu(2) + mcset_q2min
294  gennu=geny*ebeamenucl
295  genx = genq2/(2.*massp*gennu)
296  genw2 = massp**2.+(genq2*(1./genx-1.))
297  geneprim = ebeamenucl - gennu
298  if ((genq2.lt.mcset_q2min).or.(genq2.gt.(2.*gennu*massp))
299  & .or.(genw2.lt.mcset_w2min)) then
300  goto 10
301  endif
302  temp=(genq2/(2.*ebeamenucl*geneprim))-1.
303  if ((temp.le.1.).and.(temp.ge.-1.)) then
304  etheta = acos(temp)
305  genthe = sngl(etheta)
306 
307  x = genx
308  y = geny
309  q2 = genq2
310  w2 = genw2
311  u = gennu
312  lst(2) = 2
313  lst(7) = 0
314  DO i=1,mstu(4)
315  DO j=1,5
316  k(i,j)=0.
317  p(i,j)=0.
318  v(i,j)=0.
319  enddo
320  enddo
321 
322  call lepto
323  epx=p(4,1)
324  epy=p(4,2)
325  epz=p(4,3)
326  eprim=p(4,4)
327  ept=sqrt(epx**2+epy**2+epz**2)
328  eptheta=acos(epz/ept)
329  epphi=atan(p(4,2)/p(4,1))
330  genphi=epphi
331  else
332  goto 10
333  endif
334 
335 ! If doing radcor, do it, and THEN generate final state
336  call radgen_event
337 
338 c write(*,*)'0:', geny, genq2, genx, genw2, gennu
339 c write(*,*) 'after radgen', mcradcor_ytrue,
340 c & mcradcor_q2true, mcradcor_xtrue, mcradcor_w2true,
341 c & mcradcor_nutrue, mcradcor_ebrems, mcradcor_sigcor
342 
343  if (mcradcor_ebrems.gt.0.) then
344  x = mcradcor_xtrue
345  y = mcradcor_ytrue
346  q2 = mcradcor_q2true
347  w2 = mcradcor_w2true
348  u = mcradcor_nutrue
349  weight = weight * mcradcor_sigcor
350  if ((q2 .lt. cut(5)) .or. (w2 .lt. cut(7))) then
351  write(*,*)'WARNING:Radiated event beyond PEPSI range'
352  write (*,*) ' Q2 ', q2, cut(5)
353  write (*,*) ' W2 ', w2, cut(7)
354  goto 10
355  endif
356  elseif (mcradcor_ebrems.eq.0.) then
357  x = genx
358  y = geny
359  q2 = genq2
360  w2 = genw2
361  u = gennu
362  weight = weight * mcradcor_sigcor
363  if ((q2 .lt. cut(5)) .or. (w2 .lt. cut(7))) then
364  write(*,*)'WARNING: Radiated event beyond PEPSI range'
365  write (*,*) ' Q2 ', q2, cut(5)
366  write (*,*) ' W2 ', w2, cut(7)
367  goto 10
368  endif
369  endif
370  lst(2) = 2
371  lst(7) = 1
372  DO i=4,mstu(4)
373  DO j=1,5
374  k(i,j)=0.
375  p(i,j)=0.
376  v(i,j)=0.
377  enddo
378  enddo
379 
380  call lepto
381 
382  p(4,1)=epx
383  p(4,2)=epy
384  p(4,3)=epz
385  p(4,4)=eprim
386 
387  if (lst(21).ne.0) then
388  write(*,*) 'Gen Hadrons: PEPSI failed, trying again'
389  goto 10
390  endif
391  genevent=nevgen-lastgenevent
392  endif
393 
394  DO i=1,mstu(4)
395  if ((k(i,1).eq.0).and.(k(i,2).eq.0).and.(k(i,3).eq.0).and.
396  & (k(i,4).eq.0).and.(k(i,5).eq.0)) then
397  tracknr=i-1
398  if (mcradcor_ebrems.gt.0.) then
399  nrtrack=tracknr+1
400  else
401  nrtrack=tracknr
402  endif
403  goto 400
404  endif
405  ENDDO
406 400 continue
407 
408 c get some info for the event
409  call mkasym(dble(q2), dble(x), mcset_tara, mcset_tarz, da1, da2)
410  call mkasym(dble(genq2), dble(genx), mcset_tara, mcset_tarz,
411  & gena1, gena2)
412  call mkf2(dble(q2), dble(x), mcset_tara, mcset_tarz, f2, f1)
413  call mkf2(dble(genq2), dble(genx), mcset_tara, mcset_tarz,
414  & genf2, genf1)
415  call mkr(dble(q2), dble(x), r)
416  call mkr(dble(genq2), dble(genx), genr)
417  dilut = real(fdilut(dble(Q2), dble(X), mcset_tara))
418  gendilut = real(fdilut(dble(genQ2), dble(genX), mcset_tara))
419 
420 c... calculate some kinematic variables
421  gengamma2=genq2/gennu**2
422  gamma2=q2/u**2
423 
424  gendepol=(geny*(1.+gengamma2*geny/2.)*(2.-geny)-2*geny**2*
425  & masse**2/genq2)/(geny**2*(1.-2.*masse**2/genq2)*
426  & (1.+gengamma2)+2.*(1.+genr)*
427  & (1.-geny-gengamma2*geny**2/4.))
428  depol=(y*(1.+gamma2*y/2.)*(2.-y)-2*y**2*masse**2/q2)/
429  & (y**2*(1.-2.*masse**2/q2)*(1.+gamma2)+2.*(1.+r)*
430  & (1.-y-gamma2*y**2/4.))
431 
432  gend=gendepol*((sqrt(1.-geny-gengamma2*y**2/4.)*
433  & (1+gengamma2*geny/2.))/(1.-geny/2.)*(1+gengamma2*geny/2.)-
434  & geny**2*masse**2/genq2)
435  d=depol*((sqrt(1.-y-gamma2*y**2/4.)*(1+gamma2*y/2.)) /
436  & (1.-y/2.)*(1+gamma2*y/2.)-y**2*masse**2/q2)
437 
438  geneta=(sqrt(gengamma2)*(1.-geny-gengamma2*geny**2/4.-
439  & geny**2*masse**2/genq2)) / ((1.+gengamma2**2*geny/2.)*
440  & (1-geny/2.)-geny**2*masse**2/genq2)
441  eta=(sqrt(gamma2)*(1.-y-gamma2*y**2/4.-y**2*masse**2/q2)) /
442  & ((1.+gamma2**2*y/2.)*(1-y/2.)-y**2*masse**2/q2)
443 
444  genchi=(sqrt(gengamma2)*(1.-geny/2.-geny**2.*masse**2/genq2))/
445  & (1.+gengamma2*geny/2.)
446  chi=(sqrt(gamma2)*(1.-y/2.-y**2.*masse**2/q2))/(1.+gamma2*y/2.)
447 
448  geneps=1./(1.+1./2.*(1.-2.*masse**2/genq2)*
449  & ((geny**2+gengamma2*geny**2)/
450  & (1.-geny-gengamma2*geny**2/4.)))
451  eps=1./(1.+1./2.*(1.-2.*masse**2/q2)*((y**2+gamma2*y**2)/
452  & (1.-y-gamma2*y**2/4.)))
453 
454  dbrack = (4.*paru(1)*alpha**2)/genq2**2
455  dxsec = dbrack*((1.-geny-(massp*genx*geny/2./ebeame))*
456  & genf2/genx+geny**2*genf1)*genx*genq2
457  dxsec = dxsec*(1.- gendilut*(mcset_pbvalue*mcset_ptvalue)*
458  & (gendepol*(gena1+geneta*gena2)))
459 
460  if (qedrad.eq.1) weight = weight*dxsec
461 c write(*,*)'rad:', genq2, parl(23), parl(24), weight,
462 c & mcradcor_sigcor, dxsec
463 
464  if (mcradcor_ebrems.gt.0.) then
465  radgamenucl=sqrt(dplabg(1)**2+dplabg(2)**2+dplabg(3)**2)
466  radgame=pgamma*radgamenucl+pgamma*pbeta*(-dplabg(3))
467  radgamp=+pgamma*pbeta*radgamenucl+pgamma*(-dplabg(3))
468 c write(*,*)'photon:', radgamenucl, radgame, dplabg(3), radgamp
469  else
470  radgamenucl=0d0
471  radgame=0d0
472  radgamp=0d0
473  endif
474 
475  if ((parl(23).gt.0.).and.(parl(24).gt.0.)) then
476  mcextraweight = parl(24)/parl(23)
477  else
478  mcextraweight = 1.
479  endif
480 
481  write(29,32) 0, iev, genevent, lst(23), lst(24), lst(22),
482  & lst(25), lst(26), y, q2, x, w2, u, mcfixedweight, weight,
483  & mcextraweight, mcradcor_sigcor, dilut, real(F1), real(F2),
484  & real(dA1), real(dA2), real(R), real(Depol), real(d),
485  & real(eta), real(eps), real(chi), radgamenucl, nrtrack
486  32 format((i4,1x,$),(i10,1x,$),6(i6,1x,$),5(f12.6,1x,$),
487  & 2(f18.6,1x,$), 14(f12.6,1x,$),i12,/)
488  write(29,*)'============================================'
489 
490  DO i=1,tracknr
491  write(29,34) i,k(i,1),k(i,2),k(i,3),k(i,4),k(i,5),
492  & p(i,1),p(i,2),p(i,3),p(i,4),p(i,5),
493  & v(i,1),v(i,2),v(i,3)
494  ENDDO
495  if (mcradcor_ebrems.gt.0.) then
496  write(29,34) nrtrack, 55, 22, 1, 0, 0,
497  & sngl(dplabg(1)),sngl(dplabg(2)),sngl(radgamp),
498  & sngl(radgame), 0., 0., 0., 0.
499  34 format(2(i6,1x,$),i10,1x,$,3(i6,1x,$),8(f15.6,1x,$),/)
500  endif
501  write(29,*)'=============== Event finished ==============='
502 
503  lastgenevent=nevgen
504 
505  300 CONTINUE
506 
507  500 if (qedrad.eq.2) then
508  write(*,*) 'lookup table is generated;'
509  write(*,*) 'to run now PEPSI change parameter qedrad to 1'
510  endif
511 
512 c...print the pythia cross section which is needed to get an absolut
513 c normalisation the number is in pbarns
514  write(*,*)'==================================================='
515  write(*,*)'Pepsi total cross section in pb from numerical',
516  & ' integration', parl(23)
517  write(*,*)'Pepsi total cross section in pb from MC simulation',
518  & parl(24)
519  write(*,*)'mcextraweight =', mcextraweight
520  write(*,*)'Total Number of generated events', nev
521  write(*,*)'Total Number of trials', nevgen
522  write(*,*)'==================================================='
523  close(29)
524  END