FOA: Rather weird error message

Dear Afiq,

sorry if I'm a bit confused, but it seems you removed all weight informations from your
routines, you removed the arguments from analysis_begin. I guess the problem may be
due to this.
Could you confirm that if you use one of the template analyses, like analysis_root_pp_V.f,
it works technically in your setup (i.e. makes the code compile and produces some plots
at the end of the run)?
If this is the case, could you check what happens if you modify this working analysis very
minimally, just changing things relevant to physical observables, without touching the
structure of the routines (arguments and so on)?

Thanks.
Cheers.
Paolo

Indeed I did, because as I mentioned, I wanted to find a minimal case that would work. Anyway, here is one more example, based on the pp_V one, where I changed only the nexternal to match the hmumu process and the pdg ID to 25. I expected by doing this I'll end up with histos of Higgs kinematics, but nope, all I got was the error.

[1]
c
c Example analysis for "p p > w+ [QCD]" process.
c Example analysis for "p p > w- [QCD]" process.
c Example analysis for "p p > z [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 bin,xmi,xms,pi
      parameter (pi=3.14159265358979312d0)
      call open_root_file()
      nwgt_analysis=nwgt
      xmi=100.d0
      xms=150.d0
      bin=1.0d0
      do j=1,2
      do kk=1,nwgt_analysis
      l=(kk-1)*10+(j-1)*5
      call rbook(l+ 1,'V pt '//cc(j)//weights_info(kk)
     & ,2.d0,0.d0,200.d0)
      call rbook(l+ 2,'V log pt '//cc(j)//weights_info(kk)
     & ,0.05d0,0.d0,5.d0)
      call rbook(l+ 3,'V y '//cc(j)//weights_info(kk)
     & ,0.25d0,-9.d0,9.d0)
      call rbook(l+ 4,'V eta '//cc(j)//weights_info(kk)
     & ,0.25d0,-9.d0,9.d0)
      call rbook(l+ 5,'mV '//cc(j)//weights_info(kk)
     & ,bin,xmi,xms)
      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
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)*10+(i-1)*5
      do jj=1,5
        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,pv(0:3),xmv,ptv,yv,etav
      double precision getrapidity,getpseudorap,dot
      external getrapidity,getpseudorap,dot
      if (nexternal.ne.6) then
         write (*,*) 'error #1 in analysis_fill: '/
     & /'only for process "p p > V [QCD]"'
         stop 1
      endif
      if (.not. (abs(ipdg(1)).le.5 .or. ipdg(1).eq.21)) then
         write (*,*) 'error #2 in analysis_fill: '/
     & /'only for process "p p > V [QCD]"'
         stop 1
      endif
      if (.not. (abs(ipdg(2)).le.5 .or. ipdg(2).eq.21)) then
         write (*,*) 'error #3 in analysis_fill: '/
     & /'only for process "p p > V [QCD]"'
         stop 1
      endif
      if (.not. (abs(ipdg(4)).le.5 .or. ipdg(4).eq.21)) then
         write (*,*) 'error #4 in analysis_fill: '/
     & /'only for process "p p > V [QCD]"'
         stop 1
      endif
      if (abs(ipdg(3)).ne.25) then
         write (*,*) 'error #5 in analysis_fill: '/
     & /'only for process "p p > V [QCD]"'
         stop 1
      endif
C
      do i=0,3
        pv(i)=p(i,3)
      enddo
      xmv=sqrt(max(dot(pv,pv),0d0))
      ptv=sqrt(max(pv(1)**2+pv(2)**2,0d0))
      yv=getrapidity(pv(0),pv(3))
      etav=getpseudorap(pv(0),pv(1),pv(2),pv(3))
C
      do i=1,2
         do kk=1,nwgt_analysis
            www=wgts(kk)
            l=(kk-1)*10+(i-1)*5
            if (ibody.ne.3 .and.i.eq.2) cycle
            call rfill(l+1,ptv,WWW)
            if(ptv.gt.0) call rfill(l+2,log10(ptv),WWW)
            call rfill(l+3,yv,WWW)
            call rfill(l+4,etav,WWW)
            call rfill(l+5,xmv,WWW)
         enddo
      enddo
C
 999 return
      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

Hi Afiq,

by comparing your analysis and the default pp_V one, I can only
wonder if the ordering of particles in MG5_aMC corresponds to
that assumed by the checks at the beginning of analysis_fill.

I mean, you assume that the Higgs is particle number 3, and so
on: maybe what you are experiencing is just one of these checks
failing.

I would suggest to do a couple of things: comment all these checks
out (namely the first command in that routine would become 'do i=0,3’)
and see if it goes through; check particle ordering, for example by
inspecting one of the files SubProcesses/P0_*/matrix_*.ps.

Just to be sure: could you confirm that by generating p p > W+
and running analysis_root_pp_V.f on it, you get some plots out?
If this is not the case, the whole problem could be caused by something
else, and not the specific analysis used.

Let me know.
Cheers.
Paolo

On 01 Oct 2015, at 09:27, Afiq Aize <email address hidden> wrote:

> Question #271942 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/271942
>
> Status: Answered => Open
>
> Afiq Aize is still having a problem:
> Indeed I did, because as I mentioned, I wanted to find a minimal case
> that would work. Anyway, here is one more example, based on the pp_V
> one, where I changed only the nexternal to match the hmumu process and
> the pdg ID to 25. I expected by doing this I'll end up with histos of
> Higgs kinematics, but nope, all I got was the error.
>
> [1]
> c
> c Example analysis for "p p > w+ [QCD]" process.
> c Example analysis for "p p > w- [QCD]" process.
> c Example analysis for "p p > z [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 bin,xmi,xms,pi
> parameter (pi=3.14159265358979312d0)
> call open_root_file()
> nwgt_analysis=nwgt
> xmi=100.d0
> xms=150.d0
> bin=1.0d0
> do j=1,2
> do kk=1,nwgt_analysis
> l=(kk-1)*10+(j-1)*5
> call rbook(l+ 1,'V pt '//cc(j)//weights_info(kk)
> & ,2.d0,0.d0,200.d0)
> call rbook(l+ 2,'V log pt '//cc(j)//weights_info(kk)
> & ,0.05d0,0.d0,5.d0)
> call rbook(l+ 3,'V y '//cc(j)//weights_info(kk)
> & ,0.25d0,-9.d0,9.d0)
> call rbook(l+ 4,'V eta '//cc(j)//weights_info(kk)
> & ,0.25d0,-9.d0,9.d0)
> call rbook(l+ 5,'mV '//cc(j)//weights_info(kk)
> & ,bin,xmi,xms)
> 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
> 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)*10+(i-1)*5
> do jj=1,5
> 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,pv(0:3),xmv,ptv,yv,etav
> double precision getrapidity,getpseudorap,dot
> external getrapidity,getpseudorap,dot
> if (nexternal.ne.6) then
> write (*,*) 'error #1 in analysis_fill: '/
> & /'only for process "p p > V [QCD]"'
> stop 1
> endif
> if (.not. (abs(ipdg(1)).le.5 .or. ipdg(1).eq.21)) then
> write (*,*) 'error #2 in analysis_fill: '/
> & /'only for process "p p > V [QCD]"'
> stop 1
> endif
> if (.not. (abs(ipdg(2)).le.5 .or. ipdg(2).eq.21)) then
> write (*,*) 'error #3 in analysis_fill: '/
> & /'only for process "p p > V [QCD]"'
> stop 1
> endif
> if (.not. (abs(ipdg(4)).le.5 .or. ipdg(4).eq.21)) then
> write (*,*) 'error #4 in analysis_fill: '/
> & /'only for process "p p > V [QCD]"'
> stop 1
> endif
> if (abs(ipdg(3)).ne.25) then
> write (*,*) 'error #5 in analysis_fill: '/
> & /'only for process "p p > V [QCD]"'
> stop 1
> endif
> C
> do i=0,3
> pv(i)=p(i,3)
> enddo
> xmv=sqrt(max(dot(pv,pv),0d0))
> ptv=sqrt(max(pv(1)**2+pv(2)**2,0d0))
> yv=getrapidity(pv(0),pv(3))
> etav=getpseudorap(pv(0),pv(1),pv(2),pv(3))
> C
> do i=1,2
> do kk=1,nwgt_analysis
> www=wgts(kk)
> l=(kk-1)*10+(i-1)*5
> if (ibody.ne.3 .and.i.eq.2) cycle
> call rfill(l+1,ptv,WWW)
> if(ptv.gt.0) call rfill(l+2,log10(ptv),WWW)
> call rfill(l+3,yv,WWW)
> call rfill(l+4,etav,WWW)
> call rfill(l+5,xmv,WWW)
> enddo
> enddo
> C
> 999 return
> 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
>
> --
> You received this question notification because you are subscribed to
> the question.

1- Commenting all the checks: Ok! So now we can assume the checks cause the error.

2- Indeed the particle ordering is correct. Both the matrix_*.ps and the LHE file indicated Higgs is particle 3. But the error... why?

As for *_pp_V.f, yes I verified that it works. However something a bit strange, I'm getting way too many plots, I suspect this is because of these lines declaring too many histograms:

      do j=1,2
      do kk=1,nwgt_analysis
      l=(kk-1)*10+(j-1)*5

Can you please educate me what does it mean, especially that 10 and 5?

Best regards,
Afiq

Ok, disabling the checks one by one, the code works fine now.

However, I still don't understand about the too-many histograms thing. Can you please teach me about them, and how to get the simplest case of one-histogram-per-variable that correctly accounts for the weights?

Best regards,
Afiq

Hi Paolo,

Many thanks for the guidance. I still need to understand deeper the scales, but am I right that basically the histograms simply correspond to 'physics cases' (which are denoted by the weights_info() label) , say, for a lack of better phrase, i.e. for the _pp_V.f example I am getting ~100 histograms because there are ~20 different scales, so I am getting ~20 * 5 (the number of variables I asked) histograms?

Similarly, I saw in different templates different multiplication factors, like 52, 26 etc, they all serve the same purpose of ensuring the histogram number does not overlap, and the precise value is determined by understanding how many cases are there for particular processes and coming up with suitable "safety numbers"?

If I may suggest something, it would be very helpful if MG5 notifies us in cases like this i.e. the failures come from user code (something like "you might want to check SubProcesses/P0*/GF*/log.txt to understand the issue better").

Thanks again and sorry for all the troubles caused.

Best regards,
Afiq

Thanks Paolo Torrielli, that solved my question.

Hi Afiq,

> Many thanks for the guidance. I still need to understand deeper the
> scales, but am I right that basically the histograms simply correspond
> to 'physics cases' (which are denoted by the weights_info() label) ,
> say, for a lack of better phrase, i.e. for the _pp_V.f example I am
> getting ~100 histograms because there are ~20 different scales, so I am
> getting ~20 * 5 (the number of variables I asked) histograms?
The 5 is indeed the number of variables, while the 20 is a 2*10.

2 is: one for NLO, one for LO (indeed note that the array cc(j) is
data cc/'|T@NLO','|T@LO ‘/, so histograms for which j=1 are NLO,
with j=2 are LO).
Notice also another thing: in analysis_fill the condition on whether
it is Born or not is “if (ibody.ne.3 .and.i.eq.2) cycle", which you
have to adapt to your multiplicity (it will become ibody.ne.5).

The 10 comes from scales. Given an observable, out of the ten histograms
you get for it, the first corresponds to the central value for
renormalisation and factorisation scales, i.e. the one curve you would
get setting to False those reweighting parameters in the run_card.
The other nine correspond to different combinations of renormalisation
and factorisation scales: there are three values considered for
renormalisation scale (namely muR = (1, 0.5, 2)*muR_central) and three
for the factorisation scale (namely muF = (1, 0.5, 2)*muF_central),
hence nine combinations (note that muR=1*muR_central and muF=1*muF_central
correspond to the first curve).
The theoretical uncertainty band will be the envelope of these variations.

I would also suggest to use the internal plotting package (HwU) we
provide, where the combination of all these variation histograms is
automatic, and you just see one histogram per observable (with the
associated theoretical-uncertainty band).
The HwU package is activated by using analysis_HwU_pp_*.f instead of
the root analysis, and the management of plots (like combinations,
rescaling, rebinning, ..) can be performed with a flexible python
script that you can run on the plots you generated (the script is
in madgraph/various/histograms.py and you can run it with --help to
see all of the options it has.) The output is in gnuplot format.

> Similarly, I saw in different templates different multiplication
> factors, like 52, 26 etc, they all serve the same purpose of ensuring
> the histogram number does not overlap, and the precise value is
> determined by understanding how many cases are there for particular
> processes and coming up with suitable "safety numbers”?
Yes, precisely.

> If I may suggest something, it would be very helpful if MG5 notifies us
> in cases like this i.e. the failures come from user code (something like
> "you might want to check SubProcesses/P0*/GF*/log.txt to understand the
> issue better").
Thanks for the suggestion!

Cheers.
Paolo

Hi Paolo,

That is a very thorough explanation. Thanks a lot for your time.

Best regards,
Afiq

Hi Afiq,

when replying to your latest email I wrote something wrong:

> Notice also another thing: in analysis_fill the condition on whether
> it is Born or not is “if (ibody.ne.3 .and.i.eq.2) cycle", which you
> have to adapt to your multiplicity (it will become ibody.ne.5).

This is not true: ibody = 3 always means Born, it does *not*
correspond to the multiplicity you have at Born level.

So, the relevant line in the analysis routine is always

         if (ibody.ne.3 .and.i.eq.2) cycle

Sorry for the hassle.
Cheers.
Paolo

Hi Paolo,

Uh oh, okay. But I changed it to 5 and the code ran fine, can you please explain what was I doing in fact with setting?

If you don't mind, I would be delighted to learn about the possible value range and their meaning.

Best regards,
Afiq

Hi Afiq,

with

if (ibody.ne.5 .and.i.eq.2) cycle

it ran fine simply because condition ‘ibody.ne.5’ is always true, hence
the command ‘cycle’ was always activated when i=2.
Namely, you should have one NLO curve (corresponding to i=1) and no LO
one (i.e. you always skip i=2).

By setting

if (ibody.ne.3 .and.i.eq.2) cycle

the i=2 step is skipped (due to cycle) only when ibody!=3: namely you fill
the i=2 histograms only when ibody=3, i.e. in the Born configurations.
This way you should get two curves, one NLO, the other LO.

As for ibody values, the only relevant for the user should be this ibody =
3, for plotting Born quantities. ibody=1 means (n+1)-body contribution to
NLO quantities, ibody=2 means the n-body contribution to NLO quantities
that does not come from Born, while ibody=3 is Born. To get NLO one keeps
the three values, while for Born only ibody=3.

Cheers.
Paolo

On 08 Oct 2015, at 09:07, Afiq Aize <email address hidden> wrote:

> Question #271942 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/271942
>
> Afiq Aize posted a new comment:
> Hi Paolo,
>
> Uh oh, okay. But I changed it to 5 and the code ran fine, can you please
> explain what was I doing in fact with setting?
>
> If you don't mind, I would be delighted to learn about the possible
> value range and their meaning.
>
> Best regards,
> Afiq
>
> --
> You received this question notification because you are subscribed to
> the question.

Hi Paolo,

Indeed as you say with ibody = 5 I don't get the LO histograms. If I understand your answer correctly, then changing to ibody = 3 affects also NLO distributions because with 5 the Born contribution is never saved -> the NLO contribution is in fact (NLO - LO) = wrong.

If this is correct, then I guess it explains the change in scale for NLO distributions when I tried setting it to 3. Otherwise, I would like to ask the explanation for the change in scale I observed. (going from 5 to 3 increased it around 100 times larger).

Best regards,
Afiq

Hi Afiq,

> Indeed as you say with ibody = 5 I don't get the LO histograms. If I
> understand your answer correctly, then changing to ibody = 3 affects
> also NLO distributions because with 5 the Born contribution is never
> saved -> the NLO contribution is in fact (NLO - LO) = wrong.
No, the NLO contribution should not change.
By NLO we mean LO + real + virtual, not only real + virtual.
If i.eq.1 the command ‘cycle’ is never executed, hence you are filling
your histograms with all configurations you have (Born configurations,
virtual configurations, real configurations).

> If this is correct, then I guess it explains the change in scale for NLO
> distributions when I tried setting it to 3. Otherwise, I would like to
> ask the explanation for the change in scale I observed. (going from 5 to
> 3 increased it around 100 times larger).
This is unexpected. Could you check if now the normalisation of the inclusive
histograms you have corresponds to the total NLO cross section? Do you confirm
that now you have two sets of histograms, one for NLO one for LO?

Cheers.
Paolo

Yes and yes. Surprisingly, considering the y-scale is different.

Wait, no. The total integral is different from ibody 5 to ibody 3, with ibody 5 being somewhat less. Ibody 3 does seem to correspond to xsec but ibody 5... don't know the number.

Hi Afiq,
ibody = 5 is wrong, so we simply discard it.
If ibody = 3 is correctly normalised (i.e. = the total
cross section), this is the end of the story.
Cheers.
Paolo

On 08 Oct 2015, at 10:33, Afiq Aize <email address hidden> wrote:

> Question #271942 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/271942
>
> Afiq Aize posted a new comment:
> Wait, no. The total integral is different from ibody 5 to ibody 3, with
> ibody 5 being somewhat less. Ibody 3 does seem to correspond to xsec but
> ibody 5... don't know the number.
>
> --
> You received this question notification because you are subscribed to
> the question.

Alright, thanks!