scarrazza / apfel

A PDF Evolution Library
http://apfel.mi.infn.it/
GNU General Public License v3.0
12 stars 9 forks source link

DIS: CC FONLL-A scale variations #27

Open felixhekhorn opened 4 years ago

felixhekhorn commented 4 years ago
non-scale variated
           x           Q2     APFEL    yadism  yadism_error  rel_err[%]
0   0.000317    90.000000  0.001365  0.001365  1.052304e-11   -0.003618
1   0.001007    90.000000  0.003636  0.003636  2.957269e-11   -0.011747
2   0.003195    90.000000  0.009640  0.009638  1.129206e-10   -0.014704
3   0.010138    90.000000  0.025498  0.025494  3.063871e-10   -0.015907
4   0.032169    90.000000  0.067195  0.067179  8.827141e-10   -0.024268
5   0.102077    90.000000  0.173154  0.173063  2.515006e-09   -0.052814
6   0.304545    90.000000  0.376150  0.376114  4.893392e-09   -0.009677
7   0.536364    90.000000  0.460563  0.457903  7.999385e-09   -0.577581
8   0.768182    90.000000  0.354356  0.363529  5.272257e-09    2.588783
9   1.000000    90.000000  0.000000  0.000000  0.000000e+00         NaN
10  0.001000     4.000000  0.007303  0.007293  1.185785e-11   -0.130942
11  0.001000     5.348935  0.006370  0.006368  2.650673e-11   -0.030171
12  0.001000     7.152778  0.005681  0.005677  4.638424e-11   -0.057296
13  0.001000     9.564937  0.005165  0.005157  1.354225e-11   -0.145549
14  0.001000    12.790557  0.004764  0.004761  1.129663e-11   -0.057178
15  0.001000    17.103966  0.004459  0.004458  9.166199e-12   -0.012656
16  0.001000    22.872003  0.004228  0.004228  1.874557e-11    0.005051
17  0.001000    30.585218  0.004046  0.004046  1.615744e-11   -0.010532
18  0.001000    40.899589  0.003899  0.003899  1.228940e-11   -0.010231
19  0.001000    54.692316  0.003778  0.003778  1.823467e-11   -0.007239
20  0.001000    73.136417  0.003677  0.003677  1.288867e-11   -0.001218
21  0.001000    97.800494  0.003592  0.003592  1.575276e-11   -0.011739
22  0.001000   130.782133  0.003519  0.003518  1.637083e-11   -0.014816
23  0.001000   174.886299  0.003454  0.003454  1.217347e-11   -0.015418
24  0.001000   233.863883  0.003397  0.003396  1.171409e-11   -0.014889
25  0.001000   312.730705  0.003345  0.003345  1.807173e-11   -0.016899
26  0.001000   418.194092  0.003299  0.003298  1.244429e-11   -0.018640
27  0.001000   559.223306  0.003256  0.003255  8.366275e-12   -0.018543
28  0.001000   747.812346  0.003216  0.003216  1.052677e-11   -0.018331
29  0.001000  1000.000000  0.003179  0.003179  1.092562e-11   -0.018418
scale variated
           x           Q2     APFEL    yadism  yadism_error  rel_err[%]
0   0.000317    90.000000  0.001856  0.001615  6.113991e-12  -12.974212
1   0.001007    90.000000  0.004953  0.004306  2.015077e-11  -13.052670
2   0.003195    90.000000  0.013050  0.011373  5.598218e-11  -12.849692
3   0.010138    90.000000  0.033837  0.029727  1.880561e-10  -12.144662
4   0.032169    90.000000  0.085403  0.076386  5.067254e-10  -10.558479
5   0.102077    90.000000  0.202139  0.187519  1.417134e-09   -7.232290
6   0.304545    90.000000  0.372621  0.373441  2.981850e-09    0.220166
7   0.536364    90.000000  0.386012  0.418997  4.154874e-09    8.545153
8   0.768182    90.000000  0.239398  0.305492  2.403995e-09   27.608191
9   1.000000    90.000000  0.000000  0.000000  0.000000e+00         NaN
10  0.001000     4.000000  0.009186  0.009174  5.671767e-12   -0.130288
11  0.001000     5.348935  0.007951  0.007948  2.224786e-11   -0.039398
12  0.001000     7.152778  0.007037  0.007033  2.030852e-11   -0.055678
13  0.001000     9.564937  0.006348  0.006343  6.949554e-12   -0.089371
14  0.001000    12.790557  0.005820  0.005816  4.434129e-12   -0.070409
15  0.001000    17.103966  0.006214  0.005413  1.974105e-11  -12.893534
16  0.001000    22.872003  0.005873  0.005106  2.391705e-11  -13.057690
17  0.001000    30.585218  0.005600  0.004862  2.468689e-11  -13.166203
18  0.001000    40.899589  0.005374  0.004665  2.013577e-11  -13.195396
19  0.001000    54.692316  0.005185  0.004502  2.422744e-11  -13.169462
20  0.001000    73.136417  0.005024  0.004366  1.795298e-11  -13.105272
21  0.001000    97.800494  0.004886  0.004250  2.000938e-11  -13.025702
22  0.001000   130.782133  0.004765  0.004149  1.841735e-11  -12.922172
23  0.001000   174.886299  0.004657  0.004061  1.918520e-11  -12.802721
24  0.001000   233.863883  0.004560  0.003982  1.439796e-11  -12.673503
25  0.001000   312.730705  0.004472  0.003911  1.593528e-11  -12.541315
26  0.001000   418.194092  0.004392  0.003847  1.555100e-11  -12.405532
27  0.001000   559.223306  0.004317  0.003788  1.470485e-11  -12.267555
28  0.001000   747.812346  0.004248  0.003733  1.484172e-11  -12.129212
29  0.001000  1000.000000  0.004183  0.003682  1.069974e-11  -11.992034
diff
           x           Q2     APFEL    yadism  yadism_error  rel_err[%]
0   0.000317    90.000000 -0.000491 -0.000250  1.663703e-11  -49.071786
1   0.001007    90.000000 -0.001316 -0.000670  4.972346e-11  -49.083270
2   0.003195    90.000000 -0.003410 -0.001735  1.689027e-10  -49.129734
3   0.010138    90.000000 -0.008339 -0.004234  4.944432e-10  -49.229239
4   0.032169    90.000000 -0.018208 -0.009207  1.389440e-09  -49.434442
5   0.102077    90.000000 -0.028984 -0.014457  3.932141e-09  -50.122870
6   0.304545    90.000000  0.003530  0.002673  7.875242e-09  -24.274422
7   0.536364    90.000000  0.074551  0.038906  1.215426e-08  -47.813207
8   0.768182    90.000000  0.114958  0.058038  7.676251e-09  -49.513885
9   1.000000    90.000000  0.000000  0.000000  0.000000e+00    0.000000
10  0.001000     4.000000 -0.001883 -0.001881  1.752961e-11   -0.127752
11  0.001000     5.348935 -0.001581 -0.001580  4.875459e-11   -0.076566
12  0.001000     7.152778 -0.001356 -0.001355  6.669276e-11   -0.048901
13  0.001000     9.564937 -0.001183 -0.001185  2.049181e-11    0.155818
14  0.001000    12.790557 -0.001056 -0.001055  1.573076e-11   -0.130089
15  0.001000    17.103966 -0.001755 -0.000955  2.890725e-11  -45.617532
16  0.001000    22.872003 -0.001645 -0.000878  4.266262e-11  -46.627284
17  0.001000    30.585218 -0.001553 -0.000817  4.084433e-11  -47.432279
18  0.001000    40.899589 -0.001475 -0.000766  3.242517e-11  -48.050363
19  0.001000    54.692316 -0.001407 -0.000724  4.246211e-11  -48.519009
20  0.001000    73.136417 -0.001347 -0.000689  3.084166e-11  -48.878320
21  0.001000    97.800494 -0.001294 -0.000658  3.576214e-11  -49.153844
22  0.001000   130.782133 -0.001246 -0.000631  3.478819e-11  -49.362926
23  0.001000   174.886299 -0.001203 -0.000607  3.135867e-11  -49.518829
24  0.001000   233.863883 -0.001163 -0.000586  2.611204e-11  -49.635516
25  0.001000   312.730705 -0.001127 -0.000567  3.400701e-11  -49.724830
26  0.001000   418.194092 -0.001093 -0.000549  2.799528e-11  -49.791142
27  0.001000   559.223306 -0.001061 -0.000532  2.307112e-11  -49.841461
28  0.001000   747.812346 -0.001032 -0.000517  2.536849e-11  -49.878594
29  0.001000  1000.000000 -0.001004 -0.000503  2.162536e-11  -49.906513
APFEL banner
 WARNING: FONLL-A is a VFN scheme
          ... setting VFNS PDF evolution
 WARNING: FONLL-A is NLO scheme
          ... setting NLO perturbative order
 WARNING: if there are external grids they cannot be locked
          ... unlocking subgrids

 Welcome to 
      _/_/_/    _/_/_/_/   _/_/_/_/   _/_/_/_/   _/
    _/    _/   _/    _/   _/         _/         _/
   _/_/_/_/   _/_/_/_/   _/_/_/     _/_/_/     _/
  _/    _/   _/         _/         _/         _/
 _/    _/   _/         _/         _/_/_/_/   _/_/_/_/
 _____v3.0.4 A PDF Evolution Library, arXiv:1310.1394
      Authors: V. Bertone, S. Carrazza, J. Rojo

 Report of the evolution parameters:

 QCD evolution
 Space-like evolution (PDFs)
 Unpolarized evolution
 Evolution scheme: VFNS at N1LO
 Solution of the DGLAP equation: 'exactalpha' with maximum 6 active flavours
 Solution of the coupling equations: 'exact' with maximum 6 active flavours
 Coupling reference value:
 - AlphaQCD( 91.2000 GeV) =  0.118000
 Pole heavy quark masses:
 - Mc =   2.0000 GeV
 - Mb =   4.0000 GeV
 - Mt = 173.0700 GeV
 The matching thresholds coincide with the physical masses
 muR / muF =  1.0000

 Allowed evolution range [   1.0000 :  10000.0000 ] GeV
 Fast evolution enabled

 Initialization of the evolution completed in   5.254 s

 Report of the electroweak parameters:

 Mass of the Z = 91.188 GeV
 Mass of the W = 80.398 GeV
 Mass of the proton = 0.9380 GeV
 sin^2(thetaW) = 0.2313
 GFermi = 1.16638E-05
       | 0.9743 0.2253 0.0035 |
 CKM = | 0.2252 0.9735 0.0410 |
       | 0.0086 0.0403 0.9992 |
 Z propagator correction = 0.00000

 Report of the DIS parameters:

 Computation in the FONLL-A mass scheme
 Charged Current (CC) process
 Scattering electron - proton   
 muR / Q =  1.0000
 muF / Q =  0.5000
 (No scale variation in the evolution)
 Target Mass corrections disabled
 FONLL damping factor disabled for all heavy quarks
 Intrinsic charm disabled
alecandido commented 4 years ago

We forgot to specify that the issue is related only to the quarks' channels, but it is not present for the gluon channel.

vbertone commented 4 years ago

@felixhekhorn @AleCandido could you please be more specific as to where the missing factor 2 should be? Also, from the banner I see that when the evolution is initialised you have:

muR / muF =  1.0000

but then when the DIS module is initialised you have:

muR / Q =  1.0000
muF / Q =  0.5000

which are clearly incompatible. This suggests that you are somehow misusing the code because that ratio is automatically adjusted when the code is used correctly (incidentally, this may explain the discontinuity at the bottom mass). As a matter of fact, this very simple code:

#include "APFEL/APFEL.h"
int main()
{
  APFEL::SetMassScheme("FONLL-A");
  APFEL::SetProcessDIS("CC");
  APFEL::SetFacQRatio(0.5);
  APFEL::InitializeAPFEL_DIS();
  return 0;
}

produces the following banner:

 QCD evolution
 Space-like evolution (PDFs)
 Unpolarized evolution
 Evolution scheme: VFNS at N1LO
 Solution of the DGLAP equation: 'exactalpha' with maximum 6 active flavours
 Solution of the coupling equations: 'exact' with maximum 6 active flavours
 Coupling reference value:
 - AlphaQCD(  1.4142 GeV) =  0.350000
 Pole heavy quark masses:
 - Mc =   1.4142 GeV
 - Mb =   4.5000 GeV
 - Mt = 175.0000 GeV
 The matching thresholds coincide with the physical masses
 muR / muF =  2.0000

 Allowed evolution range [   1.0000 :  10000.0000 ] GeV
 The internal subgrids will be locked
 Fast evolution enabled

 Initialization of the evolution completed in   0.870 s

 Report of the electroweak parameters:

 Mass of the Z = 91.188 GeV
 Mass of the W = 80.385 GeV
 Mass of the proton = 0.9383 GeV
 sin^2(thetaW) = 0.2313
 GFermi = 1.16638E-05
       | 0.9743 0.2254 0.0036 |
 CKM = | 0.2252 0.9734 0.0414 |
       | 0.0089 0.0405 0.9991 |
 Z propagator correction = 0.00000

 Report of the DIS parameters:

 Computation in the FONLL-A mass scheme
 Charged Current (CC) process
 Scattering electron - proton   
 muR / Q =  1.0000
 muF / Q =  0.5000
 Target Mass corrections disabled
 FONLL damping factor for charm enabled with suppression power = 2
 FONLL damping factor for bottom enabled with suppression power = 2
 FONLL damping factor for top enabled with suppression power = 2
 Intrinsic charm disabled

where you can see that the ratio muR / muF is adjusted to 2 as appropriate. Can you share the code you are using to produce these predictions?

alecandido commented 4 years ago

Factor of two...

The missing factor of two is just the one detected from the numbers (see the diff spoiler in @felixhekhorn OP), there is no place in APFEL's code in which we could put it to fix the code itself, because instead we believe the scale variation to be missing for the MassiveZero charged current coefficient functions, as @felixhekhorn said.

Where actually they are missing for us it's hard to tell, also because there is not an actual place where the quark channel, cc, massive-zero coefficient functions are implemented, being the same of the neutral current zero-mass. I'm sorry we are not so familiar with the way APFEL has been designed, so we mainly use it as external users and only occasionally dive into the code, so if you are able to point us where the mentioned coefficient function should be placed and how the relative scale variations should be implemented we would be grateful.

Our use of APFEL (and banner inconsistencies)

The code used to produce the numbers is freely available in yadism repo, currently using the develop branch. The section involved in particular with APFEL loading, and even more in details with the Scale Variations settings, is this one that has been realized mirroring the loader provided in ccwrap (sorry for the actual runcard missing, if needed @felixhekhorn will be able to provide it asap).

However if the first part of the banner it's really related to the evolution I would not care so much, and also this may be the reason why the value is actually inconsistent: we are taking care of deactivating the evolution, indeed being yadism a pure DIS code we would like to test it against only the DIS part of APFEL.

Statement Proof

To conclude we are pretty sure on the statement that the mentioned scale variations are somehow missing in APFEL implementation, indeed we used the very same code for APFEL's call and modified our code deactivating the scale variations we are speaking about, and the results coincide perfectly, and we this setup we would have the very same discontinuity that we displayed in APFEL's numbers in the OP.

vbertone commented 4 years ago

@AleCandido The file src/DIS/IncludeScaleVariation.f takes care of including scale variation terms. Around line 304 the scale variation terms are included according to Eq. (3.17) of this document. I don't know why you don't get them right.

alecandido commented 4 years ago

Sorry @vbertone while I have nothing against 3.17, I really believe we are doing the same thing, (even if still at NLO), instead I'm not sure about the implementation, because to me it's quite hard to understand completely what's going on inside APFEL in details.

For this reason I tried to understand the content of the file you kindly pointed, and I made this sketch:

  1. first comes the ZM-VFNS (independently on the prDIS), that I would call zero-mass only, since independently on the scheme is affecting the zm coefficient functions
    • NNLO
      • in the first place the following grids are computed: P0P0, C12P0, C1LP0, C13P0 (I don't know what they are, even if I suppose them to be the ingredients of 3.17 RHS NNLO ingredients)
      • then they are used to include NNLO sv in the coefficient functions
    • NLO: it is directly compute using the splitting functions (SP at this point, I guess)
  2. then the massive ones
    • according to the perturbative order (modified only by FONLL-B), if FFNS or FONLL has been chosen for the NC some bits are added
      • indeed the NC have no NLO scale variations for the massive case
      • and at the same time CC is not available at NNLO
    • after the opposite happens: NLO is computed for CC only
  3. now there is what we are interested in: the massive-zero coefficient functions
  4. at the very end a further massive and massive-zero contributions for the NC are computed

Since what we are interested in is not zero-mass (part 1, as you said around line 304), but massive-zero (part 3, lines 586-601), and in particular the quark part, is it possible that something in there it is not correctly implemented?

On the other hand I think you can admit that the bump we pointed out it's quite a strange behavior, and it's really strange that it happens actually at the charm threshold (note: we are reporting values for Q2, but for FONLL threshold the relevant information is instead the value of muF, and this is consistent all over the place also in APFEL).

So: are you able to compute that quantity with APFEL (the points we reported) for CC, FONLL-A, muF = Q/2 (i.e. xiF = 0.5), F2charm, without having any bump? (notice that in the values we are displaying it is stressed by having used an sbaronly PDF scale independent, so I would suggest at least to use toyLH, otherwise we can give you the lhapdf files for our sbaronly PDF). Of course for having sbar -> cbar you also need to use an incoming electron/antineutrino, but as said if you would like to use the exact same configuration we will provide you the runcard (but probably the information I'm giving here are enough).

vbertone commented 4 years ago

Hi @AleCandido. Your understanding of the ZM sector of IncludeScaleVariation.f is correct. Just as a remark P0P0, C12P0, C1LP0, C13P0 are (numerical) convolutions between splitting functions and NLO coefficient functions relevant to NNLO scale variations. Also your understanding of the massive sectors is correct and I guess you agree with what is done. If not I can explain it.

It is surely possible that there is something that is not correctly implementend but the bump that you observe (that to my understanding is at the bottom threshold) has a possible explanation related to the fact that in APFEL, when you change the factorisation scale w.r.t. Q, you are also implicitly moving the heavy quark thresholds and this introduces a discontinuity at the threshold also at NLO, which does not automatically means a bug. Of course, if this is the issue there are ways to avoid it.

If you agree, I would do a benchmark using the toyLH set that is implemented in APFEL. I'll send you the numbers asap.

vbertone commented 4 years ago

@AleCandido while working on this I realised that, since you don't want to consider the evolution, using the toyLH set and doing the evolution internally with APFEL is probably not a good idea. I will use the CT18NLO LHAPDF set instead.

alecandido commented 4 years ago

@vbertone first of all: thank you for your help in the investigation.

Then, keeping the foces on the massive-zero:

(that to my understanding is at the bottom threshold) of course there are more thresholds around. But since we are speaking about the massive-zero the threshold we are mostly interested in is where FONLL kicks in, and this happens not when a specific value in Q2 is crossed, but rather a specific value in muF2 is crossed. Considering this, a value of Q2=16., with a xiF=0.5, correspond to a value of muF2=4., i.e. muF=2., and so the charm threshold, that is the first point in which FONLL is different from ZM-VFNS and everything else ("kics in"), in other words is the first point in which the massive-zero itself is active.

I agree with you that a bump is not a bug itself, but it is a quite clear clue that something is happening, and I would call it a bug if no one is able to explain that behaviour (and maybe you are, so let's keep going).


@AleCandido while working on this I realised that, since you don't want to consider the evolution, using the toyLH set and doing the evolution internally with APFEL is probably not a good idea. I will use the CT18NLO LHAPDF set instead.

For this investigation a Q2 independent PDF would be really better, so please manually install the LHAPDF file that I'm attaching (it is an sbar-only version of toyLH, in the LHAPDF format). You can load it simply specifying toy-sbaronly to LHAPDF.

toy-sbaronly.zip

alecandido commented 4 years ago

(in case you ever need it, even if I don't think so, the yadmark package of yadism project has an almost intuitive LHAPDF-PDFsets generator, for extracting specific flavors; even if the package is not deployed on PyPI it's still freely available on the repo and easily installable)

vbertone commented 4 years ago

@AleCandido here is a plot for F2charm as a function of Q at x = 10^{-3} using the toy-sbaronly.zip set (I'm really not sure why you are using this weird set with only sbar and without evolution, especially when doing scale variations). As you can see there is no discontinuity. Personally, I don't share the attitude of defining a bug something that you don't understand.

alecandido commented 4 years ago

Hi @vbertone, thank you for the plot you produced.

I'm really not sure why you are using this weird set with only sbar and without evolution, especially when doing scale variations I don't think it is weird at all, it is just the suitable object for making this test:

  • it is isolating a single contribution, decoupling the quark channel from everything else (maybe it would be clearer for you if it was a single PDF in the evolution basis, but we are mainly working in the flavour basis, and at the end of the day the main thing to decouple from it's the gluon)
  • it is not Q2 dependent because in this way we are not affected by the different evaluation point, being Q2 or muF != Q2, in this way only the new bit added to the coefficient it's relevant, again decoupling anything else

Personally, I don't share the attitude of defining a bug something that you don't understand.

It's not that we simply noticed something weird we didn't understand and we just claimed a bug. We worked for a full day on it, checking the implementation in yadism, and then, when everything from our side was explored and settled, we also took care of diving in the other side, and we tried to understand.

The problem we faced it's that a huge code written by someone else, with a not so familiar design, may be difficult to occasionally explore, so we tried to find the bug (or inconsistency) ourselves, we failed, and we simply reported that the numbers produced by APFEL are consistent with a wrong implementation of yadism, that we have under control.

If possible we would be pleased to find the bug in APFEL, because of two reasons:

So if not calling it "a bug" let's call it "a problem". Of course, I will never be sure of what it's right, but there is also a confidence level for the SM itself...

alecandido commented 4 years ago

There is also another thing I didn't consider in my previous comment: perspective.

From our perspective it's rather easy to change our code and silence a channel, but not being so familiar with APFEL this is not true for the other side. This is why we are using an uncommon PDF, in order to select some information without having to touch the code.

Looking at the mirror: probably for you it would be easier to be sure not to use the gluon channel running yadism with that PDF instead of diving in it and deactivating some channels from the code.

vbertone commented 4 years ago

Dear @AleCandido, I fully understand all you wrote and I would be glad to help you understand the reason of the discrepancy taking your perspective as much as possible. As a matter of fact, I did use the PDF set you provided me with for producing the plot I shared with you, which effectively turns off the gluon channel. Unfortunately, though, I cannot see any bug in APFEL related to factorisation-scale variations in CC F2c, as you claim. I'm not sure there is much more I can do on my side.

alecandido commented 4 years ago

I see your point @vbertone, and it's perfectly true that you did your best with the information we provided.

It's evident from your plot that you're not experiencing any discontinuity while we are, so it's undeniable that we are doing two different things. At this point it's our responsibility to provide more information to fully reproduce the run. If you would like to keep participating I would ask you just two things:

Actually doing the second should be sufficient, because we should be able to run it on our own, with your current setup and with the other one, but a further couple of eyes not conditioned by our experience it's valuable.

vbertone commented 4 years ago

Hi @AleCandido, here is the code I've used for the plot:

#include <APFEL/APFEL.h>
#include <LHAPDF/LHAPDF.h>
#include <vector>
#include <fstream>
int main()
{
  LHAPDF::PDF* pdfs = LHAPDF::mkPDF("toy-sbaronly");
  APFEL::SetMassScheme("FONLL-A");
  APFEL::SetProcessDIS("CC");
  APFEL::EnableDampingFONLL(false);
  APFEL::SetPDFSet("toy-sbaronly");
  APFEL::SetPoleMasses(pdfs->quarkThreshold(4), pdfs->quarkThreshold(5), pdfs->quarkThreshold(6));
  const double x = 0.001;
  const std::vector<double> kren = {1, 0.5};
  const int    nQ = 100;
  const double Qmin = 2;
  const double Qmax = 30;
  const double Qstp = exp( log( Qmax/ Qmin ) / ( nQ - 1 ) );
  std::ofstream fout("F2charm.dat");
  for (double k : kren)
    {
      APFEL::SetFacQRatio(k);
      APFEL::InitializeAPFEL_DIS();
      for (double Q = Qmin; Q <= Qmax; Q *= Qstp)
    {
      APFEL::SetAlphaQCDRef(pdfs->alphasQ(Q), Q);
      APFEL::ComputeStructureFunctionsAPFEL(k * Q, Q);
      fout << std::scientific << x << "\t" << Q << "\t" << APFEL::F2charm(x) << "\t" << pdfs->xfxQ(-3, x, k * Q) << "\t" << pdfs->alphasQ(Q) << std::endl;
    }
      fout << std::endl;
    }
  fout.close();
  return 0;
}

Let me know if you find any problem.

alecandido commented 4 years ago

Sorry not to reply soon (and not to be replying still).

We are still interested in the investigation, and at a first sight your code it's actually reproducing the scenario we described, but right now we need some more time for doing something else, so we will come back soon on this issue.

The only difference I spotted it's that you are correctly using the masses of the PDF I gave you, the only ones you had available, but instead we are actually using different values for the masses:

'mc': 2,
'Qmc': 2,
'kcThr': 1.0,
'mb': 4,
'Qmb': 4,
'kbThr': 1.0,
'mt': 173.07,
'Qmt': 173.07,
'ktThr': 1.0,

so this shift the thresholds, e.g.: you are using a value of mc = 1.3, so our Q2=16 jump should be placed at Q2=6.76 in your run, but it would be still in scope.

I'm still not giving you the full runcard only because the run was not on my machine and I need to retrieve it, maybe I could also reproduce from the runner we are always using, but before adding further entropy I would prefer to wait a little more (the excerpt I gave you belongs to the part we never change so I'm sure, but the most relevant information belongs of course to the other part).

vbertone commented 4 years ago

@AleCandido I think this explains the discontinuity that is the genuine FONLL discontinuity at the charm threshold that is usually tamed by introducing the damping factor (that is off here). The reason is simply that F2charm (as well as FLcharm and F3charm) in the ZM-VFN scheme is by assumption set to zero below the charm threshold and switches on right above it (according to the definition of ZM-VFNS given by Collins, Wilczek, and Zee (CWZ)). Since the scale at which PDFs are computed is muF, if you vary muF w.r.t. Q by a factor 0.5 the structure function will cross the threshold at Q = 2 mc (that with your parameter choice, unfortunately, coincides with the bottom mass, i.e. 4 GeV) and there the ZM components kicks in giving rise to the discontinuity. To conclude, APFEL seems to be behaving in the correct way, or at least in the way I planned it to behave, and thus I continue to see no bug in APFEL.

alecandido commented 4 years ago

I know that the choice of charm and bottom masses it's not the best one, we set it before start working on scale variations, and as soon as we started working on them intensively we realized it and thought to switch. For the moment we are keeping them not to break some tests, but I'm sorry for the confusion they generate.

Instead I don't understand the part on F2charm (e.g.): following the FONLL paper F2charm is defined as the heavy contribution all over the place plus a difference term. This last term is designed to be the difference between the F2light with one more light flavor (i.e. charm being light) and the massless limit of the heavy (massive-zero coefficient), in which all the power like occurrences of the mass are really set to 0 and only the log behavior is kept. Of course the log becomes unreliable below threshold, and that's the reason why you can not keep the term all over the place, but at threshold the log is exactly 0, being log(Q2/m2), so at this point the massive-zero coefficient should be equal to the light + 1, making the difference null, at least at NLO

Note that in the FONLL scheme continuity does not play a role in the matching conditions: at NLO continuity ensues accidentally from the matching conditions, but at higher orders subleading discontinuities may arise.

So I still don't see where the discontinuities should come from in FONLL-A.

alecandido commented 4 years ago

Actually I produced the same picture of you just changing the masses:

APFEL::SetPoleMasses(2., 4., 173.07);

And here is the picture: F2charm

Notice the discontinuity is only in the scale varied one, but also the other it's crossing all the thresholds (charm and bottom of course, not top).

vbertone commented 4 years ago

I think that the misunderstanding is araising from confusing the charm mass (mc) with the charm threshold. If Q = muF you are correct in saying that at NLO the log of mc / Q, that is the only term that does not cancel between zero-mass and massive-zero (btw this is an accident of NLO, beyond that there will also be constant non-power-suppressed terms that do not cancel), provides a natural damping at the threshold, where muF = Q = mc, reason why the curve with muF = Q is continuous at Q = 2 GeV. However, in general if you take muF / Q = k you will have that the ZM contribution, present in difference term of the FONLL structure function, will kick in at Q = mc / k (for the reason explained in my previous message) and there the log is no longer zero (it will be log(k)) and therefore you get a discontinuity. Notice that the fact that F2charm in FONLL at NLO with PDF and alpha_s matching at the charm mass and muF = Q is continuous at Q = mc is more of an accident than a rule. If you go to NNLO or simply displace the charm threshold the discontinuity will be there for any choice of muF and to get rid of that discontinuity you will need to introduce a damping factor.

alecandido commented 4 years ago

Probably you are right and having experienced only the NLO I'm relying too much on the continuity, even if I'm aware of that being an accident and not a rule, as you can deduce from previous FONLL paper citation.

But in this particular context (NLO CC F2/3 quark channel) the massive-zero it's not dependent at all on the log, making it continuous for both the non-scale-varied and the scale-varied bits (since the massive-zero coefficient is essentially the same as the NLO).

Maybe I'm still making an error somewhere, and this may explain the discrepancy, but if what I said it's correct the discontinuity in APFEL it's not explained.

vbertone commented 4 years ago

As I said, it is the zero-mass (ZM) (not the massive-zero) to be discontinuous at muF = k Q = mc and the discontinuity at NLO is entirely driven by log(k), which explains the continuity for k = 1.

alecandido commented 4 years ago

As I said before, it happens that for this particular case, since the log is missing, the massive-zero and the zero-mass are actually the same thing. So I don't see how it is possible to be the same thing and being continuous and discontinuous. Am I still missing something?

alecandido commented 4 years ago

In particular we realized that for the quark channel at NLO DIS FONLL it's trivial:

If you try to run the two options you will obtain different numbers, but the discrepancy it's only due to the different value of as (because of the FNS affecting as evolution). If you trick the program to use the same value of as in any way you will obtain the very same numbers.

alecandido commented 4 years ago

@vbertone can you point us where the Kf bits of Eq. (1.63) of DIS.pdf in the APFEL docs are implemented?

Looks like they should be the FONLL-A massive-zero quark-channel scale variations. And it should be the same one of the zero-mass.

vbertone commented 4 years ago

They are here for the M0 and here for the ZM, both in IncludeScaleVariation.f.

alecandido commented 4 years ago

Is it regular that SC2 is evaluated at a different point than SC3?

In particular it's the only occurrence I see of the left hand side being evaluated on Nf_FF in the second argument, and as we discussed should be the same as SC3zm.

alecandido commented 4 years ago

Same thing for massive SC3mm, but the massive-zero should be the same of the zero-mass, not of the massive.

vbertone commented 4 years ago

Sorry, I'm try to catch up. What do you mean by "Is it regular that SC2 is evaluated at a different point than SC3?"

alecandido commented 4 years ago

If you go to the pointed lines you will see:

                        SC2m0CC(igrid,jxi,3,1,beta,alpha) = 

and


                        SC3m0CC(igrid,Nf_FF,3,1,beta,alpha) = 

You see that the second argument it's not the same, and being the second argument Nf_FF it's a pretty isolated case, common only with the massive one.

vbertone commented 4 years ago

I see, I don't remember this by heart. I need to check.