Writing raw event data during FO integration

Hi,

I'd like to write raw event data to a text file during a MadGraph NLO FO run (4-vectors, ibody values, and weights). My goal is to read these values into root and construct my own histograms from them. Currently, I'm doing this in the analysis_fill subroutine of the analysis_root_xxx.f (where "xxx" is a user-defined name). This is a file that the user can use to define custom histogram filling. This generates four text files (one for each subprocess) which I then parse and combine in root to create histograms. The histograms I've generated are similar but not the same as those generated by the analysis_root_xxx.f file.

One issue I noticed is that the first chunk of events I was writing to my text files had weights of 0 for weights(i) where i=2:9. These events were not being used to fill the fortran-generated histograms. I corrected for this, but I'm still getting slightly different results.

My question is: do you know where, in the code, should I place a subroutine for writing the raw event data into a text file so that I can build my own histograms in root from that data and successfully reproduce the histograms given by the fortran code?

I've attached the file analysis_root_aaa.f which is the way I'm currently writing the generated events to text files.

Thanks in advance!

Sam

My parameter settings look like this:

#NLOSM_FO model lite
generate p p > h mu+ mu- [QCD]
output NLOSM_ulite
y
launch
fixed_order=ON
set ebeam1 7000.0
set ebeam2 7000.0
set pdlabel cteq6l1
set ptj 0.0
set etaj 5.0
set ptl 10.0
set etal 2.5
set drll 0.0
set mll 0.0
set ptgmin 10.0
set R0gamma 0.1
set mt 1.725000e+02
set mz 9.118760e+01
set mh 1.250900e+02
set mu- : 0.105658367
set gf 1.166370e-05
set wz 2.495200
set ww 2.088720
set wh 4.08e-03
done

I changed these settings to decrease computing time so that I can debug more quickly (taken from run_card.dat):

#***********************************************************************
# Number of points per itegration channel (ignored for aMC@NLO runs) *
#***********************************************************************
-1 = req_acc_FO ! Required accuracy (-1=ignored, and use the
                           ! number of points and iter. below)
# These numbers are ignored except if req_acc_FO is equal to -1
 50 = npoints_FO_grid ! number of points to setup grids
 1 = niters_FO_grid ! number of iter. to setup grids
 100 = npoints_FO ! number of points to compute Xsec
 1 = niters_FO ! number of iter. to compute Xsec

c
c Example analysis for "p p > h mu+ mu- [QCD]" process
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine analysis_begin(nwgt,weights_info)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer nwgt
      character*(*) weights_info(*)
      integer j,kk,l,nwgt_analysis
      common/c_analysis/nwgt_analysis
      character*6 cc(2)
      data cc/'|T@NLO','|T@LO '/
      real * 8 etaMin,etaMax,etaBin,eMin,eMax,eBin,pi
      real * 8 ptMin,ptMax,ptBin,phiMin,phiMax,phiBin
      parameter (pi=3.14159265358979312d0)
      call open_root_file()
      nwgt_analysis=nwgt

      etaBin=1.d-1
      etaMin=-5.d0
      etaMax=5.d0

      eBin=1.d1
      eMin=0.d0
      eMax=1.d3

      ptBin=2.d0
      ptMin=0.d0
      ptMax=2.5d2

      phiBin=pi/10.d0
      phiMin=-1.d0*pi
      phiMax=1.d0*pi

      do j=1,2
      do kk=1,nwgt_analysis
      l=(kk-1)*100+(j-1)*28

      call rbook(l+ 1,'eMu1 '//cc(j)//weights_info(kk),
     & eBin,eMin,eMax)
      call rbook(l+ 2,'ptMu1 '//cc(j)//weights_info(kk),
     & ptBin,ptMin,ptMax)
      call rbook(l+ 3,'etaMu1 '//cc(j)//weights_info(kk),
     & etaBin,etaMin,etaMax)
      call rbook(l+ 4,'phiMu1 '//cc(j)//weights_info(kk),
     & phiBin,phiMin,phiMax)

      call rbook(l+ 5,'eMu2 '//cc(j)//weights_info(kk),
     & eBin,eMin,eMax)
      call rbook(l+ 6,'ptMu2 '//cc(j)//weights_info(kk),
     & ptBin,ptMin,ptMax)
      call rbook(l+ 7,'etaMu2 '//cc(j)//weights_info(kk),
     & etaBin,etaMin,etaMax)
      call rbook(l+ 8,'phiMu2 '//cc(j)//weights_info(kk),
     & phiBin,phiMin,phiMax)

      call rbook(l+ 9,'eZ '//cc(j)//weights_info(kk),
     & eBin,eMin,eMax)
      call rbook(l+ 10,'ptZ '//cc(j)//weights_info(kk),
     & ptBin,ptMin,ptMax)
      call rbook(l+ 11,'yZ '//cc(j)//weights_info(kk),
     & etaBin,etaMin,etaMax)
      call rbook(l+ 12,'etaZ '//cc(j)//weights_info(kk),
     & etaBin,etaMin,etaMax)
      call rbook(l+ 13,'phiZ '//cc(j)//weights_info(kk),
     & phiBin,phiMin,phiMax)
      call rbook(l+ 14,'mZ '//cc(j)//weights_info(kk),
     & 1.d0,5.d1,1.3d2)

      call rbook(l+ 15,'eH '//cc(j)//weights_info(kk),
     & eBin,eMin,eMax)
      call rbook(l+ 16,'ptH '//cc(j)//weights_info(kk),
     & ptBin,ptMin,ptMax)
      call rbook(l+ 17,'yH '//cc(j)//weights_info(kk),
     & etaBin,etaMin,etaMax)
      call rbook(l+ 18,'etaH '//cc(j)//weights_info(kk),
     & etaBin,etaMin,etaMax)
      call rbook(l+ 19,'phiH '//cc(j)//weights_info(kk),
     & phiBin,phiMin,phiMax)
      call rbook(l+ 20,'mH '//cc(j)//weights_info(kk),
     & 1.d0,1.05d2,1.45d2)

      call rbook(l+ 21,'eZH '//cc(j)//weights_info(kk),
     & eBin,eMin,eMax)
      call rbook(l+ 22,'ptZH '//cc(j)//weights_info(kk),
     & ptBin,ptMin,ptMax)
      call rbook(l+ 23,'yZH '//cc(j)//weights_info(kk),
     & etaBin,etaMin,etaMax)
      call rbook(l+ 24,'etaZH '//cc(j)//weights_info(kk),
     & etaBin,etaMin,etaMax)
      call rbook(l+ 25,'phiZH '//cc(j)//weights_info(kk),
     & phiBin,phiMin,phiMax)
      call rbook(l+ 26,'mZH '//cc(j)//weights_info(kk),
     & 5.d0,1.55d2,10.05d2)
      call rbook(l+ 27,'Mu1 nw '//cc(j)//weights_info(kk),
     & eBin,eMin,eMax)
      call rbook(l+ 28,'Mu1 nw2 '//cc(j)//weights_info(kk),
     & eBin,eMin,eMax)

      enddo
      enddo

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine analysis_end(xnorm)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      double precision xnorm
      integer i,jj
      integer kk,l,nwgt_analysis
      common/c_analysis/nwgt_analysis
      call wxnorm(xnorm)
c Do not touch the following lines. These lines make sure that the
c histograms will have the correct overall normalisation: cross section
c (in pb) per bin.
      do i=1,2
      do kk=1,nwgt_analysis
      l=(kk-1)*100+(i-1)*28
      do jj=1,26
        call ropera(l+jj,'+',l+jj,l+jj,xnorm,0.d0)
      enddo
      enddo
      enddo
      call close_root_file
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      include 'nexternal.inc'
      integer istatus(nexternal)
      integer iPDG(nexternal)
      double precision p(0:4,nexternal)
      double precision wgts(*)
      integer ibody
      double precision wgt,var
      integer i,kk,l,nwgt_analysis
      common/c_analysis/nwgt_analysis
      double precision www
      double precision p4Mu1(0:3),p4Mu2(0:3),p4Z(0:3),p4H(0:3),p4ZH(0:3)
      double precision eMu1,ptMu1,etaMu1,phiMu1
      double precision eMu2,ptMu2,etaMu2,phiMu2
      double precision eZ,ptZ,etaZ,yZ,phiZ,mZ
      double precision eH,ptH,etaH,yH,phiH,mH
      double precision eZH,ptZH,etaZH,yZH,phiZH,mZH
      double precision getrapidity,getpseudorap,dot
      external getrapidity,getpseudorap,dot
      if (nexternal.ne.6 .and. nexternal.ne.7) then
         write (*,*) 'error #1 in analysis_fill: '/
     & /'"p p > h mu+ mu- [QCD]": external particles'
         stop 1
      endif
      if (abs(ipdg(1)).gt.5 .and. abs(ipdg(1)).ne.21) then
         write (*,*) 'error #2 in analysis_fill: '/
     & /'"p p > h mu+ mu- [QCD]": incoming parton 1'
         stop 1
      endif
      if (abs(ipdg(2)).gt.5 .and. abs(ipdg(2)).ne.21) then
         write (*,*) 'error #3 in analysis_fill: '/
     & /'"p p > h mu+ mu- [QCD]": incoming parton 2'
         stop 1
      endif
      if (abs(ipdg(3)).ne.25) then
         write (*,*) 'error #4 in analysis_fill: '/
     & /'"p p > h mu+ mu- [QCD]": Higgs'
         stop 1
      endif
      if (abs(ipdg(4)).ne.13 .or. abs(ipdg(5)).ne.13) then
         write (*,*) 'error #5 in analysis_fill: '/
     & /'"p p > h mu+ mu- [QCD]": muons'
         stop 1
      endif
C
      do i=0,3
        p4H(i)=p(i,3)
        p4Mu1(i)=p(i,4)
        p4Mu2(i)=p(i,5)
        p4Z(i)=p4Mu1(i) + p4Mu2(i)
        p4ZH(i)=p4Z(i) + p4H(i)
      enddo

      eMu1=p4Mu1(0)
      ptMu1=sqrt(max(p4Mu1(1)**2+p4Mu1(2)**2,0d0))
      etaMu1=getpseudorap(p4Mu1(0),p4Mu1(1),p4Mu1(2),p4Mu1(3))
      phiMu1=atan2(p4Mu1(2),p4Mu1(1))

      eMu2=p4Mu2(0)
      ptMu2=sqrt(max(p4Mu2(1)**2+p4Mu2(2)**2,0d0))
      etaMu2=getpseudorap(p4Mu2(0),p4Mu2(1),p4Mu2(2),p4Mu2(3))
      phiMu2=atan2(p4Mu2(2),p4Mu2(1))

      eZ=p4Z(0)
      mZ=sqrt(max(dot(p4Z,p4Z),0d0))
      ptZ=sqrt(max(p4Z(1)**2+p4Z(2)**2,0d0))
      yZ=getrapidity(p4Z(0),p4Z(3))
      etaZ=getpseudorap(p4Z(0),p4Z(1),p4Z(2),p4Z(3))
      phiZ=atan2(p4Z(2),p4Z(1))

      eH=p4H(0)
      mH=sqrt(max(dot(p4H,p4H),0d0))
      ptH=sqrt(max(p4H(1)**2+p4H(2)**2,0d0))
      yH=getrapidity(p4H(0),p4H(3))
      etaH=getpseudorap(p4H(0),p4H(1),p4H(2),p4H(3))
      phiH=atan2(p4H(2),p4H(1))

      eZH=p4ZH(0)
      mZH=sqrt(max(dot(p4ZH,p4ZH),0d0))
      ptZH=sqrt(max(p4ZH(1)**2+p4ZH(2)**2,0d0))
      yZH=getrapidity(p4ZH(0),p4ZH(3))
      etaZH=getpseudorap(p4ZH(0),p4ZH(1),p4ZH(2),p4ZH(3))
      phiZH=atan2(p4ZH(2),p4ZH(1))
C
      call wweights(wgts, nwgt_analysis, ibody, p4Mu1, p4Mu2, p4H)
      call rfill(28,eMu1,1.d0)

c call wdebug(p4Mu1, p4Mu2, p4H, mZ, mH, mZH, wgts, nwgt)
      do i=1,2
         do kk=1,nwgt_analysis
            www=wgts(kk)
            l=(kk-1)*100+(i-1)*28
            if (ibody.ne.3 .and.i.eq.2) cycle
            call rfill(l+ 1,eMu1,WWW)
            call rfill(l+ 2,ptMu1,WWW)
            call rfill(l+ 3,etaMu1,WWW)
            call rfill(l+ 4,phiMu1,WWW)

            call rfill(l+ 5,eMu2,WWW)
            call rfill(l+ 6,ptMu2,WWW)
            call rfill(l+ 7,etaMu2,WWW)
            call rfill(l+ 8,phiMu2,WWW)

            call rfill(l+ 9,eZ,WWW)
            call rfill(l+ 10,ptZ,WWW)
            call rfill(l+ 11,yZ,WWW)
            call rfill(l+ 12,etaZ,WWW)
            call rfill(l+ 13,phiZ,WWW)
            call rfill(l+ 14,mZ,WWW)

            call rfill(l+ 15,eH,WWW)
            call rfill(l+ 16,ptH,WWW)
            call rfill(l+ 17,yH,WWW)
            call rfill(l+ 18,etaH,WWW)
            call rfill(l+ 19,phiH,WWW)
            call rfill(l+ 20,mH,WWW)

            call rfill(l+ 21,eZH,WWW)
            call rfill(l+ 22,ptZH,WWW)
            call rfill(l+ 23,yZH,WWW)
            call rfill(l+ 24,etaZH,WWW)
            call rfill(l+ 25,phiZH,WWW)
            call rfill(l+ 26,mZH,WWW)
            call rfill(l+ 27,eMu1,1.d0)
         enddo
      enddo
C
 999 return
      end

      subroutine wweights(wgts, sz, ibody, p4Mu1, p4Mu2, p4H)
c This program is a prototype for writing weights to a file
c The weights are printed to a file.
      double precision wgts(*), wsum, p4Mu1(0:3), p4Mu2(0:3), p4H(0:3)
      integer k, sz, ibody
      character(len=15) fn
      logical :: exist

      wsum = 0.d0

      do i = 2,sz
          wsum = wsum + wgts(i)
      end do
c print *, "The weight sum is", wsum

      if (wsum == 0d0) then
        return
      else

c print *, "Destination is 'weights_sam' "
c read *, fn
      fn = "weights_sam"

      inquire(file=fn, exist=exist)
      if (exist) then
       open(1, file=fn, status="old", position="append", action="write")
c print *, "~File named '",fn,"' was reopened"
      else
          open(1, file=fn, status="new", action="write")
          print *, "~File named '",fn,"' was created"
      end if

      write (1,6) ibody
C write (1,9) wgts(1), wgts(2), wgts(3), wgts(4)
      k = 0
      do i = 1, sz
         write (1,8) wgts(i)
      end do
      write (1,7) "p4Mu1: ", p4Mu1(0), p4Mu1(1), p4Mu1(2), p4Mu1(3)
      write (1,7) "p4Mu2: ", p4Mu2(0), p4Mu2(1), p4Mu2(2), p4Mu2(3)
      write (1,7) "p4H: ", p4H(0), p4H(1), p4H(2), p4H(3)
6 format (i5)
7 format(a8,d26.17," ",d26.17, " ",d26.17," ",d26.17)
8 format (d26.17)
9 format (d26.17, d26.17, d26.17, d26.17)
      close(1)

      end if
      end

      subroutine wdebug(p4Mu1, p4Mu2, mZ, p4H, mH, mZH)
c This program lets me write the values of xnorm to a file
      double precision p4Mu1(0:3), p4Mu2(0:3), p4H(0:3)
      double precision mZ, mH, mZH
      character(len=15) fn
      logical :: exist

      print *, "Destination is 'debug_sam' "
      fn = "debug_sam"

      inquire(file=fn, exist=exist)
      if (exist) then
       open(1, file=fn, status="old", position="append", action="write")
       print *, "~File named '",fn,"' was reopened"
      else
          open(1, file=fn, status="new", action="write")
          print *, "~File named '",fn,"' was created"
      end if

      write (1,*) "BEGIN EVENT:"
      write (1,8) "p4Mu1: ",p4Mu1(0),p4Mu1(1),p4Mu1(2),p4Mu1(3)
      write (1,8) "p4Mu2: ",p4Mu2(0),p4Mu2(1),p4Mu2(2),p4Mu2(3)
      write (1,8) "p4H: ",p4H(0),p4H(1),p4H(2),p4H(3)
      write (1,5) "Masses: ", mZ, mH, mZH
      write (1,*)
5 format (a8,d26.17,d26.17,d26.17)
8 format (a8,d26.17,d26.17,d26.17,d26.17)

      close(1)
      end

      subroutine wxnorm(xnorm)
c This program lets me write the values of xnorm to a file
      double precision xnorm
      character(len=15) fn
      logical :: exist

      print *, "Destination is 'xnorm_sam' "
      fn = "xnorm_sam"

      inquire(file=fn, exist=exist)
      if (exist) then
       open(1, file=fn, status="old", position="append", action="write")
       print *, "~File named '",fn,"' was reopened"
      else
          open(1, file=fn, status="new", action="write")
          print *, "~File named '",fn,"' was created"
      end if

      write (1,4) xnorm
4 format (d26.17)

      close(1)
      end

      function getrapidity(en,pl)
      implicit none
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
      parameter (tiny=1.d-8)
      xplus=en+pl
      xminus=en-pl
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
         if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny)then
            y=0.5d0*log( xplus/xminus )
         else
            y=sign(1.d0,pl)*1.d8
         endif
      else
         y=sign(1.d0,pl)*1.d8
      endif
      getrapidity=y
      return
      end

      function getpseudorap(en,ptx,pty,pl)
      implicit none
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
      parameter (tiny=1.d-5)
c
      pt=sqrt(ptx**2+pty**2)
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
        eta=sign(1.d0,pl)*1.d8
      else
        th=atan2(pt,pl)
        eta=-log(tan(th/2.d0))
      endif
      getpseudorap=eta
      return
      end