bin/BCVEGPY/bcvegpy.f

0001 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0002 c                            main idea                                 c

0003 c  using vegas together with pythia to get the more precise results.   c

0004 c  by the several first runs of in the vegas, we may get the optimized c

0005 c  density distribution function which make the distribution of cross- c

0006 c  section not fluctuate too much. After the running of vegas, it will c

0007 c  call the pythia subroutines to generate events. in this way the     c

0008 c  MC efficiency is greatly improved.                                  c

0009 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0010 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0011 c    A pure linux version BCVEGPY2.1;    using GNU C compiler make     c

0012 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0013 c!!!       to save all the datas, you have to create a directory    !!!c

0014 c!!!              named ( data ) at the present directory.   !!!    !!!c

0015 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0016 c main program for seting the external process gg->bc+b+\bar{c}        c

0017 c into pythia, whose version is pythia6208.if higher version of pythia c

0018 c exists, your may directly replay the old one with the new one, but   c

0019 c remember to commoment out two subroutines upinit and upevnt there.   c

0020 c we also can provide the quark-anti-quark annihilation mechanism      c

0021 c which via the subprocess q+\bar{q}->Bc+b+\bar{c} with Bc in s-wave   c

0022 c states. Such mechanism has a quite small contribution around 1% of   c

0023 c that of gluon-gluon fusion.                                          c

0024 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0025 c  bcvegpy1.0   finished  in 10, auguest 2003                          c

0026 c               improved  in 1, november 2003                          c

0027 c copyright (c) c.h. chang, chafik driouich, paula eerola and x.g. wu  c

0028 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0029 c the original version of bcvegpy is in reference: hep-ph/0309120      c

0030 c or, in computer physics communication 159, 192-224(2004).            c

0031 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0032 c second version bcvegpy2.0, in which (p-wave) generation based on     c

0033 c the gluon-gluon fusion subprocess is given. the gauge invariance     c

0034 c for all the p-wave states have been checked exactly. because we      c

0035 c have implicitly used the propertities of poralization vector to      c

0036 c simplify the amplitude, (such (p.\epsilon(p)=0)) the gauge check can c

0037 c not be done by using the present amplitudes. the interesting reader  c

0038 c may ask the authors for the (full amplitudes) that keep the gauge    c

0039 c invariance exactly.         reference: hep-ph/0504017                c

0040 c      computer physics communication 174, 41,251(2006)                c

0041 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0042 c    bcvegpy2.0         finished in 24, febrary 2005                   c

0043 c   copyright (c)       c.h chang, j.x. wang and x.g. wu               c

0044 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0045 c problems or suggestions email to:         wuxg@itp.ac.cn             c

0046 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0047 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0048 c this is the linux version for BCVEGPY2.1, with better modularity and c

0049 c code reusability, i.e. less cross-talk among different modules.      c

0050 c this version combines all the feedback suggestions from the users.   c

0051 c thanks are given for: y.n. gao, j.b. he and z.w. yang                c

0052 c                        finished in 6, April 2006                     c

0053 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0054
0055 c...main program.

0056 c      program bcvegpy

0057       subroutine bcvegpy
0058       implicit double precision(a-h, o-z)
0059         implicit integer(i-n)
0060
0061 c...three pythia functions return integers, so need declaring.

0062       external pydata
0063
0064 c...pythia common block.

0065       common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
0066       parameter (maxnup=500)
0067       common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,idup(maxnup),
0068      &istup(maxnup),mothup(2,maxnup),icolup(2,maxnup),pup(5,maxnup),
0069      &vtimup(maxnup),spinup(maxnup)
0070       save /hepeup/
0071
0072       parameter (maxpup=100)
0073       integer pdfgup,pdfsup,lprup
0074       common/heprup/idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
0075      &idwtup,nprup,xsecup(maxpup),xerrup(maxpup),xmaxup(maxpup),
0076      &lprup(maxpup)
0077       save /heprup/
0078
0079 c...user process event common block.

0080       common/pypars/mstp(200),parp(200),msti(200),pari(200)
0081         common/counter/ibcstate,nev
0082         common/vegcross/vegsec,vegerr,iveggrade
0083         common/confine/ptcut,etacut
0084         logical generate
0085         common/genefull/generate
0086       common/totcross/appcross
0087       common/ptpass/ptmin,ptmax,crossmax,etamin,etamax,
0088      &  smin,smax,ymin,ymax,psetamin,psetamax
0089       common/loggrade/ievntdis,igenerate,ivegasopen,igrade
0090         common/mixevnt/xbcsec(8),imix,imixtype
0091         character*8 begin_time,end_time,blank
0092
0093
0094 c....Temporaty file for initialization/event output

0095         MSTP(161) = 77
0096         OPEN (77, FILE='BCVEGPY.init',STATUS='unknown')
0097         MSTP(162) = 78
0098         OPEN (78, FILE='BCVEGPY.evnt',STATUS='unknown')
0099 c....Final Les Houches Event file, obtained by combaining above two

0100         MSTP(163) = 79
0101         OPEN (79, FILE='BCVEGPY.lhe',STATUS='unknown')
0102         MSTP(164) = 1
0103
0104 c*********************************************

0105 c... setting initial parameters in parameter.F.

0106 c... User may change the values to his/her favorate one.

0107 c*********************************************

0108       call setparameter
0109
0110 c... initialization, including:

0111 c... if ivegasopen=1, then generate vegas-grade;

0112 c... open files to record data (intermediate results for vegas or

0113 c... pythia running information);

0114 c... if ivegasopen=0 and igrade=1, then initialize the grade;

0115 c... setting some initial values for pythia running.

0116         call evntinit
0117
0118 C...Fills the HEPRUP commonblock with info on incoming beams and allowed

0119 C...processes, and optionally stores that information on file.

0120         call bcvegpy_PYUPIN
0121 c...

0122 c*****************************************************************

0123 c...using pybook to record the different distributions or

0124 c...event-distributions.

0125 c*****************************************************************

0126 c*****************************************************************

0127 c...any subroutines about pybook can be conveniently commented out

0128 c...and then one can directly use his/her own way to record the data.

0129 c*****************************************************************

0130
0131 c...pybook init.

0132         call pybookinit
0133
0134 c...approximate total cross-section.

0135         appcross =0.0d0
0136 c

0137         blank  ='    '
0138         ncount =0
0139 c       call time(begin_time)

0140
0141 c*******************************************************

0142 c...there list some typical ways for recording the data.

0143 c...users may use one convenient way/or their own way.

0144 c******************************************************

0145 c***********find highest weight*******************

0146       MAXWGT=0.0D0
0147       do iev=1,10000
0148       call pyevnt
0149       if (MAXWGT.le.xwgtup) then
0150          MAXWGT=xwgtup
0151       end if
0152       end do
0153 c---------------------------------------------------

0154       do iev=1,nev
0155                 call pyevnt
0156 c     call bcvegpy_write_lhe

0157 c**************The unweighting scheme****************

0158 c******the ratio of xwgtup and the xwgtup_max setted at the beginning******

0159       xwt_r=xwgtup/MAXWGT
0160
0161       x_r = rand(0)
0162
0163       if (x_r.le.xwt_r.and.xwt_r.le.1) then
0164 c         write(*,*) "random number ",x_r," ratio ",xwt_r

0165 c         write(*,*) "max xwgtup",MAXWGT

0166          call bcvegpy_write_lhe
0167       end if
0168       if(xwt_r.gt.1) then
0169          do ic=1,int(xwt_r)
0170 c         write(*,*) "fill times",ic

0171          call bcvegpy_write_lhe
0172          end do
0173          sxwt=xwt_r-int(xwt_r)
0174          x_r = rand(0)
0175          if (x_r.le.sxwt) then
0176          call bcvegpy_write_lhe
0177          end if
0178       end if
0179
0180 c****************************************************************

0181                 if (idwtup.eq.1.and.iev.ne.1.and.generate) then
0182                 call pylist(7)
0183 c               call time(end_time)

0184                 print *, iev,blank,end_time
0185             end if
0186
0187                 if(msti(51).eq.1) go to 400
0188
0189 c*********************************************************

0190             do i=1,10
0191              if(k(i,2).eq.541) then
0192 c...pt of bc.

0193                 pt =pyp(i,10)
0194 c...true rapidity.

0195                   eta=pyp(i,17)
0196 c...pseudo-rapidity

0197                   pseta=pyp(i,19)
0198 c...these two constrain (and other) may be added here to partly compensate

0199 c...some numerical problems.

0200                           if(pt.lt.ptcut) xwgtup=0.0d0
0201                   if(abs(eta) .gt. etacut) xwgtup=0.0d0
0202 c...rapidity of the hard-interaction subsystem(ln(x1/x2)/2.0)

0203                           y=pari(37)
0204 c...the mass of the complete three- final state(\sqrt(shat))

0205                           st=pari(19)
0206
0207 c**********************************************************************

0208 c...   users may use his own way to record the data      ..............

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

0210 c...to fill the histogram. we list three methods to get the differential

0211 c...distributions. 1) idwtup=3; 2) (idwtup=1 and generate.eq.fause); 3)

0212 c...(idwtup=1 and generate.eq.true). where method 1) and 2) are the

0213 c...quickest, while the third method is slow.

0214
0215 c...we also list three ways to get the event number distributions:

0216 c...1) (idwtup=1 and generte.eq.true) and ievntdis.eq.1;

0217 c...2)  idwtup=3 and ievntdis.eq.0; 3) (idwtup=1 and generate.eq.fause)

0218 c...and ievntdis.eq.0; the method 2) and 3) are the same, both needs

0219 c...a proper normalization for numbers and at the same time

0220 c...recording every event with it's corresponding weight so as to get

0221 c...final right event number distributions. the method 1) is general

0222 c...one used by experimental, which will spend a long time. so for

0223 c...theoretical studies we suggest using method 2) or 3).

0224 c**********************************************************************

0225         call uppyfill(idwtup,generate,xwgtup,pt,eta,st,y,pseta)
0226
0227                           isub=msti(1)
0228                   ncount=ncount+1
0229                   if(ncount.le.10) then
0230                   write(*,'(a,i5)') 'following event is subprocess',isub
0231                   call pylist(7)
0232                         call pylist(1)
0233                 end if
0234                 call pyedit(2)
0235                end if
0236             end do
0237         end do
0238
0239
0240 c***************************************************************

0241 c...close grade files

0242       call upclosegradefile(iveggrade,imix,imixtype)
0243 c***************************************************************

0244
0245 c***************************************************************

0246 c...can be commentted out by user, if using his/her own way to

0247 c...record the data.

0248 c---------------------------------------------------------------

0249 c...pyfact all the pybooks.

0250 c       call uppyfact(idwtup,generate,ievntdis)

0251 c...open files to record the obtained pybook data for distributions.

0252 c       call updatafile

0253 c...dump the data into the corresponding files.

0254 c      call uppydump

0255 c...close pybook files.

0256 c      call upclosepyfile

0257 c***************************************************************

0258
0259 c       call time(end_time)

0260         write(3,'(a,d19.6,a)') 'maximum diff. cross-sec=',crossmax,'pb'
0261         write(3,'(i9,3x,a,3x,a)') nev,begin_time,end_time
0262
0263 c...when the number of sampling points are high enough, it is

0264 c...just the real value. for (idwtup.eq.1.and.ievntdis.eq.1),

0265 c...becaue of small event number the value of appcros is not

0266 c...accurate.

0267         if(idwtup.ne.1.and.ievntdis.eq.1) then
0268           write(3,'(a,d16.6,1x,a)') '!approxi. total cross-sec=',
0269      &  appcross,'nb'
0270       end if
0271
0272 c...store the approximate total cross-section.

0273 c...appcross=\sum(xwgtup*wgt)/nev. when the number of sampling points

0274 c...are high enough, it will be the true value obtained from pythia.

0275         pari(1) =appcross*1.0d-3
0276         write(3,'(a,3x,d16.6,x,a)') 'total cross-sec(from pythia)=',
0277      &  pari(1),'mb'
0278
0279 c*****************************

0280 c...histograms.

0281        call pyhist
0282 c*****************************

0283 400     continue
0284 c...type cross-section table.

0285 c      call pystat(1)

0286 c

0287 c      write(*,'(a)') 'the cross-section in PYSTAT table is nonsense'

0288 c      write(*,'(a)') 'see the true value in obtained cross-section file'

0289 c      write(*,*)

0290 c      write(3,*)

0291 c      write(3,'(a)') 'note: cross-section in PYSTAT table is nonsense,'

0292 c      write(3,'(a)') 'especially for the mixed events.'

0293 c      write(3,*)

0294        write(*,*)
0295        write(*,'(a)') '     program finished !'
0296        write(*,*)
0297
0298 c***************************

0299
0300 c***************************

0301 c...close intemediate file.

0302       close(3)
0303
0304 C....Produce final LHE file

0305 c      CALL PYUPIN

0306 c      call bcvegpy_PYUPIN

0307       CALL PYLHEF
0308
0309       end     !end of main program