Warning, /pythia6/pythiaeRHIC/erhic/pythia6.f95 is written in an unsupported language. File is not indexed.
0001 ! ======================================================================
0002 ! eRHIC implementation of PYTHIA event generator.
0003 ! This file contains routines to read in a steer file, initialise
0004 ! PYTHIA and generate events.
0005 !
0006 ! The code is written in fixed format for compatibility with
0007 ! the included PYTHIA and Radgen files.
0008 !
0009 ! There is no program defined - the routines are called
0010 ! from C++ code, which provides the 'main' function.
0011 ! We therefore wrap routines and variables in a module to ease
0012 ! passing data back and forth between Fortran and C++.
0013 !
0014 ! A note on variable naming: any variables, routines or types that
0015 ! are visible to C++ code are written with underscores between words
0016 ! rather than using mixed case i.e.
0017 ! names_like_this rather than NamesLikeThis etc.
0018 ! This is because Fortran is case-insensitive so all variable
0019 ! and routine names appear as all lower case in C++. Using underscores
0020 ! to separate words gives greater readability on the C++ side.
0021 !
0022 ! \todo Port radiative photon output to C++ code.
0023 ! ======================================================================
0024 module pythia6
0025
0026 private ! Make things private by default
0027
0028 ! These variables are all public so they can be
0029 ! shared with C++ code.
0030
0031 integer, public:: printascii ! Flag to print ASCII output
0032
0033 integer, public:: genevent ! Number of generations to get the current event
0034 double precision, public:: trueX ! Bjorken-x
0035 double precision, public:: trueW2 ! Invariant mass of hadronic final state
0036 double precision, public:: trueNu ! Photon energy in proton rest frame
0037 double precision, public:: F2
0038 double precision, public:: F1
0039 double precision, public:: R
0040 double precision, public:: sigma_rad
0041 double precision, public:: SigRadCor
0042 double precision, public:: EBrems
0043 double precision, public:: radgamE ! Energy of radiated photon
0044 double precision, dimension(3), public:: radgamp ! Momentum of radiated photon
0045
0046 ! These need to be accessed by the module
0047 ! subroutines, but aren't needed externally.
0048
0049 ! Logical unit to which ASCII output is written,
0050 ! is ascii output is requested.
0051 integer, parameter:: asciiLun = 29
0052
0053 integer:: NEV ! Number of events to generate
0054 integer:: IEV ! Counts the number of the current event
0055 double precision:: pbeam ! Proton beam momentum (GeV/c)
0056 double precision:: ebeam ! Lepton beam momentum (GeV/c)
0057 double precision:: massp ! Proton mass
0058 double precision:: masse ! Lepton mass
0059 double precision:: pgamma ! Gamma factor (E/m) of beam proton
0060 double precision:: pbeta ! Beta (p/E) of beam proton
0061
0062 ! The following subroutines are defined:
0063 ! initialise()
0064 ! finish()
0065 ! reset()
0066 ! The following functions are defined:
0067 ! generate()
0068
0069 public:: initialise, finish, generate
0070 private:: reset
0071
0072 contains
0073
0074 ! ================================================================
0075 ! Initialisation routine.
0076 ! Sets module variables to default values.
0077 ! Reads input from steer file and directs PYTHIA to process it.
0078 ! Opens output ASCII file and prints header if ASCII output
0079 ! has been requested.
0080 ! Returns an integer with the following meaning:
0081 ! 0: Success
0082 ! <0: An error was encountered
0083 ! >0: Success, user requested only initialisation steps in steer
0084 ! If a value >0 is returned then event generation should not
0085 ! be performed - this indicates e.g. that radgen table
0086 ! initialisation was requested.
0087 ! ================================================================
0088 function initialise()
0089
0090 include 'pythia.inc' ! All PYTHIA commons blocks
0091 include "mc_set.inc"
0092 include "py6strf.inc"
0093 include "mcRadCor.inc"
0094 include "radgen.inc"
0095 include "phiout.inc"
0096
0097 ! Added by liang 1/6/12
0098 ! Switches for nuclear correction
0099 common/PYNUCL/INUMOD, CHANUM, ORDER
0100 save/PYNUCL/
0101
0102 ! Function variable
0103 integer:: initialise
0104
0105 double precision INUMOD,CHANUM
0106 ! -------------------------------------------------------------
0107 ! Arrays storing times used in setting seed.
0108 ! -------------------------------------------------------------
0109 integer, dimension(3):: today
0110 integer, dimension(3):: now
0111 integer:: ORDER
0112 integer:: NPRT, I, ltype
0113 integer:: idum1, idum2, initseed
0114 integer:: un,istat
0115 double precision:: sqrts
0116 double precision:: pbeamE, ebeamE, epznucl
0117 logical:: UseLut, GenLut
0118 ! -------------------------------------------------------------
0119 ! ASCII output file
0120 ! -------------------------------------------------------------
0121 character(len=256):: outname
0122 character(len=100) :: PARAM ! Parameters passed to pygive()
0123
0124 ! Set some default values before reading anything.
0125 initialise = -1
0126 genevent = 0
0127 pbeam = 100.
0128 ebeam = 4.
0129 masse = PYMASS(11)
0130 massp = PYMASS(2212)
0131 NEV = 0
0132 IEV = 0
0133 iModel=0
0134 ltype=11
0135
0136 ! Read the steer file passed on the command line.
0137
0138 ! Read output file name
0139 read(*, *) outname
0140 ! Read lepton beam type
0141 read(*, *) ltype
0142 ! Read parameters for PYINIT call (beam and target particle energy).
0143 read(*, *) pbeam, ebeam
0144 ! Read number of events to generate, and to print.
0145 read(*, *) NEV,NPRT
0146 ! Read min/max x of radgen lookup table
0147 read(*, *) mcSet_XMin, mcSet_XMax
0148 ! Read min/max y of generation range
0149 read(*, *) mcSet_YMin, mcSet_YMax
0150 ! Read min/max Q2 of generation range
0151 read(*, *) mcSet_Q2Min, mcSet_Q2Max
0152 ! Read information for cross section used in radgen
0153 read(*, *) genSet_FStruct, genSet_R
0154 ! Read parameters of radcorr:
0155 ! do radcorr (1), generate look-up table (2)
0156 read(*, *) qedrad
0157 ! Read parameters for PYTHIA-Model = which generation is done
0158 read(*, *) iModel
0159 ! Read target type mass and charge
0160 read(*, *) mcSet_TarA, mcSet_TarZ
0161 ! Read nuclear pdf parameter mass number A, charge number Z
0162 read(*, *) INUMOD, CHANUM
0163 ! Read nuclear pdf correction order
0164 read(*, *) ORDER
0165 ! Read information for cross section used in radgen
0166 100 read(*, '(A)', end=200) PARAM
0167 call PYGIVE(PARAM)
0168 goto 100
0169
0170 ! -------------------------------------------------------------
0171 ! Initialize PYTHIA.
0172 ! -------------------------------------------------------------
0173 200 write(*, *) '*********************************************'
0174 write(*, *) 'NOW all parameters are read by PYTHIA'
0175 write(*, *) '*********************************************'
0176
0177
0178 ! Getting the date and time of the event generation
0179 call idate(today) ! today(1)=day, (2)=month, (3)=year
0180 call itime(now) ! now(1)=hour, (2)=minute, (3)=second
0181 open(newunit=un, file="/dev/urandom", access="stream",
0182 + form="unformatted", action="read", status="old", iostat=istat)
0183 if (istat == 0) then
0184 read(un) initseed
0185 close(un)
0186 else
0187 ! Take date as the SEED for the random number generation
0188 ! Use constant seed for debug purposes
0189 ! \todo Need to change this back for release!!!
0190 initseed = today(1) + 10*today(2) + today(3) + now(1)
0191 + + 5*now(3)
0192 endif
0193 write(6,*) 'SEED = ', abs(initseed)
0194 call rndmq(idum1,idum2,abs(initseed),' ')
0195
0196 sqrts = sqrt(4. * pbeam * ebeam)
0197 write(*, *) '*********************************************'
0198 write(*, *) 'proton beam energy:', pbeam, 'GeV'
0199 write(*, *) 'lepton beam energy:', ebeam, 'GeV'
0200 write(*, *) 'resulting sqrt(s):', sqrts, 'GeV'
0201 write(*, *) '*********************************************'
0202
0203 ! proton is defined in positive z and as target
0204 P(2,1)=0.0
0205 P(2,2)=0.0
0206 P(2,3)=pbeam
0207 ! lepton is defined in negative z and as beam
0208 P(1,1)=0.0
0209 P(1,2)=0.0
0210 P(1,3)=-ebeam
0211 ! Find the masses of the beam particles
0212 if(mcSet_TarZ.eq.0) then
0213 massp=PYMASS(2112)
0214 else
0215 massp=PYMASS(2212)
0216 endif
0217 masse=PYMASS(ltype)
0218 ! Compute beam-dependent properties
0219 pbeamE=sqrt(pbeam**2+massp**2)
0220 pbeta=pbeam/pbeamE
0221 pgamma=pbeamE/massp
0222 ebeamE=sqrt(ebeam**2+masse**2)
0223 ebeamEnucl=pgamma*ebeamE-pgamma*pbeta*(-ebeam)
0224 epznucl=-pgamma*pbeta*(ebeamE)+pgamma*(-ebeam)
0225 write(*, *) ebeamEnucl, ebeamE, epznucl, -ebeam
0226 mcSet_EneBeam=sngl(ebeamEnucl)
0227
0228 ! Take actions dependent on radgen options
0229 if(iModel.eq.0) then
0230 UseLUT=.false.
0231 GenLUT=.false.
0232 qedrad=0
0233 MSTP(199)=0
0234 mcRadCor_EBrems=0.
0235 elseif(iModel.eq.1) then
0236 if(qedrad.eq.0) then
0237 mcRadCor_EBrems=0.
0238 UseLUT=.false.
0239 GenLUT=.false.
0240 MSTP(199)=1
0241 elseif(qedrad.eq.1) then
0242 mcRadCor_EBrems=0.
0243 UseLUT=.true.
0244 GenLUT=.false.
0245 MSTP(199)=1
0246 call radgen_init(UseLUT,GenLUT)
0247 write(*, *) 'I have initialized radgen'
0248 elseif(qedrad.eq.2) then
0249 write(*, *) 'radgen lookup table will be generated'
0250 mcRadCor_EBrems=0.
0251 UseLUT=.true.
0252 GenLUT=.true.
0253 MSTP(199)=1
0254 call radgen_init(UseLUT,GenLUT)
0255 write(*, *) 'lookup table is generated;'
0256 write(*, *) 'to run now pythia change parameter qedrad to 1'
0257 ! Set the return value to greater than zero to
0258 ! indiate that event generation should not proceed.
0259 initialise = 1
0260 return
0261 endif
0262 endif
0263
0264 ! Initialise PYTHIA beam species and momenta
0265 if((mcSet_TarZ.eq.1).and.(ltype.eq.11)) then
0266 call pyinit ('3MOM','gamma/e-','p+',WIN)
0267 elseif((mcSet_TarZ.eq.1).and.(ltype.eq.-11)) then
0268 call pyinit ('3MOM','gamma/e+','p+',WIN)
0269 elseif((mcSet_TarZ.eq.0).and.(ltype.eq.-11)) then
0270 call pyinit ('3MOM','gamma/e+','n0',WIN)
0271 elseif((mcSet_TarZ.eq.0).and.(ltype.eq.11)) then
0272 call pyinit ('3MOM','gamma/e-','n0',WIN)
0273 endif
0274
0275 ! If we ever want to simulate fixed target we need to change this
0276 ! win=ebeam
0277 ! call pyinit('fixt','gamma/e-','p+', WIN)
0278
0279 ! -------------------------------------------------------------
0280 ! Open ascii output file
0281 ! -------------------------------------------------------------
0282 if(printascii .ne. 0) then
0283 open(asciiLun, file=outname)
0284 write(*, *) 'the outputfile will be named:', outname
0285 ! This is what we write in the ascii-file
0286 write(asciilun,*)' PYTHIA EVENT FILE '
0287 write(asciilun,*)'============================================'
0288 write(asciilun,*) 'I, ievent, genevent, subprocess, nucleon,'//
0289 + ' targetparton, xtargparton, beamparton, xbeamparton,'//
0290 + ' thetabeamprtn, truey, trueQ2, truex, trueW2, trueNu, leptonphi,'//
0291 + ' s_hat, t_hat, u_hat, pt2_hat, Q2_hat, F2, F1, R, sigma_rad,'//
0292 + ' SigRadCor, EBrems, photonflux, t-diff, nrTracks'
0293 write(asciilun,*)'============================================'
0294
0295 write(asciilun,*)' I K(I,1) K(I,2) K(I,3) K(I,4) K(I,5)'//
0296 + ' P(I,1) P(I,2) P(I,3) P(I,4) P(I,5) V(I,1) V(I,2) V(I,3)'
0297 write(asciilun,*)'============================================'
0298 endif ! printascii
0299
0300 initialise = 0
0301
0302 end function initialise
0303 ! ================================================================
0304
0305
0306 ! ================================================================
0307 ! Generates a single event.
0308 ! Returns 0 upon success.
0309 ! Returns 1 upon an error, or if the total number of events
0310 ! (counted by NEV) requested has been reached.
0311 ! ================================================================
0312 function generate()
0313
0314 include 'pythia.inc' ! All PYTHIA commons blocks
0315 include "mc_set.inc"
0316 include "py6strf.inc"
0317 include "mcRadCor.inc"
0318 include "radgen.inc"
0319 include "phiout.inc"
0320
0321 integer:: generate ! Function return value
0322 integer:: nrtrack
0323 integer:: i ! loop variable
0324 ! Tracks the number of trials between successful event
0325 ! generations
0326 integer, save:: lastgenevent = 0
0327
0328 ! Return if we have reached the maximum event count
0329 if(.not. IEV .lt. NEV) then
0330 generate = 1
0331 return
0332 endif
0333
0334 ! Clear previous contents of publicly available variables
0335 call reset()
0336
0337 ! call PYTHIA event generation
0338 999 call PYEVNT
0339 if(MSTI(61).eq.1) then
0340 write(*, *) 'go back to PYEVNT call'
0341 goto 999
0342 endif
0343
0344 genevent = NGEN(0, 3) - lastgenevent
0345
0346 trueX = VINT(307)/VINT(309)/(4*pbeam*ebeam)
0347 trueW2 = massp**2 + VINT(307)*(1/trueX-1)
0348 trueNu = (trueW2 + VINT(307) - massp**2)/(2.*massp)
0349
0350 ! Calculate values for radiated photon
0351 if(mcRadCor_EBrems.gt.0.) then
0352 ! dplabg is from common/phiout/ in phiout.inc.
0353 EBrems=sqrt(dplabg(1)**2+dplabg(2)**2+dplabg(3)**2)
0354 ! Boost energy and momentum from lab frame.
0355 radgamE=pgamma*EBrems-pgamma*pbeta*dplabg(3)
0356 radgamp(3)=-pgamma*pbeta*EBrems+pgamma*dplabg(3)
0357 ! Copy dplabg(1) and (2) components to radgamp unchanged
0358 ! to avoid having to interface to the entire phiout
0359 ! common block in C++ code.
0360 do i = 1, 2
0361 radgamp(i) = dplabg(i)
0362 enddo
0363 else
0364 EBrems=0D0
0365 radgamE=0D0
0366 do i = 1, 3
0367 radgamp(i)=0D0
0368 enddo
0369 endif
0370
0371 ! N is from the pyjets common block
0372 ! Determine the number of tracks to print.
0373 ! If there was radiative emission there is one
0374 ! more track than from the pythia record.
0375 if(mcRadCor_EBrems.gt.0.) then
0376 nrtrack = N + 1
0377 else
0378 nrtrack = N
0379 endif
0380
0381 if((msti(1).ge.91).and.(msti(1).le.94)) then
0382 msti(16)=0
0383 endif
0384
0385 F2 = py6f2
0386 F1 = py6f1
0387 R = py6r
0388 sigma_rad = mcRadCor_Sigrad
0389 SigRadCor = mcRadCor_sigcor
0390
0391 ! Print ASCII output if requested
0392 if(printascii .ne. 0) then
0393 ! Write the event header
0394 write(asciilun, 32) 0, IEV, genevent, msti(1), msti(12),
0395 + msti(16), pari(34), msti(15), pari(33), pari(53),
0396 + VINT(309), VINT(307), trueX, trueW2, trueNu,
0397 + VINT(313), pari(14), pari(15), pari(16),
0398 + pari(18), pari(22), sngl(py6f2), sngl(py6f1),
0399 + py6r, mcRadCor_Sigrad, mcRadCor_sigcor, EBrems,
0400 + VINT(319), VINT(45), nrtrack
0401 32 format((I4,1x),(I10,1x),3(I4,1x),(I10,1x),f9.6,1x,
0402 + I12,1x,
0403 + 2(f12.6,1x),7(f18.11,3x),12(f19.9,3x),I12,$/)
0404 write(asciilun, *)'============================================'
0405 ! Write track records
0406 ! N is from the pyjets common block
0407 do i = 1, N
0408 ! Check that the parent index K(I,3) of this particle is in
0409 ! a sensible range.
0410 if(K(I,3).le.nrtrack) then
0411 write(asciilun,34) I,K(I,1),K(I,2),K(I,3),K(I,4),K(I,5),
0412 + P(I,1),P(I,2),P(I,3),P(I,4),P(I,5),
0413 + V(I,1),V(I,2),V(I,3)
0414 endif
0415 enddo
0416 ! Write radiative photon
0417 if(mcRadCor_EBrems.gt.0.) then
0418 write(asciilun,34) nrtrack, 55, 22, 1, 0, 0,
0419 + sngl(dplabg(1)),sngl(dplabg(2)),sngl(-radgamp),
0420 + sngl(radgamE), 0., 0., 0., 0.
0421 endif
0422 34 format(2(I6,1x),I10,1x,3(I8,1x),8(f15.6,1x),$/)
0423 write(asciilun,*)'=============== Event finished ==============='
0424 endif ! printascii
0425 lastgenevent=NGEN(0,3)
0426
0427 ! Increment event count and return
0428 IEV = IEV + 1
0429 generate = 0
0430
0431 end function generate
0432 ! ================================================================
0433
0434
0435 ! ================================================================
0436 ! Perform end-of-run actions.
0437 ! Print output to screen.
0438 ! Close ASCII file if it was requested.
0439 ! ================================================================
0440 subroutine finish()
0441
0442 include 'pythia.inc'
0443
0444 ! Print cross sections.
0445 call PYSTAT(1)
0446 call PYSTAT(4)
0447
0448 write(*, *)"The charm mass used is: ", PMAS(4,1)
0449
0450 ! Print the Pythia cross section which is needed to get an absolut
0451 ! normalisation the number is in microbarns
0452 write(*, *)'==================================================='
0453 write(*, *)'Pythia total cross section normalisation:',
0454 + pari(1)*1000, ' microbarn'
0455 write(*, *)'Total Number of generated events', MSTI(5)
0456 write(*, *)'Total Number of trials', NGEN(0,3)
0457 write(*, *)'==================================================='
0458
0459 if(printascii .ne. 0) then
0460 close(asciilun)
0461 endif
0462
0463 ! Check pdf status
0464 call PDFSTA
0465
0466 end subroutine finish
0467 ! ================================================================
0468
0469
0470 ! ================================================================
0471 ! Reset publicly accessible values to default values.
0472 ! ================================================================
0473 subroutine reset()
0474 genevent = 0
0475 trueX = 0.D0
0476 trueW2 = 0.D0
0477 trueNu = 0.D0
0478 F2 = 0.D0
0479 F1 = 0.D0
0480 R = 0.D0
0481 sigma_rad = 0.D0
0482 SigRadCor = 0.D0
0483 EBrems = 0.D0
0484 radgamE = 0.D0
0485 radgamp = 0.D0
0486 end subroutine reset
0487 ! ================================================================
0488
0489 end module pythia6
0490 ! ======================================================================
