Modifying exisiting ej cut yields 0 cross section

Hi,

> For example, if I simply don't want to cut on, say final state particles 3 and 5 (this is just to give an example, it's not meant to be physical at this stage), and add the piece
> i .eq. VBFj1 .or. i .eq. VBFj2
> to the definition of 'notgood' in the standard emin cut below:

I guess that you miss some parenthesis to resolve the priority ordering between or and and operation.
Currently as soon as "i" is equal to "3" then notgood is "True" and therefore you remove the phase-space point (independently of the kinematics).

> P.S.: Olivier told me in a question I asked recently that I should modify the phase space integrator for my own cuts to keep efficiencies high. Is this maybe the case here again/as well (even though I modified an existing cut)? Unfortunately, I didn't really find a way on doing that, but would very much appreciate any help on that as well, if that could help preventing these kind of problems in the future!

Actually weaken an existing cut is even worse than adding one, since in this case we have learn from the "ej" value what is the minimum energy for each invariant mass in the system (including sqrts).
Now that we weaken this cut, all those relation are not valid anymore but the code will still assume that those minimum energy exists, creating potential bias in the generation.

For the optimisation, the main file to modify is setcuts.f where you define the minimum energy (sqrts) passing the cuts.
Then in myamp.f you have the typical energy for each resonances (with the definition of the xo variable) which control on how soft/strong the variable need to be generated.
This one is less critical and need to be modified only in presence of bad efficiency.

Cheers,

Olivier

> On 4 Oct 2017, at 13:04, Patrick Schaefers <email address hidden> wrote:
>
> New question #658919 on MadGraph5_aMC@NLO:
> https://answers.launchpad.net/mg5amcnlo/+question/658919
>
> Dear Developers,
>
> I am considering VBF processes and would like to cut on the energy of just the two VBF jets instead of cutting on all jets (i.e. in 'generate u u~ > u u~ g g', I'd like to apply the cut usually just on u and u~ and not on the gluons).
>
> While trying a couple of things, however, I noticed that even a very simple change in the default emin cut in cuts.f yields a zero cross section and I am not exactly sure why.
>
> For example, if I simply don't want to cut on, say final state particles 3 and 5 (this is just to give an example, it's not meant to be physical at this stage), and add the piece
> i .eq. VBFj1 .or. i .eq. VBFj2
> to the definition of 'notgood' in the standard emin cut below:
>
> integer VBFj1,VBFj2 ! added at the top of cuts.f
>
> VBFj1=3
> VBFj2=5
> c
> c E min & max cuts
> c
> do i=nincoming+1,nexternal
> if(debug) write (*,*) 'p(0,',i,')=',p(0,i),' ',emin(i),':',emax(i)
> notgood=(i .eq. VBFj1 .or. i .eq. VBFj2 .or. p(0,i) .le. emin(i)).or.
> & (emax(i).ge.0d0 .and. p(0,i) .gt. emax(i))
> if (notgood) then
> if(debug) write (*,*) i,' -> fails'
> passcuts=.false.
> return
> endif
> enddo
>
> Then I get the zero cross section error, even if ej is set to 0 in the run_card.
> On the other hand, however, if I just set the ej cut in the run_card and have it applied on all particles (i.e. with the default, unmodied cut), everything works fine, of course.
> I would have assumed that applying the cut on even fewer particles should actually work - do you know what I am missing here?
>
>
> Just for completeness, here is how I tried to physically correctly implement the cut (this worked for as long as I don't modify the ej cut as above - setting a value on the default cut works, however):
>
> C----------------------------------------------------------------------C
> C- MULTIBOSONS BACKGROUND INV_MASS CUT -C
> C----------------------------------------------------------------------C
>
> C maxDeta=0d0
> C VBFj1=0
> C VBFj2=0
> C MBEfailed=0
> C mInvMB=1000d0
> C firstPset=0
>
> C c-- RESET ptemp(0:3)
> C do j=0,3
> C ptemp(j)=0 ! for the first center particle
> C ptemp2(j)=0 ! for the second center particle
> C enddo
>
> C C-- FIND THE MAXIMUM DELTA ETA AND ITS CORRESPONDING JETS
>
> C do i=nincoming+1,nexternal
> C do j=i+1,nexternal
> C if( is_a_j(i) .and. is_a_j(j) .and.
> C & abs(rap(p(0,i))-rap(p(0,j))) .gt. maxDeta ) then
> C maxDeta=abs(rap(p(0,i))-rap(p(0,j)))
> C VBFj1=i
> C VBFj2=j
> C endif
> C enddo
> C enddo
>
> C-- NOW APPLY THE DELTA ETA CUT AND THE ADDITIONAL M_INV CUT
>
> C if (maxDeta .ge. deltaeta) then
> C do i=nincoming+1,nexternal
> C if (i .eq. VBFj1 .or. i .eq. VBFj2) cycle
> C do j=0,3
> C if (firstPset .eq. 0 ) then
> C ptemp(j) = ptemp(j) + p(j,i)
> C firstPset=1
> C elseif ( firstPset .eq. 1 ) then
> C ptemp2(j) = ptemp2(j) + p(j,i)
> C endif
> C enddo
> C enddo
> C if ( !ptemp(0) .lt. 150d0 .or. ptemp2(0) .lt. 150d0 .or.
> C & dsqrt( ((ptemp(0)+ptemp2(0)) * (ptemp(0)+ptemp2(0)))
> C & - ((ptemp(1)+ptemp2(1)) * (ptemp(1)+ptemp2(1)))
> C & - ((ptemp(2)+ptemp2(2)) * (ptemp(2)+ptemp2(2)))
> C & - ((ptemp(3)+ptemp2(3)) * (ptemp(3)+ptemp2(3))) )
> C & .lt. mInvMB ) then
> C passcuts=.false.
> C return
> C endif
> C elseif (maxDeta .lt. deltaeta) then
> C passcuts=.false.
> C return
> C endif
>
> So I am first finding the two VBF jets via the maximum Delta eta and label them VBFj1 and VBFj2. Then I set a cut on the invariant mass on everything but the two VBF jets.
> Again, this cut worked just fine for me, even if I set a value on ej in the run_card, but only if I do not touch the emin cut in cuts.f (i.e. if I do not add the above piece ' i .eq. VBFj1 .or. i .eq. VBFj2' in the definition of 'notgood'.)
>
>
> Thank you very much for your help in advance!
>
> Cheers,
> Patrick
>
>
> P.S.: Olivier told me in a question I asked recently that I should modify the phase space integrator for my own cuts to keep efficiencies high. Is this maybe the case here again/as well (even though I modified an existing cut)? Unfortunately, I didn't really find a way on doing that, but would very much appreciate any help on that as well, if that could help preventing these kind of problems in the future!
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi Olivier,

thank you for your fast reply!

> I guess that you miss some parenthesis to resolve the priority ordering between or and and operation.
> Currently as soon as "i" is equal to "3" then notgood is "True" and therefore you remove the phase-space point (independently of the kinematics).

That is exactly what I would like to do. If, as in this example, I would check particle 3, I do not want to apply the cut, regardless of whether it would have passed the energy condition or not. The same should be true once I check particle 5.

> Actually weaken an existing cut is even worse than adding one, since in this case we have learn from the "ej" value what is the minimum energy for each invariant mass in the system (including sqrts).

Okay, I think I understood the underlying mechanics now.
Would this in my case mean that it would be better to just write a new cut for ej, where I take into account that it should not trigger for all jets, but only those passing my VBF criteria? The VBF criteria would probably need to be defined in setcuts.f then, right? In order to avoid the problem you mentioned above. But would that be possible?

The criterion for when I want to cut on ej is:
Find the two jets with the largest Delta eta and then apply/check for the ej cut.
I wrote the code to find the two jets as (just as above)

  do i=nincoming+1,nexternal
    do j=i+1,nexternal
      if( is_a_j(i) .and. is_a_j(j) .and.
        & abs(rap(p(0,i))-rap(p(0,j))) .gt. maxDeta ) then
        maxDeta=abs(rap(p(0,i))-rap(p(0,j)))
        VBFj1=i
        VBFj2=j
      endif
    enddo
  enddo

Could I just plug this into setcuts.f and transfer the VBFj1 and VBFj2 I found into cuts.f into the new ej cut?

Sorry for this wall of text again, I just want to try to be as precise as possible.

Cheers,
Patrick

Hi,

> That is exactly what I would like to do. If, as in this example, I would
> check particle 3, I do not want to apply the cut, regardless of whether
> it would have passed the energy condition or not. The same should be
> true once I check particle 5.

Yes but what you wrote is

> do i=1,5
> if (i.eq.3) then
> pass_cuts = .false.
> return
> endif
> enddo

This function always return .false. whatever the kinematics is (and this is not what you want to do)

you should rather do something like
> do i=1,5
> if(i.eq.3.or.i.eq.5)then
> if (fail some kinematic cut)then
> pass_cuts = .false.
> endif
> endif
> enddo

You can obviously do it in a single line as you tried to do it, but you have to do it carefully (and this single line syntax is not very human readable --the reason why you failed to write it correctly--).

> do i=nincoming+1,nexternal
> do j=i+1,nexternal
> if( is_a_j(i) .and. is_a_j(j) .and.
> & abs(rap(p(0,i))-rap(p(0,j))) .gt. maxDeta ) then
> maxDeta=abs(rap(p(0,i))-rap(p(0,j)))
> VBFj1=i
> VBFj2=j
> endif
> enddo
> enddo
>
> Could I just plug this into setcuts.f and transfer the VBFj1 and VBFj2 I
> found into cuts.f into the new ej cut?

I should say that you should do that directly in cuts.f and not in setcuts.f
Note that we have some default cut already existing for VBF.
I would suggest to look at those.

> On 4 Oct 2017, at 17:22, Patrick Schaefers <email address hidden> wrote:
>
> Question #658919 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/658919
>
> Status: Answered => Open
>
> Patrick Schaefers is still having a problem:
> Hi Olivier,
>
> thank you for your fast reply!
>
>> I guess that you miss some parenthesis to resolve the priority ordering between or and and operation.
>> Currently as soon as "i" is equal to "3" then notgood is "True" and therefore you remove the phase-space point (independently of the kinematics).
>
> That is exactly what I would like to do. If, as in this example, I would
> check particle 3, I do not want to apply the cut, regardless of whether
> it would have passed the energy condition or not. The same should be
> true once I check particle 5.
>
>> Actually weaken an existing cut is even worse than adding one, since
> in this case we have learn from the "ej" value what is the minimum
> energy for each invariant mass in the system (including sqrts).
>
> Okay, I think I understood the underlying mechanics now.
> Would this in my case mean that it would be better to just write a new cut for ej, where I take into account that it should not trigger for all jets, but only those passing my VBF criteria? The VBF criteria would probably need to be defined in setcuts.f then, right? In order to avoid the problem you mentioned above. But would that be possible?
>
> The criterion for when I want to cut on ej is:
> Find the two jets with the largest Delta eta and then apply/check for the ej cut.
> I wrote the code to find the two jets as (just as above)
>
> do i=nincoming+1,nexternal
> do j=i+1,nexternal
> if( is_a_j(i) .and. is_a_j(j) .and.
> & abs(rap(p(0,i))-rap(p(0,j))) .gt. maxDeta ) then
> maxDeta=abs(rap(p(0,i))-rap(p(0,j)))
> VBFj1=i
> VBFj2=j
> endif
> enddo
> enddo
>
> Could I just plug this into setcuts.f and transfer the VBFj1 and VBFj2 I
> found into cuts.f into the new ej cut?
>
> Sorry for this wall of text again, I just want to try to be as precise
> as possible.
>
> Cheers,
> Patrick
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.