bin/BCVEGPY/parameter.F

0001 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0002 c...this subroutine is used to set the necessary parameters for      c

0003 c...the initialization for hadronic production of bc meson.          c

0004 c...to use the program youd need to make a directory: (data) to      c

0005 c...save all the obtained data-files.                                c

0006 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0007 c   to have a better understanding of setting the parameters         c

0008 c   you may see the README file to get more detailed information.    c

0009 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0010 c copyright (c) z.x zhang, chafik driouich, paula eerola and x.g. wu c

0011 c reference: hep-ph/0309120                                          c

0012 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0013
0014         subroutine setparameter
0015 c...preamble: declarations.

0016         implicit double precision(a-h, o-z)
0017         implicit integer(i-n)
0018
0019 #include "bcvegpy_set_par.inc"
0020
0021 c...user process event common block.

0022       common/pypars/mstp(200),parp(200),msti(200),pari(200)
0023       common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
0024       common/pydatr/mrpy(6),rrpy(100)
0025
0026 c...transform of running information.

0027         common/confine/ptcut,etacut
0028 c...parameters transformtion.

0029         double complex colmat,bundamp
0030       common/upcom/ecm,pmbc,pmb,pmc,fbcc,pmomup(5,8),
0031      &  colmat(10,64),bundamp(4),pmomzero(5,8)
0032       common/funtrans/nq2,npdfu
0033         common/rconst/pi
0034         common/usertran/ishower,idpp
0035         common/vegcross/vegsec,vegerr,iveggrade
0036 c...to transform the subprocess cross-section.

0037       common/subopen/subfactor,subenergy,isubonly
0038 c...transform of pdf information

0039         common/outpdf/ioutpdf,ipdfnum
0040 c...transform the subprocess information

0041       common/qqbar/iqqbar,iqcode
0042 c...transform of the vegas information

0043       common/vegasinf/number,nitmx
0044 c...transform the events number and bc state.

0045         common/counter/ibcstate,nev
0046 c...transform some variables

0047       common/loggrade/ievntdis,igenerate,ivegasopen,igrade
0048 c...ioctet--whether getting the color-octet component contributions.

0049 c...here only for gg->(c\bar{b})+b+~c, (c\bar{b}) in color-octet

0050 c...s-waee states. coeoct--coefficient for color-octet

0051 c...matrixment: relation between the color-octet matrixment and the

0052 c...color singlet matrixment.

0053       common/coloct/ioctet
0054         common/octmatrix/coeoct
0055 c...transform the first derivative of the wavefunction at the zero or the

0056 c...wavefunction at the zero.

0057         common/wavezero/fbc
0058 c...xbcsec(8) records the total differential cross-sections for different

0059 c...states: 1---singlet 1s0; 2---singlet 3s1; 7---octet 1s0; 8---octet 3s1;

0060 c...3---singlet 1p1; 4---singlet 3p0; 5---singlet 3p1; 6---singlet 3p2.

0061         common/mixevnt/xbcsec(8),imix,imixtype
0062
0063         logical wronginput
0064 c

0065 c... Read parameter setting from file bcvegpy_set_par.nam

0066 c

0067         Call Read_parameter_settings
0068 c

0069 c... change the intitial state of the random number

0070 c        mrpy(1) = 19780503             ! default value

0071          mrpy(1) = ibcrandom
0072          write( 6, * ) ' '
0073          write( 6, * ) ' Change default value of random, mrpy(1), to ',
0074      +   mrpy(1)
0075          write( 6, * ) ' '
0076 c... end random change

0077
0078         pi = dacos(-1.0d0)
0079
0080 c*************************************************

0081 c** some commen parameters that are used for both

0082 c** imix=0 or imix=1 are listed here for convenience.

0083 c*************************************************

0084
0085 c...setting the basic parameters. note pmbc should be equal to

0086 c...(pmb+pmc) to ensure the gauge invariance of the amplitude.

0087 c      pmb  =4.9d0    !b quark mass (gev)

0088 c      pmc  =1.386d0    !~c quark mass (gev)

0089       pmb  = pmassb    !b quark mass (gev)
0090       pmc  = pmassc    !~c quark mass (gev)
0091       pmbc = pmb+pmc   !~bc mass (gev)
0092
0093 c------------------------------------------

0094 c...kinematical variables

0095       ptcut =0.0d+0   !bc p_t cut (gev)
0096       etacut=1.0d+9   !bc rapidity (y) cut
0097
0098 c---------------------------------------

0099 c...npdfu=2:lhc=1.40d+4GeV; npdfu=1:tevatron=1.96d+3GeV.

0100 c---------------------------------------

0101
0102 c...hadron information

0103       nq2   =3      !types of q^2 (seven typical energy scales are
0104 c                     ! designed in program)

0105 c      npdfu =2      !types of incoming beams. =1, tevatron; =2, lhc.

0106       npdfu = naccel
0107       if(npdfu.eq.1) ecm=ENERGYOFTEVA  !defined in run.F
0108       if(npdfu.eq.2) ecm=ENERGYOFLHC   !defined in run.F
0109       mstp(51) =7   !types of parton distribution functions. used for
0110 c                     !pythia inner pdf and under the condition ioutpdf=0.

0111       idpp     =3      !=idwtup: types of pythia generates events
0112       mstu(111)=1      !setting the value of alphas. pythia parameter
0113       paru(111)=0.2d0  !constant value of alphas, used when mstu(111)=0.
0114       ievntdis =0      !switch on/off to get the event number
0115 c                        !distributions for p_t, rap, shat, pseta, y, z.

0116
0117 c------------------------------------------

0118 c...event information

0119       igenerate =0 !whether generating full events must be with idwtup=1
0120 c      ishower=0    !switch off aspects: initial and final state showers,

0121 c                    !multiple interactions and hadronization.

0122       ishower = i_shower
0123
0124 c...using pythia pdf or not.

0125       ioutpdf  =1
0126 c...=100, grv98lo; =200, mrst2001l; =300, cteq6l1

0127       ipdfnum  =300
0128       if(ioutpdf.eq.1 .and. ipdfnum.eq.300) call setctq6(4)   !cteq6l1
0129
0130 c**************************************************

0131 c**************************************************

0132 c...whether to get the mixed events for bc meson.

0133 c...here, we only consider the mixed results for

0134 c...the dominant gluon-gluon fusion mechanism.

0135 c--------------------------------------------------

0136 c...if(imix==1), then we can generate the mixed events

0137 c...up to the eight states (1s0,3s1,1p1,3p0,3p1,3p2,1s0_8,3s1_8),

0138 c...which are produced through the gluon-gluon fusion mechanism.

0139 c...if(imix==0), then one can generate every states separately.

0140 c**************************************************

0141 c       imix=1

0142         imix = i_mix
0143
0144 c--------------------------------------------------

0145 c***** set the parameters for mixing results

0146 c--------------------------------------------------

0147
0148         if (imix.eq.1) then
0149 c*************************************************

0150 c...setting fixed values for mixing processes.

0151 c...imixtype: set how many states need to be mixed:

0152 c...imixtype=1: all the eight states need to be mixed.

0153 c...imixtype=2: the mixed events for two color-singlet s-wave states.

0154 c...imixtype=3: the mixed events for the four color-singlet

0155 c...            p-wave plus two color-octet s-wave states,

0156 c************************************************

0157         imixtype=1
0158
0159 c*************************************************

0160 c...some of the parameters are fixed only for simplicity.

0161 c...here, to make our program short, we only provide

0162 c...three possible way to get the mixed results are programmed:

0163 c... 1) ivegasopen=1, no matter what value of igrade.

0164 c... 2) ivegasopen=0 and igrade=1.

0165 c... it is not an appreciable way:

0166 c... 3) ivegasopen=0 and igrade=0.

0167 c*************************************************

0168         mixmethod=1
0169
0170         if(mixmethod.eq.1) then
0171              ivegasopen=1
0172            igrade    =0   !or 1, no use.
0173         end if
0174         if(mixmethod.eq.2) then
0175            ivegasopen=0
0176            igrade    =1
0177         end if
0178         if(mixmethod.eq.3) then
0179            ivegasopen=0
0180            igrade    =0
0181         end if
0182
0183 c----------------------------

0184 c... the following value should be increased

0185 c... to achieve the users precision goal.

0186 c----------------------------

0187           number=NVEGCALL
0188           nitmx =NVEGITMX
0189           nev   =NUMOFEVENTS
0190         igrade   =1   ! no use when ivegasopen=1
0191         iveggrade=0   ! should be fixed to 0
0192
0193 c---------------------------------------------------

0194 c-- fixed parameter in the case of imix=1

0195         isubonly =0
0196           subenergy=100.0d0  ! no use for present purpose
0197         iqqbar   =0        ! should be fixed to 0.
0198         iqcode   =2        ! no use for present purpose
0199
0200           return
0201         end if
0202 c--------------------------------------------------

0203 c***** end of the initialization for mixing results

0204 c--------------------------------------------------

0205
0206 c******************************************************

0207 c******************************************************

0208 c***** The following is for generating single state

0209 c******************************************************

0210 c******************************************************

0211
0212 c*****************************************

0213 c*** this is the previous way for easy running. with ireaddata=1,

0214 c*** it will read the all the necessary parameter values from the

0215 c*** datafile (totput.dat).

0216 c***

0217 c*** since now we take ( makefile ) to do the compiling,

0218 c*** such method is not very useful and now is commented out.

0219 c***

0220 c***    ireaddata=0

0221 c***    if(ireaddata.eq.1) goto 1001

0222 c*****************************************

0223 c*****************************************

0224 c...note: the following initial value is only for checking that

0225 c...the program is ok. to do numerical calculations/event generation,

0226 c...one needs to be careful about all these values, especailly

0227 c...those parameters that can improve the precision.

0228 c*****************************************

0229
0230 c****************************************

0231 c...to choose what bc state we need to generate.

0232 c...=1, color-singlet 1s0 ;=2, color-singlet 3s1 ;=3, 1p1 ;=4, 3p0 ;

0233 c...=5, 3p1 ;=6, 3p2 ;=7, color-octet 1s0 ;=8, color-octet 3s1.

0234         ibcstate  =4
0235
0236 c*****************************************

0237 c...for s-wave, r(0)=1.241  --->  \psi(0)=(0.35).

0238       ioctet=0
0239       if(ibcstate.eq.7 .or. ibcstate.eq.8) then
0240         ioctet=1
0241         coeoct=0.2d0
0242         if(ibcstate.eq.7) ibcstate=1
0243         if(ibcstate.eq.8) ibcstate=2
0244       end if
0245
0246 c******************************************

0247         if(ibcstate.eq.1 .or. ibcstate.eq.2) then
0248           fbc =1.241
0249           fbcc=dsqrt(3.0d0*fbc**2/pi/pmbc)
0250         else
0251 c...for p-wave, r'(0)=0.44833  --->  \psi'(0)=(0.219); r'(0)**2=(0.201)

0252           fbc =0.44833
0253 c...the value of p-wave matrix element <0|p|0>.

0254           fbcc=fbc**2*9.0d0/(2.0d0*pi)
0255         end if
0256
0257 c****************************************

0258 c...vegas parameters

0259       isingmethod=1
0260       if(isingmethod.eq.1) then
0261           ivegasopen=1   !whether using vegas.
0262           igrade    =0   !or 1, no use.
0263       end if
0264       if(isingmethod.eq.2) then
0265           ivegasopen=0
0266           igrade    =1   ! whether using existed grade,
0267       end if
0268       if(isingmethod.eq.3) then
0269           ivegasopen=0
0270           igrade    =0
0271       end if
0272
0273         number=NVEGCALL !total number of calling fxn in each iteration in vegas
0274         nitmx =NVEGITMX     !maximum iteration used in vegas
0275         nev   =NUMOFEVENTS   !number of events to be generated
0276
0277 c****************************************

0278 c...improve vegas grade

0279       iveggrade=0  ! switch on/off to use the existed grade before
0280 c                  ! running vegas, the initial condition for the running

0281 c                  ! and that for getting the existed grade should be the

0282 c                  ! same. this is used to get a more precise sampling

0283 c                  ! grade. used together with ivegasopen=1.

0284
0285 c***************************************

0286 c...subprocess information

0287         isubonly =0    !switch on whether only to get the information,
0288 c                      !for the subprocess including cross-section, pt,

0289 c                      !rapidity distributions and so on.

0290 c...c.m. energy(gev) of the subprocess gg->bc+b+c~

0291       if(isubonly.eq.1) subenergy=100.0d0
0292
0293 c***************************************

0294 c...mechanism through subprocess q~q->bc

0295       iqqbar   =0      !switch on/off whether using the subprocess
0296 c                        q+\bar{q}-> to generate the bc events.

0297 c...all the quarks are massless.

0298 c...=1, u; =2, d; =3, s

0299       if(iqqbar.eq.1) iqcode=2
0300
0301 c  1001  continue

0302
0303 c****************************************

0304 c****************************************

0305 c****************************************

0306 c**** comment out the method of reading data from outer

0307 c**** data file.

0308 c**** this way is the old way to save compile time and is

0309 c**** not necessary for the present linux version, since

0310 c**** now we use ( make ) to compile the program.

0311 c****************************************

0312 c...open input data file.

0313 c       if(ireaddata.eq.1) then

0314 c       open(unit=1,file='totput.dat',status='old',access='sequential',

0315 c     &    form='formatted')

0316 c...read the (zero point wave function) -----> initial for fbcc.

0317 c       read(1,*) pmbc,pmb,pmc,fbcc

0318 c*** here the input ibcstate should be 1,2,3,4,5,6.

0319 c       read(1,*) ptcut,etacut,ecm,ibcstate,igenerate,ivegasopen

0320 c       read(1,*) number,nitmx

0321 c       read(1,*) nq2,npdfu,nev

0322 c       read(1,*) ishower,mstp51,idpp,mstu111,paru111

0323 c       read(1,*) isubonly,subenergy,igrade

0324 c       read(1,*) inumevnt,iveggrade,iqqbar,iqcode

0325 c       read(1,*) ioutpdf,ipdfnum,ioctet,coeoct

0326 c

0327 c****************************************

0328 c         if(ibcstate.eq.7 .or. ibcstate.eq.8) then

0329 c          if(ibcstate.eq.7) ibcstate=1

0330 c          if(ibcstate.eq.8) ibcstate=2

0331 c          ioctet=1

0332 c          coeoct=0.2d0

0333 c         end if

0334
0335 c****************************************

0336 c        if(ibcstate.eq.1 .or. ibcstate.eq.2) then

0337 c...decay constant f_{bc} of s-wave.

0338 c          fbc=fbcc

0339 c          fbcc=dsqrt(fbcc**2*3.0d0/(pmbc*pi))

0340 c        else

0341 c...the value of p-wave matrix element <0|p|0>.

0342 c          fbc=fbcc

0343 c          fbcc=fbcc**2*9.0d0/(2.0d0*pi)

0344 c        end if

0345 c

0346 c        mstp(51)=mstp51

0347 c        mstu(111)=mstu111

0348 c        paru(111)=paru111

0349 c

0350 c        close(1)

0351 c       end if

0352 c***************************************

0353 c***************************************

0354 c***************************************

0355
0356 c----------------------------------------------

0357 c...error message.

0358       wronginput=.false.
0359         call uperror(wronginput)
0360         if(wronginput) stop '-----input error! stop the program !!!'
0361
0362         call parameters()
0363       call dparameters()
0364       call coupling()
0365
0366         return
0367         end
0368 c***************************************

0369 c***************************************

0370         Subroutine Read_parameter_settings
0371 c

0372 c... Get parameters from namelist

0373         implicit double precision(a-h, o-z)
0374         implicit integer(i-n)
0375
0376         Namelist / bcvegpy_set_par / ENERGYOFTEVA, ENERGYOFLHC,
0377      +   pmassb,  pmassc, naccel, i_shower, i_mix,
0378      +   NUMOFEVENTS, NVEGCALL, NVEGITMX, ibcrandom
0379 #include "bcvegpy_set_par.inc"
0380 c

0381 c-------------------------------------------------------------------------------

0382 c

0383         open( unit=1, file='bcvegpy_set_par.nam',Status='Old',Err=99)
0384         read( 1, nml=bcvegpy_set_par, err=90)
0385         write( 6, * ) ' '
0386         write( 6, * ) ' Contents of namelist *bcvegpy_set_par*: '
0387         write( 6, nml=bcvegpy_set_par)
0388         write( 6, * ) ' '
0389         Close( 1 )
0390         Return
0391 c

0392   90    Write( 6, * ) ' !!!!! Unable to read namelist bcvegpy_set_par '
0393         Call Exit
0394   99    Write( 6, * ) ' !!!!! Unable to open bcvegpy_set_par.nam'
0395         Call Exit
0396         End