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