DIS: CC FONLL-A scale variations

[felixhekhorn]

• We believe to have found a bug in APFEL:
• when computing F2charm/F3charm in FONLL-A scheme (i.e. NLO) for xiF!=0 there is the scale variations for the massive-zero coefficient missing
• the actual coefficient is the same as in the NC zero-mass (light) case - so it amounts to a missing 2 in the difference of the scale-variated and non-scale-variated terms (when using a fact-scale independent PDF)
• note that our result is continuous in Q2 (at ~16GeV^2), but APFEL not

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]

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

[vbertone]

@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]

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]

@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]

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
   • as before first NC NNLO
   • then CC NLO
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]

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]

@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]

@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]

(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]

@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]

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:

• we are contributing to check a code that physicists are actually using, and even if sometimes may happen that we misunderstand something, sometimes happen that we actually find a bug, so better to raise false negatives as well than remaining always silent, don't you think so?
• finding the actual source in APFEL would give us an answer, confirming that yadism it's correctly implemented, or telling us what we should fix (hopefully), and we can not be 100% sure of our implementation until the source of difference it's not completely understood

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]

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]

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]

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:

• rerun with the setup we will provide
• share the code you used for producing the numbers

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]

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.