MadGraph5_aMC@NLO

Merge lp:~maddevelopers/mg5amcnlo/RPA into lp:mg5amcnlo

  1. RPA
  2. Merge into 3.x
Proposed by Anton Safronov
Status: Needs review
Proposed branch: lp:~maddevelopers/mg5amcnlo/RPA
Merge into: lp:mg5amcnlo
Diff against target: 2049 lines (+1445/-75)
10 files modified
Template/NLO/Cards/run_card.dat (+4/-0)
Template/NLO/Source/run.inc (+6/-1)
Template/NLO/SubProcesses/fks_singular.f (+235/-25)
Template/NLO/SubProcesses/genps_fks.f (+5/-1)
Template/NLO/SubProcesses/madfks_plot.f (+194/-6)
Template/NLO/SubProcesses/setcuts.f (+5/-1)
madgraph/iolibs/export_fks.py (+28/-3)
madgraph/iolibs/template_files/parton_lum_n_fks.inc (+6/-1)
madgraph/various/banner.py (+7/-4)
madgraph/various/histograms.py (+955/-33)
To merge this branch: bzr merge lp:~maddevelopers/mg5amcnlo/RPA
Reviewer Review Type Date Requested Status
Olivier Mattelaer Pending
Review via email: mp+406325@code.launchpad.net

Commit message

I have added the pA(Ap) mode to the calculation of NLO processes.

Description of the change

I have added the new flag inside run_card, that with "True" value will run pA collision calculation.
Also, a lot of files have been changed in order to achieve this goal.
As an output one may have an HwU file, with additional crossections, and one will be able to plot RpA or RAp using this new data. Also, this code calculates uncertainties for the cross sections "on fly".

In order to test, you should just change the rpa_choice inside the run_card to "True", then to add "True" to the reweight_PDF and "False" to reweight_scale. You have two add two PDFs from lhapdf. This could be any PDF (proton-proton, proton-nucleus, nucleus-proton, nucleus-nucleus).

Then, also you have to choose Fixed-order calculation, NLO or LO, and "without parton-shower"

To post a comment you must log in.
Revision history for this message
Olivier Mattelaer (olivier-mattelaer) wrote :

Hi,

Here is my review,

I did not review histograms.py since you seems to try to force to remove all changes done in the official version (to make it py3 ready) which is not acceptable.

So I will need a second review for that (important) file. For the rest, I have made my comments and that should not be too complex to make the modifications.

Cheers,

Olivier

Reply
lp:~maddevelopers/mg5amcnlo/RPA updated
977. By Anton Safronov <email address hidden>

Added fixes according to Oliver questions, new version of Histograms.py, still problem with fks_sigular

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) wrote :

Any update here on fks_ssingular?
Do you want me to check this last revision already or do you need help on the remaining issue?

Cheers,

Olivier

Reply
Revision history for this message
Anton Safronov (safronovanton96) wrote :

> Any update here on fks_ssingular?
> Do you want me to check this last revision already or do you need help on the
> remaining issue?
>
> Cheers,
>
> Olivier

Good evening, Olivier

I think yes, I need a little help here, or a fresh opinion.
The problem is the following:

My arrays f1_p(0:100,0:MAXPROC),f2_p(0:100,0:MAXPROC),f3_p(0:100,0:MAXPROC) are multidimensional.
And I need it for the "for" loop
You're right, that "0:100" is not smart move, but the problem is, that in Fortran I can't define varibles inside the code, I should do it in the very beginning of the function. MAXPROC is global varibles, and MG5 always know what is the value. And "incontr" or "max_contr" is unknown before the line:

call weight_lines_allocated(nexternal,max_contr,iwgt,max_iproc)

This line is inside the function, and I can't normally control it from fks_singular.
I have tried to to something in "weight_lines_allocated" - but it didn't work

So, if you have any ideas - could you please share with me?

Best regards,
Anton Safronov

Reply
Revision history for this message
Olivier Mattelaer (olivier-mattelaer) wrote :

The solution is then to use allocatable memory rather than static memory allocation.
Which is part of fortran 90 convention.
Maybe the best is to contact Hua-sheng here since he has more experience than me in f90.

Cheers,

Olivier

Reply
lp:~maddevelopers/mg5amcnlo/RPA updated
978. By Anton Safronov <email address hidden>

Added allocatable varibles to fks_singular, stiil need to be checked

Revision history for this message
Anton Safronov (safronovanton96) wrote :

> The solution is then to use allocatable memory rather than static memory
> allocation.
> Which is part of fortran 90 convention.
> Maybe the best is to contact Hua-sheng here since he has more experience than
> me in f90.
>
> Cheers,
>
> Olivier

Good evening, Olivier

I have sent a letter to Hua-sheng
Mean while I have tried to do something with this, using allocatable varibles. You could find a new version of RPA program, with this fixю

And right now I have the following issure:

DEBUG: Combining HwU plots.
DEBUG: histograms.py :: =======
DEBUG: histograms.py :: Loading histograms from '/projet/pth/safronov/MG5/Anton-Process/newwithallo/SubProcesses/P0_gg_ttx/all_G1/MADatNLO.HwU'.
DEBUG: Traceback (most recent call last):
DEBUG: File "/projet/pth/safronov/MG5/Anton-Process/newwithallo/bin/internal/histograms.py", line 3996, in <module>
DEBUG: consider_reweights=consider_reweights, **file_options)
DEBUG: File "/projet/pth/safronov/MG5/Anton-Process/newwithallo/bin/internal/histograms.py", line 1833, in __init__
DEBUG: selected_central_weight=selected_label)
DEBUG: File "/projet/pth/safronov/MG5/Anton-Process/newwithallo/bin/internal/histograms.py", line 734, in __init__
DEBUG: raw_labels=raw_labels):
DEBUG: File "/projet/pth/safronov/MG5/Anton-Process/newwithallo/bin/internal/histograms.py", line 1199, in parse_one_histo_from_stream
DEBUG: " specified than expected (%i)"%len(weight_header))
DEBUG: __main__.ParseError: There is more bin weights specified than expected (200)

What is that could be?

Best regards,
Anton Safronov

Reply
Revision history for this message
Valentin Hirschi (valentin-hirschi) wrote :
Download full text (3.2 KiB)

Hi Anton,

Just a quick note regarding allocatable arrays: these are fairly simple to
use in fortran90 and you should find plenty of documentation online.
I you want to have some example of reasonably clean fortran90 code with
allocatable arrays, you can look at my DiscreteSampler implementation in:

<MG_folder>/vendor/DiscreteSampler/DiscreteSampler.f

where you can find examples of code using allocation, for instance:

```
! An allocatable to mark initialization
      logical, dimension(:), allocatable :: DS_isInitialized
[...]
          if (allocated(DS_isInitialized)) then
            write(*,*) "DiscreteSampler:: Error: The DiscreteSampler"//
     & " module can only be initialized once."
            stop 1
          else
            allocate(DS_isInitialized(1))
          endif
[...]
```

And ideally, you want to explicitly use the call to "deallocate" in order
to free the memory that you dynamically allocated with "allocate".

Regarding the histograms error you are facing, the ParseError given to you
is fairly instructive, and you should look into your HwU file being loaded
and verify that the number of weights specified in the header of this file
(i.e the first line) indeed matches the number of floating point weights
indicated for each bin entry of each histogram in that file.

I hope this helps,

Best,

On Thu, Aug 26, 2021 at 6:11 PM Anton Safronov <email address hidden>
wrote:

> > The solution is then to use allocatable memory rather than static memory
> > allocation.
> > Which is part of fortran 90 convention.
> > Maybe the best is to contact Hua-sheng here since he has more experience
> than
> > me in f90.
> >
> > Cheers,
> >
> > Olivier
>
> Good evening, Olivier
>
> I have sent a letter to Hua-sheng
> Mean while I have tried to do something with this, using allocatable
> varibles. You could find a new version of RPA program, with this fixю
>
> And right now I have the following issure:
>
> DEBUG: Combining HwU plots.
> DEBUG: histograms.py :: =======
> DEBUG: histograms.py :: Loading histograms from
> '/projet/pth/safronov/MG5/Anton-Process/newwithallo/SubProcesses/P0_gg_ttx/all_G1/MADatNLO.HwU'.
>
> DEBUG: Traceback (most recent call last):
> DEBUG: File
> "/projet/pth/safronov/MG5/Anton-Process/newwithallo/bin/internal/histograms.py",
> line 3996, in <module>
> DEBUG: consider_reweights=consider_reweights, **file_options)
> DEBUG: File
> "/projet/pth/safronov/MG5/Anton-Process/newwithallo/bin/internal/histograms.py",
> line 1833, in __init__
> DEBUG: selected_central_weight=selected_label)
> DEBUG: File
> "/projet/pth/safronov/MG5/Anton-Process/newwithallo/bin/internal/histograms.py",
> line 734, in __init__
> DEBUG: raw_labels=raw_labels):
> DEBUG: File
> "/projet/pth/safronov/MG5/Anton-Process/newwithallo/bin/internal/histograms.py",
> line 1199, in parse_one_histo_from_stream
> DEBUG: " specified than expected (%i)"%len(weight_header))
> DEBUG: __main__.ParseError: There is more bin weights specified than
> expected (200)
>
>
> What is that could be?
>
> Best regards,
> Anton Safronov
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/RPA/+merge/406325
> Your team MadDev...

Read more...

Reply
lp:~maddevelopers/mg5amcnlo/RPA updated
979. By Anton Safronov <email address hidden>

Added allocatable varibles to fks_singular, now fixed

Revision history for this message
Anton Safronov (safronovanton96) wrote :
Download full text (3.6 KiB)

> Hi Anton,
>
> Just a quick note regarding allocatable arrays: these are fairly simple to
> use in fortran90 and you should find plenty of documentation online.
> I you want to have some example of reasonably clean fortran90 code with
> allocatable arrays, you can look at my DiscreteSampler implementation in:
>
> <MG_folder>/vendor/DiscreteSampler/DiscreteSampler.f
>
> where you can find examples of code using allocation, for instance:
>
> ```
> ! An allocatable to mark initialization
> logical, dimension(:), allocatable :: DS_isInitialized
> [...]
> if (allocated(DS_isInitialized)) then
> write(*,*) "DiscreteSampler:: Error: The DiscreteSampler"//
> & " module can only be initialized once."
> stop 1
> else
> allocate(DS_isInitialized(1))
> endif
> [...]
> ```
>
> And ideally, you want to explicitly use the call to "deallocate" in order
> to free the memory that you dynamically allocated with "allocate".
>
> Regarding the histograms error you are facing, the ParseError given to you
> is fairly instructive, and you should look into your HwU file being loaded
> and verify that the number of weights specified in the header of this file
> (i.e the first line) indeed matches the number of floating point weights
> indicated for each bin entry of each histogram in that file.
>
> I hope this helps,
>
> Best,
>
> On Thu, Aug 26, 2021 at 6:11 PM Anton Safronov <email address hidden>
> wrote:
>
> > > The solution is then to use allocatable memory rather than static memory
> > > allocation.
> > > Which is part of fortran 90 convention.
> > > Maybe the best is to contact Hua-sheng here since he has more experience
> > than
> > > me in f90.
> > >
> > > Cheers,
> > >
> > > Olivier
> >
> > Good evening, Olivier
> >
> > I have sent a letter to Hua-sheng
> > Mean while I have tried to do something with this, using allocatable
> > varibles. You could find a new version of RPA program, with this fixю
> >
> > And right now I have the following issure:
> >
> > DEBUG: Combining HwU plots.
> > DEBUG: histograms.py :: =======
> > DEBUG: histograms.py :: Loading histograms from
> > '/projet/pth/safronov/MG5/Anton-
> Process/newwithallo/SubProcesses/P0_gg_ttx/all_G1/MADatNLO.HwU'.
> >
> > DEBUG: Traceback (most recent call last):
> > DEBUG: File
> > "/projet/pth/safronov/MG5/Anton-
> Process/newwithallo/bin/internal/histograms.py",
> > line 3996, in <module>
> > DEBUG: consider_reweights=consider_reweights, **file_options)
> > DEBUG: File
> > "/projet/pth/safronov/MG5/Anton-
> Process/newwithallo/bin/internal/histograms.py",
> > line 1833, in __init__
> > DEBUG: selected_central_weight=selected_label)
> > DEBUG: File
> > "/projet/pth/safronov/MG5/Anton-
> Process/newwithallo/bin/internal/histograms.py",
> > line 734, in __init__
> > DEBUG: raw_labels=raw_labels):
> > DEBUG: File
> > "/projet/pth/safronov/MG5/Anton-
> Process/newwithallo/bin/internal/histograms.py",
> > line 1199, in parse_one_histo_from_stream
> > DEBUG: " specified than expected (%i)"%len(weight_header))
> > DEBUG: __main__.ParseError: There is more bin weights specified ...

Read more...

Reply
lp:~maddevelopers/mg5amcnlo/RPA updated
980. By Anton Safronov <email address hidden>

min.fixes

Revision history for this message
Anton Safronov (safronovanton96) wrote :

Good evening, Olivier

I need your help and your experience in the question of scales
I can see the function reweight_scale in the fks_singular.f

In this function, there is no loop over PDFs
Seems like MG5 takes the central values if the first PDF and doing rescaling
Do you know, in what function or file MG5 calls this function and tells that "scaling should be done only for the first pdf"?

I have tried to add loop over PDFs by hands, and make MG5 calculate new xlum varibles and then rescale new values (similarly as in reweight_pdf), but it didn't work
I still have same values in the HwU as for the first PDF

Best regards,
Anton Safronov

Reply
Revision history for this message
Olivier Mattelaer (olivier-mattelaer) wrote :

That's not the correct place for such discussion. We do have a gitlab to discuss that.
Can you create an issue there? (and ping me or assign it to me)

Olivier

Reply
lp:~maddevelopers/mg5amcnlo/RPA updated
981. By Anton Safronov <email address hidden>

New version, with working plots, scale variation,29.11.2021

982. By Anton Safronov <email address hidden>

Merging with new version of MadGraph, 01.12.21

983. By Anton Safronov <email address hidden>

Fixed scale issure, 01.12.21

984. By Anton Safronov <email address hidden>

update, 03.03.2022

985. By Anton Safronov <email address hidden>

update, 01.04.2022. Updated the scale and PDF variation loops, now you could use any number of PDFs in the run_card. And you will get for all of them PDF variation + dcale variation. But still - check needed, if everything is working fine

986. By Anton Safronov <email address hidden>

update, 16.10.2023

Unmerged revisions

986. By Anton Safronov <email address hidden>

update, 16.10.2023

985. By Anton Safronov <email address hidden>

update, 01.04.2022. Updated the scale and PDF variation loops, now you could use any number of PDFs in the run_card. And you will get for all of them PDF variation + dcale variation. But still - check needed, if everything is working fine

984. By Anton Safronov <email address hidden>

update, 03.03.2022

983. By Anton Safronov <email address hidden>

Fixed scale issure, 01.12.21

982. By Anton Safronov <email address hidden>

Merging with new version of MadGraph, 01.12.21

981. By Anton Safronov <email address hidden>

New version, with working plots, scale variation,29.11.2021

980. By Anton Safronov <email address hidden>

min.fixes

979. By Anton Safronov <email address hidden>

Added allocatable varibles to fks_singular, now fixed

978. By Anton Safronov <email address hidden>

Added allocatable varibles to fks_singular, stiil need to be checked

977. By Anton Safronov <email address hidden>

Added fixes according to Oliver questions, new version of Histograms.py, still problem with fks_sigular

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== modified file 'Template/NLO/Cards/run_card.dat'
2--- Template/NLO/Cards/run_card.dat 2021-11-23 21:22:11 +0000
3+++ Template/NLO/Cards/run_card.dat 2023-10-16 13:26:44 +0000
4@@ -61,6 +61,10 @@
5 %(lpp2)s = lpp2 ! beam 2 type (0 = no PDF)
6 %(ebeam1)s = ebeam1 ! beam 1 energy in GeV
7 %(ebeam2)s = ebeam2 ! beam 2 energy in GeV
8+#*****************************************************
9+# RPA choice: this flag switches to the pA(Ap) mode *
10+#*****************************************************
11+ %(rpa_choice)s = rpa_choice ! choice of the regime
12 #***********************************************************************
13 # PDF choice: this automatically fixes also alpha_s(MZ) and its evol. *
14 #***********************************************************************
15
16=== modified file 'Template/NLO/Source/run.inc'
17--- Template/NLO/Source/run.inc 2020-10-19 21:03:09 +0000
18+++ Template/NLO/Source/run.inc 2023-10-16 13:26:44 +0000
19@@ -88,4 +88,9 @@
20 double precision ptmax4pdg(0:25)
21 double precision mxxmin4pdg(0:25)
22 logical mxxpart_antipart(1:25)
23- common/TO_PDG_SPECIFIC_CUT/pdg_cut, ptmin4pdg,ptmax4pdg, mxxmin4pdg, mxxpart_antipart
24\ No newline at end of file
25+ common/TO_PDG_SPECIFIC_CUT/pdg_cut, ptmin4pdg,ptmax4pdg, mxxmin4pdg, mxxpart_antipart
26+c
27+c RpA-choice
28+c
29+ logical rpa_choice
30+ common/RPA/ rpa_choice
31
32=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
33--- Template/NLO/SubProcesses/fks_singular.f 2021-12-03 15:40:10 +0000
34+++ Template/NLO/SubProcesses/fks_singular.f 2023-10-16 13:26:44 +0000
35@@ -2117,7 +2117,7 @@
36 include 'run.inc'
37 include 'timing_variables.inc'
38 include 'genps.inc'
39- integer i,kr,kf,iwgt_save,dd
40+ integer i,kr,kf,iwgt_save,dd,nn
41 double precision xlum(maxscales),dlum,pi,mu2_r(maxscales),c_mu2_r
42 $ ,c_mu2_f,mu2_f(maxscales),mu2_q,alphas,g(maxscales)
43 $ ,rwgt_muR_dep_fac,conv
44@@ -2129,14 +2129,48 @@
45 INTEGER IPROC
46 DOUBLE PRECISION PD(0:MAXPROC)
47 COMMON /SUBPROC/ PD, IPROC
48+
49+ DOUBLE PRECISION PD1(0:MAXPROC), PD2(0:MAXPROC)
50+ COMMON /PDFVALUES/ PD1, PD2
51+
52+ double precision, allocatable :: f1_p(:,:,:),f2_p(:,:,:)
53+
54+ DOUBLE PRECISION f3(0:MAXPROC)
55+
56+ DOUBLE PRECISION xlum_mod(1:nint(scalevarF(0)))
57+ integer ii,bb,k,j,jmax
58+
59+
60 parameter (conv=389379660d0) ! conversion to picobarns
61 call cpu_time(tBefore)
62- if (icontr.eq.0) return
63+
64+ do nn=1,lhaPDFid(0)
65+
66+ if (icontr.eq.0) return
67+
68+ if (nn.eq.1) then
69+ jmax=1
70+ else if (rpa_choice.eqv..true.) then
71+ jmax=3
72+ else if (rpa_choice.eqv..false.) then
73+ jmax=1
74+ endif
75+
76+
77+ do j=1,jmax
78+
79+ if (nn.eq.1.and.rpa_choice.eqv..true.) then
80+ allocate(f1_p(nint(scalevarF(0)),icontr,MAXPROC))
81+ allocate(f2_p(nint(scalevarF(0)),icontr,MAXPROC))
82+ endif
83+
84+
85 c currently we have 'iwgt' weights in the wgts() array.
86- iwgt_save=iwgt
87+ iwgt_save=iwgt
88 c loop over all the contributions in the weight lines module
89- do i=1,icontr
90+ do i=1,icontr
91 iwgt=iwgt_save
92+
93 nFKSprocess=nFKS(i)
94 xbk(1) = bjx(1,i)
95 xbk(2) = bjx(2,i)
96@@ -2151,20 +2185,94 @@
97 mu2_r(kr)=c_mu2_r*scalevarR(kr)**2
98 g(kr)=sqrt(4d0*pi*alphas(sqrt(mu2_r(kr))))
99 enddo
100+
101 c factorisation scale variation (require recomputation of the PDFs)
102 do kf=1,nint(scalevarF(0))
103 if ((.not. lscalevar(dd)) .and. kf.ne.1) exit
104+
105 mu2_f(kf)=c_mu2_f*scalevarF(kf)**2
106 q2fact(1)=mu2_f(kf)
107 q2fact(2)=mu2_f(kf)
108- xlum(kf) = dlum()
109- if (separate_flavour_configs .and. ipr(i).ne.0) then
110- if (nincoming.eq.2) then
111- xlum(kf)=pd(ipr(i))*conv
112- else
113- xlum(kf)=pd(ipr(i))
114- endif
115+
116+ if (nn.EQ.1) then! ---> central proton PDFs to be stored;
117+ xlum(kf) = dlum()
118+
119+ if (rpa_choice.eqv..true.) then
120+ do ii=1,IPROC
121+ f1_p(kf,i,ii)=PD1(ii)
122+ f2_p(kf,i,ii)=PD2(ii)
123+ enddo
124+ endif
125+
126+ if (separate_flavour_configs .and. ipr(i).ne.0) then
127+ if (nincoming.eq.2) then
128+ xlum(kf)=pd(ipr(i))*conv
129+ else
130+ xlum(kf)=pd(ipr(i))
131+ endif
132+ endif
133+ else
134+ if (j.eq.1) then
135+ xlum_mod(kf)=0D0
136+ call InitPDFm(nn,0)
137+
138+ xlum(kf) = dlum()
139+
140+ f3(0)=0
141+
142+ do ii=1,IPROC
143+ f3(ii)=PD2(ii)*PD1(ii)
144+ enddo
145+
146+ xlum(kf)=0
147+
148+ do bb=1,IPROC
149+ xlum(kf) = xlum(kf) + f3(bb)*conv
150+ enddo
151+
152+ else if (j.eq.2) then
153+
154+ xlum_mod(kf)=0D0
155+ call InitPDFm(nn,0)
156+
157+ xlum(kf) = dlum()
158+
159+ f3(0)=0
160+
161+ do ii=1,IPROC
162+ f3(ii)=PD2(ii)*f1_p(kf,i,ii)
163+ enddo
164+
165+ xlum(kf)=0
166+
167+ do bb=1,IPROC
168+ xlum(kf) = xlum(kf) + f3(bb)*conv
169+ enddo
170+
171+ else if (j.eq.3) then
172+
173+ xlum_mod(kf)=0D0
174+ call InitPDFm(nn,0)
175+
176+ xlum(kf) = dlum()
177+
178+ f3(0)=0
179+
180+ do ii=1,IPROC
181+ f3(ii)=PD1(ii)*f2_p(kf,i,ii)
182+ enddo
183+
184+ xlum(kf)=0
185+
186+ do bb=1,IPROC
187+ xlum(kf) = xlum(kf) + f3(bb)*conv
188+ enddo
189+
190+ endif
191+
192+
193 endif
194+
195 enddo
196 do kf=1,nint(scalevarF(0))
197 if ((.not. lscalevar(dd)) .and. kf.ne.1) exit
198@@ -2182,7 +2290,15 @@
199 enddo
200 enddo
201 enddo
202+ enddo
203+ enddo
204 enddo
205+
206+ if (rpa_choice.eqv..true.) then
207+ deallocate(f1_p)
208+ deallocate(f2_p)
209+ endif
210+
211 call cpu_time(tAfter)
212 tr_s=tr_s+(tAfter-tBefore)
213 return
214@@ -2306,7 +2422,7 @@
215 include 'run.inc'
216 include 'timing_variables.inc'
217 include 'genps.inc'
218- integer n,izero,i,nn
219+ integer n,izero,i,nn,counter
220 parameter (izero=0)
221 double precision xlum,dlum,pi,mu2_r,mu2_f,mu2_q,rwgt_muR_dep_fac,g
222 & ,alphas,conv
223@@ -2318,18 +2434,43 @@
224 INTEGER IPROC
225 DOUBLE PRECISION PD(0:MAXPROC)
226 COMMON /SUBPROC/ PD, IPROC
227- parameter (conv=389379660d0) ! conversion to picobarns
228- call cpu_time(tBefore)
229+
230+ DOUBLE PRECISION PD1(0:MAXPROC), PD2(0:MAXPROC)
231+ COMMON /PDFVALUES/ PD1, PD2
232+
233+ double precision, allocatable :: f1_p(:,:),f2_p(:,:)
234+ DOUBLE PRECISION f3(0:MAXPROC)
235+
236+ DOUBLE PRECISION xlum_mod(1:3)
237+ parameter (conv=389379660d0)
238+ integer jmax,j,ii,bb,k
239 if (icontr.eq.0) return
240+
241 do nn=1,lhaPDFid(0)
242-c Use as external loop the one over the PDF sets and as internal the one
243-c over the icontr. This reduces the number of calls to InitPDF and
244-c allows for better caching of the PDFs
245+
246+ if (nn.eq.1) then
247+ jmax=1
248+ else if (rpa_choice.eqv..true.) then
249+ jmax=3
250+ else if (rpa_choice.eqv..false.) then
251+ jmax=1
252+ endif
253+
254 do n=0,nmemPDF(nn)
255- iwgt=iwgt+1
256+ iwgt=iwgt+jmax
257+ counter = iwgt
258 call weight_lines_allocated(nexternal,max_contr,iwgt
259 $ ,max_iproc)
260 call InitPDFm(nn,n)
261+
262+ if (jmax==3) then
263+ iwgt=iwgt-2
264+ endif
265+
266+ if (nn.EQ.1 .and.n.EQ.0.and.rpa_choice.eqv..true.) then
267+ allocate(f1_p(icontr,MAXPROC))
268+ allocate(f2_p(icontr,MAXPROC))
269+ endif
270 do i=1,icontr
271 nFKSprocess=nFKS(i)
272 xbk(1) = bjx(1,i)
273@@ -2348,19 +2489,88 @@
274 xlum=pd(ipr(i))
275 endif
276 endif
277+ xlum_mod(2)=0D0
278+ xlum_mod(3)=0D0
279+
280+ if (nn.EQ.1 .and. n.EQ.0 .and.rpa_choice.eqv..true.) then
281+ do ii=1,IPROC
282+ f1_p(i,ii)=PD1(ii)
283+ f2_p(i,ii)=PD2(ii)
284+ enddo
285+ endif
286 c Recompute the strong coupling: alpha_s in the PDF might change
287 g=sqrt(4d0*pi*alphas(sqrt(mu2_r)))
288 c add the weights to the array
289- wgts(iwgt,i)=xlum * (wgt(1,i) + wgt(2,i)*log(mu2_r/mu2_q)
290- & +wgt(3,i)*log(mu2_f/mu2_q))*g**QCDpower(i)
291- wgts(iwgt,i)=wgts(iwgt,i)*
292- & rwgt_muR_dep_fac(sqrt(mu2_r),sqrt(mu2_r),cpower(i))
293- enddo
294- enddo
295- enddo
296+ do j=1,jmax
297+
298+ if (j==1) then! pp or AA case
299+
300+ wgts(iwgt,i)= xlum * (wgt(1,i) + wgt(2,i)*log(mu2_r/mu2_q)
301+ $ +wgt(3,i)*log(mu2_f/mu2_q))*g**QCDpower(i)
302+
303+ wgts(iwgt,i)=wgts(iwgt,i)*
304+ & rwgt_muR_dep_fac(sqrt(mu2_r),sqrt(mu2_r),cpower(i))
305+
306+ else if (j==2) then! pA case
307+ iwgt=iwgt+1
308+
309+ f3(0)=0
310+ do ii=1,IPROC
311+ f3(ii)=f1_p(i,ii)*PD2(ii)
312+ enddo
313+
314+ do bb=1,IPROC
315+ xlum_mod(2)=xlum_mod(2) + f3(bb)*conv
316+ enddo
317+
318+ wgts(iwgt,i)=xlum_mod(2) * (wgt(1,i) + wgt(2,i)*log(mu2_r/mu2_q)
319+ $ +wgt(3,i)*log(mu2_f/mu2_q))*g**QCDpower(i)
320+
321+ wgts(iwgt,i)=wgts(iwgt,i)*
322+ & rwgt_muR_dep_fac(sqrt(mu2_r),sqrt(mu2_r),cpower(i))
323+
324+ else if (j==3) then! Ap case
325+ iwgt=iwgt+1
326+
327+ f3(0)=0
328+ do ii=1,IPROC
329+ f3(ii)=PD1(ii)*f2_p(i,ii)
330+ enddo
331+
332+ do bb=1,IPROC
333+ xlum_mod(3) = xlum_mod(3) + f3(bb)*conv
334+ enddo
335+
336+ wgts(iwgt,i)=xlum_mod(3) * (wgt(1,i) + wgt(2,i)*log(mu2_r/mu2_q)
337+ $ +wgt(3,i)*log(mu2_f/mu2_q))*g**QCDpower(i)
338+
339+ wgts(iwgt,i)=wgts(iwgt,i)*
340+ & rwgt_muR_dep_fac(sqrt(mu2_r),sqrt(mu2_r),cpower(i))
341+
342+ endif
343+ enddo ! i loop
344+
345+ if (jmax.eq.3.and.icontr.gt.1) then
346+ iwgt=iwgt-2
347+ endif
348+
349+ enddo ! n loop
350+
351+ if (iwgt.ne.counter) then
352+ iwgt=counter
353+ endif
354+
355+ enddo ! j loop
356+ enddo ! nn loop
357+
358 call InitPDFm(1,0)
359 call cpu_time(tAfter)
360 tr_pdf=tr_pdf+(tAfter-tBefore)
361+
362+ if (rpa_choice.eqv..true.) then
363+ deallocate(f1_p)
364+ deallocate(f2_p)
365+ endif
366 return
367 end
368
369
370=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
371--- Template/NLO/SubProcesses/genps_fks.f 2021-11-02 08:58:23 +0000
372+++ Template/NLO/SubProcesses/genps_fks.f 2023-10-16 13:26:44 +0000
373@@ -752,7 +752,11 @@
374 enddo
375 if( firsttime .or. iconfig0.ne.iconfigsave ) then
376 if (nincoming.eq.2) then
377- stot = 4d0*ebeam(1)*ebeam(2)
378+ if (ebeam(1).le.0.938.or.ebeam(2).le.0.938) then
379+ stot=2d0*ebeam(1)*ebeam(2)
380+ else
381+ stot=4d0*ebeam(1)*ebeam(2)
382+ endif
383 else
384 stot=pmass(1)**2
385 endif
386
387=== modified file 'Template/NLO/SubProcesses/madfks_plot.f'
388--- Template/NLO/SubProcesses/madfks_plot.f 2021-05-18 11:30:50 +0000
389+++ Template/NLO/SubProcesses/madfks_plot.f 2023-10-16 13:26:44 +0000
390@@ -4,9 +4,9 @@
391 implicit none
392 include 'run.inc'
393 include "nexternal.inc"
394- integer nwgt
395- character(len=50),allocatable :: weights_info(:),ctemp(:)
396- character*13 temp
397+ integer nwgt,j,jmax
398+ character(len=100),allocatable :: weights_info(:),ctemp(:)
399+ character*100 temp
400 integer i,npdfs,ii,jj,n,kk,nn
401 double precision xsecScale_acc(maxscales,maxscales,maxdynscales)
402 $ ,xsecPDFr_acc(0:maxPDFs,maxPDFsets)
403@@ -17,7 +17,20 @@
404 nwgt=1
405 if (.not.allocated(weights_info)) allocate(weights_info(1))
406 weights_info(nwgt)="central value "
407+
408 if (do_rwgt_scale) then
409+
410+ do nn=1,lhaPDFid(0)
411+
412+ if (nn.eq.1) then
413+ jmax=1
414+ else if (nn.ne.1.and.rpa_choice.eqv..true.) then
415+ jmax=3
416+ else if (nn.ne.1.and.rpa_choice.eqv..false.) then
417+ jmax=1
418+ endif
419+
420+ if (nn.eq.1) then
421 do kk=1,dyn_scale(0)
422 c set the weights_info string for scale variations
423 if (lscalevar(kk)) then
424@@ -56,7 +69,145 @@
425 endif
426 endif
427 enddo
428+
429+
430+ else
431+
432+ do j=1,jmax
433+
434+ if (j.eq.1) then
435+
436+ do kk=1,dyn_scale(0)
437+c set the weights_info string for scale variations
438+ if (lscalevar(kk)) then
439+ do ii=1,nint(scalevarF(0))
440+ do jj=1,nint(scalevarR(0))
441+ nwgt=nwgt+1
442+ allocate(ctemp(nwgt))
443+ ctemp(1:nwgt-1)=weights_info
444+ call move_alloc(ctemp,weights_info)
445+ if (ickkw.ne.-1) then
446+ write(weights_info(nwgt),
447+ & '(a4,i4,x,a8,f6.3,x,a8,f6.3,x,i4)')
448+ $ "dyn=",dyn_scale(kk),"muR(AA)=",scalevarR(jj)
449+ $ ,"muF(AA)=",scalevarF(ii),nn
450+ else
451+ write(weights_info(nwgt),
452+ & '(a4,i4,x,a4,f6.3,x,a4,f6.3)')
453+ $ "dyn=",dyn_scale(kk),"muS=",scalevarR(jj)
454+ $ ,"muH=",scalevarF(ii)
455+ endif
456+ enddo
457+ enddo
458+ else
459+ nwgt=nwgt+1
460+ allocate(ctemp(nwgt))
461+ ctemp(1:nwgt-1)=weights_info
462+ call move_alloc(ctemp,weights_info)
463+ if (ickkw.ne.-1) then
464+ write(weights_info(nwgt),'(a4,i4,x,a8,f6.3,x,a8,f6.3,x,i4)')
465+ $ "dyn=",dyn_scale(kk),"muR(AA)=",scalevarR(1)
466+ $ ,"muF(AA)=",scalevarF(1),nn
467+ else
468+ write(weights_info(nwgt),'(a4,i4,x,a4,f6.3,x,a4,f6.3)')
469+ $ "dyn=",dyn_scale(kk),"muS=",scalevarR(1)
470+ $ ,"muH=",scalevarF(1)
471+ endif
472+ endif
473+ enddo
474+
475+ else if (j.eq.2) then
476+
477+ do kk=1,dyn_scale(0)
478+c set the weights_info string for scale variations
479+ if (lscalevar(kk)) then
480+ do ii=1,nint(scalevarF(0))
481+ do jj=1,nint(scalevarR(0))
482+ nwgt=nwgt+1
483+ allocate(ctemp(nwgt))
484+ ctemp(1:nwgt-1)=weights_info
485+ call move_alloc(ctemp,weights_info)
486+ if (ickkw.ne.-1) then
487+ write(weights_info(nwgt),
488+ & '(a4,i1,x,a8,f6.3,x,a8,f6.3,x,i4)')
489+ $ "dyn=",dyn_scale(kk),"muR(pA)=",scalevarR(jj)
490+ $ ,"muF(pA)=",scalevarF(ii),nn
491+ else
492+ write(weights_info(nwgt),
493+ & '(a4,i4,x,a4,f6.3,x,a4,f6.3)')
494+ $ "dyn=",dyn_scale(kk),"muS=",scalevarR(jj)
495+ $ ,"muH=",scalevarF(ii)
496+ endif
497+ enddo
498+ enddo
499+ else
500+ nwgt=nwgt+1
501+ allocate(ctemp(nwgt))
502+ ctemp(1:nwgt-1)=weights_info
503+ call move_alloc(ctemp,weights_info)
504+ if (ickkw.ne.-1) then
505+ write(weights_info(nwgt),'(a4,i1,x,a8,f6.3,x,a8,f6.3,x,i4)')
506+ $ "dyn=",dyn_scale(kk),"muR(pA)=",scalevarR(1)
507+ $ ,"muF(pA)=",scalevarF(1),nn
508+ else
509+ write(weights_info(nwgt),'(a4,i4,x,a4,f6.3,x,a4,f6.3)')
510+ $ "dyn=",dyn_scale(kk),"muS=",scalevarR(1)
511+ $ ,"muH=",scalevarF(1)
512+ endif
513+ endif
514+ enddo
515+
516+ else if (j.eq.3) then
517+
518+ do kk=1,dyn_scale(0)
519+c set the weights_info string for scale variations
520+ if (lscalevar(kk)) then
521+ do ii=1,nint(scalevarF(0))
522+ do jj=1,nint(scalevarR(0))
523+ nwgt=nwgt+1
524+ allocate(ctemp(nwgt))
525+ ctemp(1:nwgt-1)=weights_info
526+ call move_alloc(ctemp,weights_info)
527+ if (ickkw.ne.-1) then
528+ write(weights_info(nwgt),
529+ & '(a4,i4,x,a8,f6.3,x,a8,f6.3,x,i4)')
530+ $ "dyn=",dyn_scale(kk),"muR(Ap)=",scalevarR(jj)
531+ $ ,"muF(Ap)=",scalevarF(ii),nn
532+ else
533+ write(weights_info(nwgt),
534+ & '(a4,i4,x,a4,f6.3,x,a4,f6.3)')
535+ $ "dyn=",dyn_scale(kk),"muS=",scalevarR(jj)
536+ $ ,"muH=",scalevarF(ii)
537+ endif
538+ enddo
539+ enddo
540+ else
541+ nwgt=nwgt+1
542+ allocate(ctemp(nwgt))
543+ ctemp(1:nwgt-1)=weights_info
544+ call move_alloc(ctemp,weights_info)
545+ if (ickkw.ne.-1) then
546+ write(weights_info(nwgt),'(a4,i4,x,a8,f6.3,x,a8,f6.3,x,i4)')
547+ $ "dyn=",dyn_scale(kk),"muR(Ap)=",scalevarR(1)
548+ $ ,"muF(Ap)=",scalevarF(1),nn
549+ else
550+ write(weights_info(nwgt),'(a4,i4,x,a4,f6.3,x,a4,f6.3)')
551+ $ "dyn=",dyn_scale(kk),"muS=",scalevarR(1)
552+ $ ,"muH=",scalevarF(1)
553+ endif
554+ endif
555+ enddo
556+
557+ endif
558+ enddo
559+ endif
560+
561+ enddo
562+
563 endif
564+
565+
566+
567 if (do_rwgt_pdf) then
568 do nn=1,lhaPDFid(0)
569 if (lpdfvar(nn)) then
570@@ -87,15 +238,51 @@
571 endif
572 c set the weights_info string for PDF variation
573 if (lpdfvar(nn)) then
574+ if (nn.eq.1) then
575+ jmax=1
576+ else if (nn.ne.1.and.rpa_choice.eqv..true.) then
577+ jmax=3
578+ else if (nn.ne.1.and.rpa_choice.eqv..false.) then
579+ jmax=1
580+ endif
581+
582+
583+ !do j=1,jmax
584+
585 do n=0,nmemPDF(nn)
586+
587+ do j=1,jmax
588 nwgt=nwgt+1
589 allocate(ctemp(nwgt))
590 ctemp(1:nwgt-1)=weights_info
591 call move_alloc(ctemp,weights_info)
592+
593+ if (j==1) then !pp
594+
595 write(temp,'(a4,i8)') "PDF=",lhaPDFid(nn)+n
596 write(weights_info(nwgt),'(a)') trim(adjustl(temp))/
597- $ /' '//trim(adjustl(lhaPDFsetname(nn)))
598- enddo
599+ $ /' '//trim(adjustl(lhaPDFsetname(nn)))
600+ !enddo
601+
602+ else if (j==2) then !pA
603+
604+ write(temp,'(a4,i8)') "pA=",lhaPDFid(nn)+n !Here we fill HwU file with pA crossections only if we have the right flag in run_card dat
605+ write(weights_info(nwgt),'(a)') trim(adjustl(temp))/
606+ $ /' '//trim(adjustl(lhaPDFsetname(1)))//' with '/
607+ $ /trim(adjustl(lhaPDFsetname(nn)))
608+
609+ else if (j==3) then !Ap
610+
611+ write(temp,'(a4,i8)') "Ap=",lhaPDFid(nn)+n !Here we fill HwU file with Ap crossections only if we have the right flag in run_card dat
612+ write(weights_info(nwgt),'(a)') trim(adjustl(temp))/
613+ $ /' '//trim(adjustl(lhaPDFsetname(nn)))//' with '/
614+ $ /trim(adjustl(lhaPDFsetname(1)))
615+
616+ endif
617+ enddo
618+ enddo
619+
620+
621 else
622 nwgt=nwgt+1
623 allocate(ctemp(nwgt))
624@@ -104,7 +291,8 @@
625 write(temp,'(a4,i8)') "PDF=",lhaPDFid(nn)
626 write(weights_info(nwgt),'(a)') trim(adjustl(temp))/
627 $ /' '//trim(adjustl(lhaPDFsetname(nn)))
628- endif
629+
630+ endif
631 enddo
632 c start with central member of the first set
633 call InitPDFm(1,0)
634
635=== modified file 'Template/NLO/SubProcesses/setcuts.f'
636--- Template/NLO/SubProcesses/setcuts.f 2021-09-09 15:51:40 +0000
637+++ Template/NLO/SubProcesses/setcuts.f 2023-10-16 13:26:44 +0000
638@@ -416,7 +416,11 @@
639 endif
640 xw(i)=0d0
641 enddo
642- stot = 4d0*ebeam(1)*ebeam(2)
643+ if (ebeam(1).le.0.938.or.ebeam(2).le.0.938) then
644+ stot=2d0*ebeam(1)*ebeam(2)
645+ else
646+ stot=4d0*ebeam(1)*ebeam(2)
647+ endif
648 tau_Born_lower_bound=taumin(iFKS,ichan)**2/stot
649 tau_lower_bound=taumin_j(iFKS,ichan)**2/stot
650 c
651
652=== modified file 'madgraph/iolibs/export_fks.py'
653--- madgraph/iolibs/export_fks.py 2022-04-15 10:09:24 +0000
654+++ madgraph/iolibs/export_fks.py 2023-10-16 13:26:44 +0000
655@@ -3749,9 +3749,7 @@
656 ("c settings other partons flavours outside quark, gluon, photon to 0d0\n" + \
657 "%s%d=0d0\n") % \
658 (pdf_codes[initial_state],i + 1)
659-
660 pdf_lines = pdf_lines + "ENDIF\n"
661-
662 # Add up PDFs for the different initial state particles
663 pdf_lines = pdf_lines + "PD(0) = 0d0\nIPROC = 0\n"
664 for proc in processes:
665@@ -3771,17 +3769,44 @@
666
667 # Remove last "*" from pdf_lines
668 pdf_lines = pdf_lines[:-1] + "\n"
669-
670+
671 # this is for the lepton collisions with electron luminosity
672 # put here "%s%d_components(i_ee)*%s%d_components(i_ee)"
673 pdf_lines += "if (ABS(LPP(1)).EQ.ABS(LPP(2)).and. (ABS(LPP(1)).EQ.3.or.ABS(LPP(1)).EQ.4).and.pdlabel.ne.'none')" + \
674 "PD(IPROC)=ee_comp_prod(%s_components,%s_components)\n" % \
675 tuple(comp_list)
676
677+ pdf_lines = pdf_lines + "\nPD2(0) = 0d0\nIPROS = 0\n"
678+ for proc in processes:
679+ process_line = proc.base_string()
680+ pdf_lines = pdf_lines + "\nIPROS=IPROS+1 ! " + process_line
681+ pdf_lines = pdf_lines + "\nPD2(IPROS) = "
682+ for ibeam in [2]:
683+ initial_state = proc.get_initial_pdg(ibeam)
684+ if initial_state in pdf_codes.keys():
685+ pdf_lines = pdf_lines + "%s%d" % \
686+ (pdf_codes[initial_state], ibeam) + "\n"
687+ else:
688+ pdf_lines = pdf_lines + "1d0*"
689+ pdf_lines = pdf_lines + "\nPD1(0) = 0d0\nIPROSS = 0\n"
690+ for proc in processes:
691+ process_line = proc.base_string()
692+ pdf_lines = pdf_lines + "\nIPROSS=IPROSS+1 ! " + process_line
693+ pdf_lines = pdf_lines + "\nPD1(IPROSS) = "
694+ for ibeam in [1]:
695+ initial_state = proc.get_initial_pdg(ibeam)
696+ if initial_state in pdf_codes.keys():
697+ pdf_lines = pdf_lines + "%s%d" % \
698+ (pdf_codes[initial_state], ibeam) + "\n"
699+ #pdf_lines = pdf_lines[:-1] + "\n"
700+ else:
701+ pdf_lines = pdf_lines + "1d0*"
702+
703 # Remove last line break from pdf_lines
704 return pdf_definition_lines[:-1], pdf_data_lines[:-1], pdf_lines[:-1], ee_pdf_definition_lines
705
706
707+
708 #test written
709 def get_color_data_lines_from_color_matrix(self, color_matrix, n=6):
710 """Return the color matrix definition lines for the given color_matrix. Split
711
712=== modified file 'madgraph/iolibs/template_files/parton_lum_n_fks.inc'
713--- madgraph/iolibs/template_files/parton_lum_n_fks.inc 2021-04-16 18:57:14 +0000
714+++ madgraph/iolibs/template_files/parton_lum_n_fks.inc 2023-10-16 13:26:44 +0000
715@@ -32,7 +32,12 @@
716 C
717 INTEGER IPROC
718 DOUBLE PRECISION PD(0:MAXPROC)
719- COMMON /SubProc/ PD, IPROC
720+ COMMON /SubProc/ PD, IPROC
721+
722+ INTEGER IPROS, IPROSS
723+ DOUBLE PRECISION PD1(0:MAXPROC), PD2(0:MAXPROC)
724+ COMMON /PDFVALUES/ PD1, PD2
725+
726 include "coupl.inc"
727 include "run.inc"
728 integer imirror
729
730=== modified file 'madgraph/various/banner.py'
731--- madgraph/various/banner.py 2022-04-15 11:56:56 +0000
732+++ madgraph/various/banner.py 2023-10-16 13:26:44 +0000
733@@ -4821,9 +4821,12 @@
734 self.add_param('lpp1', 1, fortran_name='lpp(1)')
735 self.add_param('lpp2', 1, fortran_name='lpp(2)')
736 self.add_param('ebeam1', 6500.0, fortran_name='ebeam(1)')
737- self.add_param('ebeam2', 6500.0, fortran_name='ebeam(2)')
738+ self.add_param('ebeam2', 6500.0, fortran_name='ebeam(2)')
739+##############Anton
740+ self.add_param('rpa_choice', False, fortran_name='rpa_choice',comment='some help here')
741+##############Anton
742 self.add_param('pdlabel', 'nn23nlo', allowed=['lhapdf', 'cteq6_m','cteq6_d','cteq6_l','cteq6l1', 'nn23lo','nn23lo1','nn23nlo','ct14q00','ct14q07','ct14q14','ct14q21'] +\
743- sum(self.allowed_lep_densities.values(),[]) )
744+ sum(self.allowed_lep_densities.values(),[]) )
745 self.add_param('lhaid', [244600],fortran_name='lhaPDFid')
746 self.add_param('lhapdfsetname', ['internal_use_only'], system=True)
747 #shower and scale
748@@ -4914,8 +4917,8 @@
749
750 # for lepton-lepton collisions, ignore 'pdlabel' and 'lhaid'
751 if abs(self['lpp1'])!=1 or abs(self['lpp2'])!=1:
752- if self['lpp1'] == 1 or self['lpp2']==1:
753- raise InvalidRunCard('Process like Deep Inelastic scattering not supported at NLO accuracy.')
754+ #if self['lpp1'] == 1 or self['lpp2']==1:
755+ #raise InvalidRunCard('Process like Deep Inelastic scattering not supported at NLO accuracy.')
756
757 if self['lpp1'] == 0 == self['lpp2']:
758 if self['pdlabel']!='nn23nlo' or self['reweight_pdf']:
759
760=== modified file 'madgraph/various/histograms.py'
761--- madgraph/various/histograms.py 2022-04-04 19:27:18 +0000
762+++ madgraph/various/histograms.py 2023-10-16 13:26:44 +0000
763@@ -662,6 +662,24 @@
764 '\s*muf\s*=\s*(?P<muf_fact>%s)\s*$'%a_float_re,re.IGNORECASE)
765 weight_label_PDF_adv = re.compile('^\s*PDF\s*=\s*(?P<PDF_set>\d+)\s+(?P<PDF_set_cen>\S+)\s*$')
766
767+############S.A
768+############ We read two additional new flags for pA and Ap collisions
769+ weight_label_PDF_pA = re.compile('^\s*pA\s*=\s*(?P<PDF_set1>\d+)\s+(?P<PDF_set_1>\S+)\s+\s*with\s+(?P<PDF_set_2>\S+)\s*$')
770+
771+ weight_label_PDF_Ap = re.compile('^\s*Ap\s*=\s*(?P<PDF_set2>\d+)\s+(?P<PDF_set_3>\S+)\s+\s*with\s+(?P<PDF_set_4>\S+)\s*$')
772+
773+ weight_label_scale_pA_adv = re.compile('^\s*dyn\s*=\s*(?P<dyn_choice1>%s)'%a_int_re+\
774+ '\s*muR\(pA\)\s*=\s*(?P<mur_fact1>%s)'%a_float_re+\
775+ '\s*muF\(pA\)\s*=\s*(?P<muf_fact1>%s)\s*$'%a_float_re,re.IGNORECASE)
776+
777+ weight_label_scale_Ap_adv = re.compile('^\s*dyn\s*=\s*(?P<dyn_choice2>%s)'%a_int_re+\
778+ '\s*muR\(Ap\)\s*=\s*(?P<mur_fact2>%s)'%a_float_re+\
779+ '\s*muF\(Ap\)\s*=\s*(?P<muf_fact2>%s)\s*$'%a_float_re,re.IGNORECASE)
780+
781+
782+
783+############S.A
784+
785
786 class ParseError(MadGraph5Error):
787 """a class for histogram data parsing errors"""
788@@ -906,6 +924,23 @@
789 others.append('PDF=%i %s'%(label[1],label[2]))
790 elif label_type == 'alpsfact':
791 others.append('alpsfact=%d'%label[1])
792+
793+#############S.A
794+#############Creating additional two labels
795+ #elif label_type == 'pdf_pA':
796+ #others.append('PDF=%i %s multiply by %s'%(label[1],label[2],label[3]))
797+ elif label_type == 'pdf_pA':
798+ others.append('pA=%i %s with %s'%(label[1],label[2],label[3]))
799+
800+ elif label_type == 'pdf_Ap':
801+ others.append('Ap=%i %s with %s'%(label[1],label[2],label[3]))
802+
803+ elif label_type == 'scale_adv_pA':
804+ others.append('dynam=%i muR(pA)=%6.3f muF(pA)=%6.3f'%(label[1],label[2],label[3]))
805+
806+ elif label_type == 'scale_adv_Ap':
807+ others.append('dynam=%i muR(Ap)=%6.3f muF(Ap)=%6.3f'%(label[1],label[2],label[3]))
808+#############S.A
809
810 return res+' & '.join(others)
811
812@@ -1022,6 +1057,13 @@
813 alpsfact_wgt = HwU.weight_label_alpsfact.match(h)
814 scale_wgt_adv = HwU.weight_label_scale_adv.match(h)
815 PDF_wgt_adv = HwU.weight_label_PDF_adv.match(h)
816+#######S.A
817+ PDF_wgt_pA = HwU.weight_label_PDF_pA.match(h)
818+ PDF_wgt_Ap = HwU.weight_label_PDF_Ap.match(h)
819+ scale_wgt_pA = HwU.weight_label_scale_pA_adv.match(h)
820+ scale_wgt_Ap = HwU.weight_label_scale_Ap_adv.match(h)
821+#######S.A
822+
823 if scale_wgt_adv:
824 header[i] = ('scale_adv',
825 int(scale_wgt_adv.group('dyn_choice')),
826@@ -1040,7 +1082,33 @@
827 elif Merging_wgt:
828 header[i] = ('merging_scale',float(Merging_wgt.group('Merging_scale')))
829 elif alpsfact_wgt:
830- header[i] = ('alpsfact',float(alpsfact_wgt.group('alpsfact')))
831+ header[i] = ('alpsfact',float(alpsfact_wgt.group('alpsfact')))
832+
833+##########S.A
834+##########Creating additional two arrys of data, that will be inside our new HwU files
835+ elif PDF_wgt_pA:
836+ header[i] = ('pdf_pA',
837+ int(PDF_wgt_pA.group('PDF_set1')),
838+ PDF_wgt_pA.group('PDF_set_1'),
839+ PDF_wgt_pA.group('PDF_set_2'))
840+ elif PDF_wgt_Ap:
841+ header[i] = ('pdf_Ap',
842+ int(PDF_wgt_Ap.group('PDF_set2')),
843+ PDF_wgt_Ap.group('PDF_set_3'),
844+ PDF_wgt_Ap.group('PDF_set_4'))
845+
846+ elif scale_wgt_pA:
847+ header[i] = ('scale_adv_pA',
848+ int(scale_wgt_pA.group('dyn_choice1')),
849+ float(scale_wgt_pA.group('mur_fact1')),
850+ float(scale_wgt_pA.group('muf_fact1')))
851+
852+ elif scale_wgt_Ap:
853+ header[i] = ('scale_adv_Ap',
854+ int(scale_wgt_Ap.group('dyn_choice2')),
855+ float(scale_wgt_Ap.group('mur_fact2')),
856+ float(scale_wgt_Ap.group('muf_fact2')))
857+##########S.A
858
859 return header
860
861@@ -1224,7 +1292,19 @@
862 scale_position = -1
863 elif type.upper()=='PDF':
864 new_wgt_label = 'delta_pdf'
865- scale_position = -2
866+ scale_position = -2
867+
868+#################S.A
869+
870+ elif type.upper()=='RPA':
871+ new_wgt_label = 'delta_pA'
872+ scale_position = -3
873+ elif type.upper()=='RAP':
874+ new_wgt_label = 'delta_Ap'
875+ scale_position = -4
876+
877+#################S.A
878+
879 elif type.upper()=='MERGING':
880 new_wgt_label = 'delta_merging'
881 elif type.upper()=='ALPSFACT':
882@@ -1298,13 +1378,46 @@
883 if wgts:
884 wgts_to_consider.append(wgts)
885 label_to_consider.append('none')
886+
887+##############S.A
888+ elif scale_position == -3:
889+ pdf_sets=[label[3] for label in self.bins.weight_labels if \
890+ HwU.get_HwU_wgt_label_type(label)=='pdf_pA']
891+
892+
893+ # remove doubles in list but keep the order!
894+ pdf_sets=[ii for n,ii in enumerate(pdf_sets) if ii not in pdf_sets[:n]]
895+ for pdf_set in pdf_sets:
896+ wgts=[label for label in self.bins.weight_labels if \
897+ HwU.get_HwU_wgt_label_type(label)=='pdf_pA' and label[3]==pdf_set]
898+ if wgts:
899+ wgts_to_consider.append(wgts)
900+ label_to_consider.append(pdf_set)
901+
902+ elif scale_position == -4:
903+ pdf_sets=[label[2] for label in self.bins.weight_labels if \
904+ HwU.get_HwU_wgt_label_type(label)=='pdf_Ap']
905+
906+
907+ # remove doubles in list but keep the order!
908+ pdf_sets=[ii for n,ii in enumerate(pdf_sets) if ii not in pdf_sets[:n]]
909+ for pdf_set in pdf_sets:
910+ wgts=[label for label in self.bins.weight_labels if \
911+ HwU.get_HwU_wgt_label_type(label)=='pdf_Ap' and label[2]==pdf_set]
912+ if wgts:
913+ wgts_to_consider.append(wgts)
914+ label_to_consider.append(pdf_set)
915+
916+##############S.A
917+
918
919 if len(wgts_to_consider)==0 or all(len(wgts)==0 for wgts in wgts_to_consider):
920 # No envelope can be constructed, it is not worth adding the weights
921 return (None,[None])
922
923 # find and import python version of lhapdf if doing PDF uncertainties
924- if type=='PDF':
925+################################ S.A
926+ if type=='PDF' or type=='RPA' or type=='RAP':
927 use_lhapdf=False
928 try:
929 lhapdf_libdir=subprocess.Popen([lhapdfconfig,'--libdir'],\
930@@ -1362,7 +1475,15 @@
931
932 if type=='PDF' and use_lhapdf:
933 lhapdf.setVerbosity(0)
934-
935+
936+###################S.A
937+ if type=='RPA' and use_lhapdf:
938+ lhapdf.setVerbosity(0)
939+
940+ if type=='RAP' and use_lhapdf:
941+ lhapdf.setVerbosity(0)
942+
943+##################S.A
944 # Place the new weight label last before the first tuple
945 position=[]
946 labels=[]
947@@ -1391,27 +1512,61 @@
948
949 if type=='PDF' and use_lhapdf and label != 'none':
950 p=lhapdf.getPDFSet(label)
951-
952+
953+###################S.A
954+ if type=='RPA' and use_lhapdf and label != 'none':
955+ p=lhapdf.getPDFSet(label)
956+
957+ if type=='RAP' and use_lhapdf and label != 'none':
958+ p=lhapdf.getPDFSet(label)
959+
960+###################S.A
961 # Now add the corresponding weight to all Bins
962 for bin in self.bins:
963- if type!='PDF':
964+###################S.A
965+ if type!='PDF' and type!='RPA' and type!='RAP':
966+###################S.A
967+
968 bin.wgts[new_wgt_labels[0]] = bin.wgts[wgts[0]]
969 bin.wgts[new_wgt_labels[1]] = min(bin.wgts[label] \
970 for label in wgts)
971 bin.wgts[new_wgt_labels[2]] = max(bin.wgts[label] \
972 for label in wgts)
973+ #a=bin.wgts[new_wgt_labels[0]]
974+ #f1.write(str(wgts))
975+ #f1.write('\n')
976+
977+###################S.A
978+ elif type=='RPA' and use_lhapdf and label != 'none' and len(wgts) > 1:
979+ pdfs = [bin.wgts[rpa] for rpa in sorted(wgts)]
980+ ep=p.uncertainty(pdfs,-1)
981+
982+ bin.wgts[new_wgt_labels[0]] = ep.central
983+ bin.wgts[new_wgt_labels[1]] = ep.central-ep.errminus
984+ bin.wgts[new_wgt_labels[2]] = ep.central+ep.errplus
985+
986+ elif type=='RAP' and use_lhapdf and label != 'none' and len(wgts) > 1:
987+ pdfs = [bin.wgts[rap] for rap in sorted(wgts)]
988+ ep=p.uncertainty(pdfs,-1)
989+ bin.wgts[new_wgt_labels[0]] = ep.central
990+ bin.wgts[new_wgt_labels[1]] = ep.central-ep.errminus
991+ bin.wgts[new_wgt_labels[2]] = ep.central+ep.errplus
992+
993+###################S.A
994+
995 elif type=='PDF' and use_lhapdf and label != 'none' and len(wgts) > 1:
996- pdfs = [bin.wgts[pdf] for pdf in sorted(wgts)]
997+ pdfs = [bin.wgts[pdf] for pdf in sorted(wgts)]
998 ep=p.uncertainty(pdfs,-1)
999 bin.wgts[new_wgt_labels[0]] = ep.central
1000 bin.wgts[new_wgt_labels[1]] = ep.central-ep.errminus
1001 bin.wgts[new_wgt_labels[2]] = ep.central+ep.errplus
1002+
1003 elif type=='PDF' and use_lhapdf and label != 'none' and len(bin.wgts) == 1:
1004 bin.wgts[new_wgt_labels[0]] = bin.wgts[wgts[0]]
1005 bin.wgts[new_wgt_labels[1]] = bin.wgts[wgts[0]]
1006 bin.wgts[new_wgt_labels[2]] = bin.wgts[wgts[0]]
1007 else:
1008- pdfs = [bin.wgts[pdf] for pdf in sorted(wgts)]
1009+ pdfs = [bin.wgts[pdf] for pdf in sorted(wgts)]
1010 pdf_up = 0.0
1011 pdf_down = 0.0
1012 cntrl_val = bin.wgts['central']
1013@@ -1845,12 +2000,25 @@
1014 # and 'MERGING' defined. If absent we specify '-1' which implies that
1015 # the 'default' value was used (whatever it was).
1016 # Also cast them in the proper type
1017- for wgt_label in all_weights:
1018- for mandatory_attribute in ['PDF','MUR','MUF','MERGING','ALPSFACT']:
1019+ for wgt_label in all_weights:
1020+ for mandatory_attribute in ['PDF','MUR','MUF','MERGING','ALPSFACT','RPA','RAP']:##############S.A
1021 if mandatory_attribute not in wgt_label:
1022 wgt_label[mandatory_attribute] = '-1'
1023+
1024+##############S.A
1025+
1026 if mandatory_attribute=='PDF':
1027 wgt_label[mandatory_attribute] = int(wgt_label[mandatory_attribute])
1028+
1029+##############S.A
1030+ if mandatory_attribute=='RPA':
1031+ wgt_label[mandatory_attribute] = int(wgt_label[mandatory_attribute])
1032+
1033+ if mandatory_attribute=='RAP':
1034+ wgt_label[mandatory_attribute] = int(wgt_label[mandatory_attribute])
1035+
1036+##############S.A
1037+
1038 elif mandatory_attribute in ['MUR','MUF','MERGING','ALPSFACT']:
1039 wgt_label[mandatory_attribute] = float(wgt_label[mandatory_attribute])
1040
1041@@ -1861,9 +2029,17 @@
1042 merging_scale_chosen = all_weights[2]['MERGING']
1043 else:
1044 merging_scale_chosen = merging_scale
1045-
1046+
1047+ #central_pA_scale = all_weights[2]['SCALE_PA']
1048+
1049 # Central weight parameters are enforced to be those of the third weight
1050 central_PDF = all_weights[2]['PDF']
1051+
1052+##############S.A
1053+ central_RPA = all_weights[2]['RPA']
1054+ central_RAP = all_weights[2]['RAP']
1055+##############S.A
1056+
1057 # Assume central scale is one, unless specified.
1058 central_MUR = all_weights[2]['MUR'] if all_weights[2]['MUR']!=-1.0 else 1.0
1059 central_MUF = all_weights[2]['MUF'] if all_weights[2]['MUF']!=-1.0 else 1.0
1060@@ -1897,8 +2073,18 @@
1061 if weight['MUR'] not in [central_MUR, -1.0] or \
1062 weight['MUF'] not in [central_MUF, -1.0]:
1063 differences.append('mur_muf_scale')
1064+
1065+
1066 if weight['PDF'] not in [central_PDF,-1]:
1067 differences.append('pdf')
1068+
1069+##############S.A
1070+ if weight['RPA'] not in [central_RPA,-1]:
1071+ differences.append('rpa')
1072+ if weight['RAP'] not in [central_RAP,-1]:
1073+ differences.append('rap')
1074+##############S.A
1075+
1076 if weight['ALPSFACT'] not in [central_alpsfact, -1]:
1077 differences.append('ALPSFACT')
1078 return set(differences)
1079@@ -1918,8 +2104,8 @@
1080 all_properties = [property for property in all_properties if
1081 not weight[property] is None]
1082
1083- # then add PDF, MUR, MUF and MERGING if present
1084- for property in ['PDF','MUR','MUF','ALPSFACT','MERGING']:
1085+ # then add PDF, MUR, MUF and MERGING if present
1086+ for property in ['PDF','MUR','MUF','ALPSFACT','MERGING','RPA','RAP']:#########S.A
1087 all_properties.pop(all_properties.index(property))
1088 if weight[property]!=-1:
1089 ordered_properties.append(property)
1090@@ -1978,8 +2164,19 @@
1091 wgt_label = ('scale',weight['MUR'],weight['MUF'])
1092 if variations in [set(['ALPSFACT']),set(['pdf','ALPSFACT'])]:
1093 wgt_label = ('alpsfact',weight['ALPSFACT'])
1094+
1095+############S.A
1096+
1097 if variations == set(['pdf']):
1098- wgt_label = ('pdf',weight['PDF'])
1099+ wgt_label = ('pdf',weight['PDF'])
1100+
1101+############S.A
1102+ if variations == set(['rpa']):
1103+ wgt_label = ('rpa',weight['RPA'])
1104+ if variations == set(['rap']):
1105+ wgt_label = ('rap',weight['RAP'])
1106+############S.A
1107+
1108 if variations == set([]):
1109 # Unknown weight (might turn out to be taken as a merging variation weight below)
1110 wgt_label = format_weight_label(weight)
1111@@ -2239,25 +2436,25 @@
1112 # color-blind individuals with either protanopia or deuteranopia. Bang
1113 # Wong [2011] Nature Methods 8, 441.
1114
1115-set style line 1 lt 1 lc rgb "%(col1)s" lw 2.5
1116+set style line 1 lt 1 lc rgb "blue" lw 2.5
1117 set style line 11 lt 2 lc rgb "%(col1)s" lw 2.5
1118 set style line 21 lt 4 lc rgb "%(col1)s" lw 2.5
1119 set style line 31 lt 6 lc rgb "%(col1)s" lw 2.5
1120 set style line 41 lt 8 lc rgb "%(col1)s" lw 2.5
1121
1122-set style line 2 lt 1 lc rgb "%(col2)s" lw 2.5
1123+set style line 2 lt 1 lc rgb "red" lw 2.5
1124 set style line 12 lt 2 lc rgb "%(col2)s" lw 2.5
1125 set style line 22 lt 4 lc rgb "%(col2)s" lw 2.5
1126 set style line 32 lt 6 lc rgb "%(col2)s" lw 2.5
1127 set style line 42 lt 8 lc rgb "%(col2)s" lw 2.5
1128
1129-set style line 3 lt 1 lc rgb "%(col3)s" lw 2.5
1130+set style line 3 lt 1 lc rgb "green" lw 2.5
1131 set style line 13 lt 2 lc rgb "%(col3)s" lw 2.5
1132 set style line 23 lt 4 lc rgb "%(col3)s" lw 2.5
1133 set style line 33 lt 6 lc rgb "%(col3)s" lw 2.5
1134 set style line 43 lt 8 lc rgb "%(col3)s" lw 2.5
1135
1136-set style line 4 lt 1 lc rgb "%(col4)s" lw 2.5
1137+set style line 4 lt 1 lc rgb "purple" lw 2.5
1138 set style line 14 lt 2 lc rgb "%(col4)s" lw 2.5
1139 set style line 24 lt 4 lc rgb "%(col4)s" lw 2.5
1140 set style line 34 lt 6 lc rgb "%(col4)s" lw 2.5
1141@@ -2327,28 +2524,28 @@
1142 # color-blind individuals with either protanopia or deuteranopia. Bang
1143 # Wong [2011] Nature Methods 8, 441.
1144
1145-set style line 1 lt 1 lc rgb "%(col1)s" lw 1.3
1146+set style line 1 lt 1 lc rgb "blue" lw 1.7
1147 set style line 101 lt 1 lc rgb "%(col1)s" lw 1.3 dt (6,3)
1148 set style line 11 lt 2 lc rgb "%(col1)s" lw 1.3 dt (6,3)
1149 set style line 21 lt 4 lc rgb "%(col1)s" lw 1.3 dt (3,2)
1150 set style line 31 lt 6 lc rgb "%(col1)s" lw 1.3 dt (2,1)
1151 set style line 41 lt 8 lc rgb "%(col1)s" lw 1.3 dt (4,3)
1152
1153-set style line 2 lt 1 lc rgb "%(col2)s" lw 1.3
1154+set style line 2 lt 1 lc rgb "red" lw 1.3
1155 set style line 102 lt 1 lc rgb "%(col2)s" lw 1.3 dt (6,3)
1156 set style line 12 lt 2 lc rgb "%(col2)s" lw 1.3 dt (6,3)
1157 set style line 22 lt 4 lc rgb "%(col2)s" lw 1.3 dt (3,2)
1158 set style line 32 lt 6 lc rgb "%(col2)s" lw 1.3 dt (2,1)
1159 set style line 42 lt 8 lc rgb "%(col2)s" lw 1.3 dt (4,3)
1160
1161-set style line 3 lt 1 lc rgb "%(col3)s" lw 1.3
1162+set style line 3 lt 1 lc rgb "green" lw 1.7
1163 set style line 103 lt 1 lc rgb "%(col3)s" lw 1.3 dt (6,3)
1164 set style line 13 lt 2 lc rgb "%(col3)s" lw 1.3 dt (6,3)
1165 set style line 23 lt 4 lc rgb "%(col3)s" lw 1.3 dt (3,2)
1166 set style line 33 lt 6 lc rgb "%(col3)s" lw 1.3 dt (2,1)
1167 set style line 43 lt 8 lc rgb "%(col3)s" lw 1.3 dt (4,3)
1168
1169-set style line 4 lt 1 lc rgb "%(col4)s" lw 1.3
1170+set style line 4 lt 1 lc rgb "purple" lw 1.3
1171 set style line 104 lt 1 lc rgb "%(col4)s" lw 1.3 dt (6,3)
1172 set style line 14 lt 2 lc rgb "%(col4)s" lw 1.3 dt (6,3)
1173 set style line 24 lt 4 lc rgb "%(col4)s" lw 1.3 dt (3,2)
1174@@ -2383,7 +2580,6 @@
1175 set style line 38 lt 6 lc rgb "%(col8)s" lw 1.3 dt (2,1)
1176 set style line 48 lt 8 lc rgb "%(col8)s" lw 1.3 dt (4,3)
1177
1178-
1179 set style line 999 lt 1 lc rgb "gray" lw 1.3
1180
1181 safe(x,y,a) = (y == 0.0 ? a : x/y)
1182@@ -2446,8 +2642,8 @@
1183 "now be rendered by invoking gnuplot.")
1184
1185 def output_group(self, HwU_out, gnuplot_out, block_position, HwU_name,
1186- number_of_ratios = -1,
1187- uncertainties = ['scale','pdf','statitistical','merging_scale','alpsfact'],
1188+ number_of_ratios = -1,
1189+ uncertainties = ['scale','pdf','statitistical','merging_scale','alpsfact','rpa','rap'],############S.A
1190 use_band = None,
1191 ratio_correlations = True,
1192 jet_samples_to_keep=None,
1193@@ -2639,11 +2835,28 @@
1194 (mu_var_pos,mu) = self[0].set_uncertainty(type='all_scale')
1195 else:
1196 (mu_var_pos,mu) = (None,[None])
1197-
1198+
1199+#########S.A
1200+
1201+
1202 if 'pdf' in uncertainties:
1203 (PDF_var_pos,pdf) = self[0].set_uncertainty(type='PDF',lhapdfconfig=lhapdfconfig)
1204 else:
1205 (PDF_var_pos,pdf) = (None,[None])
1206+
1207+#########S.A
1208+
1209+ if 'rpa' in uncertainties:
1210+ (RPA_var_pos,rpa) = self[0].set_uncertainty(type='RPA',lhapdfconfig=lhapdfconfig)
1211+ else:
1212+ (RPA_var_pos,rpa) = (None,[None])
1213+
1214+ if 'rap' in uncertainties:
1215+ (RAP_var_pos,rap) = self[0].set_uncertainty(type='RAP',lhapdfconfig=lhapdfconfig)
1216+ else:
1217+ (RAP_var_pos,rap) = (None,[None])
1218+
1219+#########S.A
1220
1221 if 'merging_scale' in uncertainties:
1222 (merging_var_pos,merging) = self[0].set_uncertainty(type='merging')
1223@@ -2654,9 +2867,25 @@
1224 else:
1225 (alpsfact_var_pos,alpsfact) = (None,[None])
1226
1227+
1228+
1229+#########S.A
1230+
1231+
1232 uncertainties_present = list(uncertainties)
1233 if PDF_var_pos is None and 'pdf' in uncertainties_present:
1234 uncertainties_present.remove('pdf')
1235+
1236+#########S.A
1237+ uncertainties_present = list(uncertainties)
1238+ if RPA_var_pos is None and 'rpa' in uncertainties_present:
1239+ uncertainties_present.remove('rpa')
1240+
1241+ uncertainties_present = list(uncertainties)
1242+ if RAP_var_pos is None and 'rap' in uncertainties_present:
1243+ uncertainties_present.remove('rap')
1244+#########S.A
1245+
1246 if mu_var_pos is None and 'scale' in uncertainties_present:
1247 uncertainties_present.remove('scale')
1248 if merging_var_pos is None and 'merging' in uncertainties_present:
1249@@ -2688,12 +2917,34 @@
1250 raise MadGraph5Error('Not all histograms in this group specify'+\
1251 ' scale uncertainties. It is required to be able to output them'+\
1252 ' together.')
1253+
1254+###########S.A
1255+
1256 if (not PDF_var_pos is None) and\
1257 PDF_var_pos != histo.set_uncertainty(type='PDF',\
1258 lhapdfconfig=lhapdfconfig)[0]:
1259 raise MadGraph5Error('Not all histograms in this group specify'+\
1260 ' PDF uncertainties. It is required to be able to output them'+\
1261 ' together.')
1262+
1263+
1264+###########S.A
1265+ if (not RPA_var_pos is None) and\
1266+ RPA_var_pos != histo.set_uncertainty(type='RPA',\
1267+ lhapdfconfig=lhapdfconfig)[0]:
1268+ raise MadGraph5Error('Not all histograms in this group specify'+\
1269+ ' PDF uncertainties. It is required to be able to output them'+\
1270+ ' together.')
1271+
1272+
1273+ if (not RAP_var_pos is None) and\
1274+ RAP_var_pos != histo.set_uncertainty(type='RAP',\
1275+ lhapdfconfig=lhapdfconfig)[0]:
1276+ raise MadGraph5Error('Not all histograms in this group specify'+\
1277+ ' PDF uncertainties. It is required to be able to output them'+\
1278+ ' together.')
1279+###########S.A
1280+
1281 if (not merging_var_pos is None) and\
1282 merging_var_pos != histo.set_uncertainty(type='merging')[0]:
1283 raise MadGraph5Error('Not all histograms in this group specify'+\
1284@@ -2748,6 +2999,25 @@
1285 %(set_ylabel)s
1286 %(set_histo_label)s
1287 plot \\"""
1288+
1289+
1290+
1291+ subhistogram_header1 = \
1292+"""#-- rendering subhistograms '%(subhistogram_type)s'
1293+%(unset label)s
1294+%(set_format_y)s
1295+unset yrange
1296+set offsets 0 ,0 , 0.1, 0.1
1297+set origin %(origin_x).4e, %(origin_y).4e
1298+set size %(size_x).4e, %(size_y).4e
1299+set mytics %(mytics)d
1300+%(set_ytics)s
1301+%(set_format_x)s
1302+%(set_yscale)s
1303+%(set_ylabel)s
1304+%(set_histo_label)s
1305+plot \\"""
1306+
1307 replacement_dic = {}
1308
1309 replacement_dic['title'] = self[0].get_HwU_histogram_name(format='human-no_type')
1310@@ -2763,7 +3033,24 @@
1311 for PDF_var in PDF_var_pos:
1312 wgts_to_consider.append(self[0].bins.weight_labels[PDF_var])
1313 wgts_to_consider.append(self[0].bins.weight_labels[PDF_var+1])
1314- wgts_to_consider.append(self[0].bins.weight_labels[PDF_var+2])
1315+ wgts_to_consider.append(self[0].bins.weight_labels[PDF_var+2])
1316+
1317+#############S.A
1318+ if not RPA_var_pos is None:
1319+ for RPA_var in RPA_var_pos:
1320+ wgts_to_consider.append(self[0].bins.weight_labels[RPA_var])
1321+ wgts_to_consider.append(self[0].bins.weight_labels[RPA_var+1])
1322+ wgts_to_consider.append(self[0].bins.weight_labels[RPA_var+2])
1323+
1324+ if not RAP_var_pos is None:
1325+ for RAP_var in RAP_var_pos:
1326+ wgts_to_consider.append(self[0].bins.weight_labels[RAP_var])
1327+ wgts_to_consider.append(self[0].bins.weight_labels[RAP_var+1])
1328+ wgts_to_consider.append(self[0].bins.weight_labels[RAP_var+2])
1329+
1330+#############S.A
1331+
1332+
1333 if not merging_var_pos is None:
1334 for merging_var in merging_var_pos:
1335 wgts_to_consider.append(self[0].bins.weight_labels[merging_var])
1336@@ -2889,7 +3176,7 @@
1337 '%s, scale variation'%title, band='scale' in use_band)
1338 else:
1339 uncertainty_plot_lines[-1]['scale'] = \
1340- ["sqrt(-1) ls %d title '%s'"%(color_index+10,'%s, scale variation'%title)]
1341+ ["sqrt(-1) ls %d title '%s'"%(color_index+10,'%s, scale variation'%title)]
1342 # And now PDF_variation if available
1343 if not PDF_var_pos is None and len(PDF_var_pos)>0:
1344 if 'pdf' in use_band:
1345@@ -2899,6 +3186,34 @@
1346 else:
1347 uncertainty_plot_lines[-1]['pdf'] = \
1348 ["sqrt(-1) ls %d title '%s'"%(color_index+20,'%s, PDF variation'%title)]
1349+
1350+
1351+##############S.A
1352+
1353+
1354+ if not RPA_var_pos is None and len(RPA_var_pos)>0:
1355+ if 'pdf' in use_band:
1356+ uncertainty_plot_lines[-1]['rpa'] = get_uncertainty_lines(
1357+ HwU_name,block_position+i, RPA_var_pos[0]+4, color_index+20,
1358+ '%s, RPA variation'%title, band='rpa' in use_band)
1359+ else:
1360+ uncertainty_plot_lines[-1]['rpa'] = \
1361+ ["sqrt(-1) ls %d title '%s'"%(color_index+20,'%s, RPA variation'%title)]
1362+
1363+
1364+ if not RAP_var_pos is None and len(RAP_var_pos)>0:
1365+ if 'pdf' in use_band:
1366+ uncertainty_plot_lines[-1]['rap'] = get_uncertainty_lines(
1367+ HwU_name,block_position+i, RAP_var_pos[0]+4, color_index+20,
1368+ '%s, RAP variation'%title, band='rap' in use_band)
1369+ else:
1370+ uncertainty_plot_lines[-1]['rap'] = \
1371+ ["sqrt(-1) ls %d title '%s'"%(color_index+20,'%s, RAP variation'%title)]
1372+
1373+###############S.A
1374+
1375+
1376+
1377 # And now merging variation if available
1378 if not merging_var_pos is None and len(merging_var_pos)>0:
1379 if 'merging_scale' in use_band:
1380@@ -2950,6 +3265,34 @@
1381 %(HwU_name,block_position+i,PDF_var+3,color_index,\
1382 '%s PDF=%s' % (title,pdf[j].replace('_','\_'))))
1383
1384+
1385+##############S.A
1386+ if not RPA_var_pos is None:
1387+ for j,RPA_var in enumerate(RPA_var_pos):
1388+ if j!=0:
1389+ n=n+1
1390+ color_index = n%self.number_line_colors_defined+1
1391+ plot_lines.append(
1392+"'%s' index %d using (($1+$2)/2):%d ls %d title '%s'"\
1393+%(HwU_name,block_position+i,RPA_var+3,color_index,\
1394+'%s RPA=%s' % (title,rpa[j].replace('_','\_'))))
1395+
1396+
1397+
1398+ if not RAP_var_pos is None:
1399+ for j,RAP_var in enumerate(RAP_var_pos):
1400+ if j!=0:
1401+ n=n+1
1402+ color_index = n%self.number_line_colors_defined+1
1403+ plot_lines.append(
1404+"'%s' index %d using (($1+$2)/2):%d ls %d title '%s'"\
1405+%(HwU_name,block_position+i,RAP_var+3,color_index,\
1406+'%s RPA=%s' % (title,rap[j].replace('_','\_'))))
1407+
1408+##############S.A
1409+
1410+
1411+
1412 # Now add the uncertainty lines, those not using a band so that they
1413 # are not covered by those using a band after we reverse plo_lines
1414 for one_plot in uncertainty_plot_lines:
1415@@ -3084,6 +3427,51 @@
1416 uncertainty_plot_lines[-1]['pdf'] = get_uncertainty_lines(
1417 HwU_name, block_position+i, PDF_var+4, color_index+20,'',
1418 ratio=True, band='pdf' in use_band)
1419+
1420+
1421+##############S.A
1422+
1423+ if not RPA_var_pos is None:
1424+ for j,RPA_var in enumerate(RPA_var_pos):
1425+ uncertainty_plot_lines.append({})
1426+ if j==0:
1427+ color_index = k%self.number_line_colors_defined+1
1428+ else:
1429+ n=n+1
1430+ color_index = n%self.number_line_colors_defined+1
1431+ # Add the central line only if advanced pdf variation
1432+ if j>0 or rpa[j]!='none':
1433+ plot_lines.append(
1434+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)-1.0) ls %d title ''"\
1435+ %(HwU_name,block_position+i,RPA_var+3,color_index))
1436+ uncertainty_plot_lines[-1]['rpa'] = get_uncertainty_lines(
1437+ HwU_name, block_position+i, RPA_var+4, color_index+20,'',
1438+ ratio=True, band='rpa' in use_band)
1439+
1440+
1441+
1442+ if not RAP_var_pos is None:
1443+ for j,RAP_var in enumerate(RAP_var_pos):
1444+ uncertainty_plot_lines.append({})
1445+ if j==0:
1446+ color_index = k%self.number_line_colors_defined+1
1447+ else:
1448+ n=n+1
1449+ color_index = n%self.number_line_colors_defined+1
1450+ # Add the central line only if advanced pdf variation
1451+ if j>0 or rap[j]!='none':
1452+ plot_lines.append(
1453+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)-1.0) ls %d title ''"\
1454+ %(HwU_name,block_position+i,RAP_var+3,color_index))
1455+ uncertainty_plot_lines[-1]['rap'] = get_uncertainty_lines(
1456+ HwU_name, block_position+i, RAP_var+4, color_index+20,'',
1457+ ratio=True, band='rap' in use_band)
1458+
1459+
1460+##############S.A
1461+
1462+
1463+
1464 if not merging_var_pos is None:
1465 for j,merging_var in enumerate(merging_var_pos):
1466 uncertainty_plot_lines.append({})
1467@@ -3198,6 +3586,18 @@
1468 if not PDF_var_pos is None:
1469 for j,PDF_var in enumerate(PDF_var_pos):
1470 if j!=0: n=n+1
1471+
1472+###########S.A
1473+ if not RPA_var_pos is None:
1474+ for j,RPA_var in enumerate(RPA_var_pos):
1475+ if j!=0: n=n+1
1476+
1477+
1478+ if not RAP_var_pos is None:
1479+ for j,RAP_var in enumerate(RAP_var_pos):
1480+ if j!=0: n=n+1
1481+###########S.A
1482+
1483 if not merging_var_pos is None:
1484 for j,merging_var in enumerate(merging_var_pos):
1485 if j!=0: n=n+1
1486@@ -3252,6 +3652,49 @@
1487 uncertainty_plot_lines[-1]['pdf'] = get_uncertainty_lines(
1488 HwU_name, block_ratio_pos, PDF_var+4, color_index+20,'',
1489 band='pdf' in use_band)
1490+
1491+###############S.A
1492+
1493+ if not RPA_var_pos is None:
1494+ for j,RPA_var in enumerate(RPA_var_pos):
1495+ uncertainty_plot_lines.append({})
1496+ if j==0:
1497+ color_index = k%self.number_line_colors_defined+1
1498+ else:
1499+ n=n+1
1500+ color_index = n%self.number_line_colors_defined+1
1501+ # Only print out the additional central value for advanced pdf variation
1502+ if j>0 or pdf[j]!='none':
1503+ plot_lines.append(
1504+ "'%s' index %d using (($1+$2)/2):%d ls %d title ''"\
1505+ %(HwU_name,block_ratio_pos,RPA_var+3,color_index))
1506+ uncertainty_plot_lines[-1]['rpa'] = get_uncertainty_lines(
1507+ HwU_name, block_ratio_pos, RPA_var+4, color_index+20,'',
1508+ band='rpa' in use_band)
1509+
1510+
1511+
1512+
1513+ if not RAP_var_pos is None:
1514+ for j,RAP_var in enumerate(RAP_var_pos):
1515+ uncertainty_plot_lines.append({})
1516+ if j==0:
1517+ color_index = k%self.number_line_colors_defined+1
1518+ else:
1519+ n=n+1
1520+ color_index = n%self.number_line_colors_defined+1
1521+ # Only print out the additional central value for advanced pdf variation
1522+ if j>0 or pdf[j]!='none':
1523+ plot_lines.append(
1524+ "'%s' index %d using (($1+$2)/2):%d ls %d title ''"\
1525+ %(HwU_name,block_ratio_pos,RAP_var+3,color_index))
1526+ uncertainty_plot_lines[-1]['rap'] = get_uncertainty_lines(
1527+ HwU_name, block_ratio_pos, RAP_var+4, color_index+20,'',
1528+ band='rap' in use_band)
1529+
1530+
1531+###############S.A
1532+
1533 if not merging_var_pos is None:
1534 for j,merging_var in enumerate(merging_var_pos):
1535 uncertainty_plot_lines.append({})
1536@@ -3300,15 +3743,494 @@
1537 # Reverse so that bands appear first
1538 plot_lines.reverse()
1539 # Add the plot lines
1540+ if not no_uncertainties:
1541+ gnuplot_out.append(',\\\n'.join(plot_lines))
1542+ # We finish here when no ratio plot are asked for.
1543+
1544+
1545+
1546+ Reg=0
1547+ path_to_card = os.path.abspath('Cards/run_card.dat')
1548+ sub_string = 'True = rpa_choice'
1549+ with open(str(path_to_card)) as file:
1550+ lines = file.readlines()
1551+ for line in lines:
1552+ if sub_string in line:
1553+ Reg=Reg+1
1554+ else:
1555+ Reg=Reg+0
1556+
1557+
1558+ if (Reg!=1):
1559+ # Now add the tail for this group
1560+ gnuplot_out.extend(['','unset label','unset multiplot','unset xlabel','set key spacing 1','',
1561+'################################################################################'])
1562+
1563+ # Return the starting data_block position for the next histogram group
1564+ return block_position+len(self)
1565+ else:
1566+ if len(self)-n_histograms==0:
1567+ # Now add the tail for this group
1568+ gnuplot_out.extend(['','unset label','unset multiplot','',
1569+'################################################################################'])
1570+ # Return the starting data_block position for the next histogram group
1571+ return block_position+len(self)
1572+################################################################################################################
1573+################################################################################################################
1574+################################################################################################################
1575+################################################################################################################
1576+ # We can finally add the last subhistograms for the ratios.
1577+ for i, histo in enumerate(self[:n_histograms]):
1578+ if i==0: continue
1579+
1580+ PDFs=[]
1581+ PDFL=[]
1582+ PDFs=['','']
1583+ PDFL=['','']
1584+ for i in range(2):
1585+ if pdf[i].find("N1")!=-1 or pdf[i].find("1_1")!=-1 or pdf[i].find("proton")!=-1:
1586+ PDFs[i]='p'
1587+ PDFL[i]='p'
1588+ elif pdf[i].find("D2")!=-1 or pdf[i].find("2_1")!=-1:
1589+ PDFs[i]='D'
1590+ PDFL[i]='A'
1591+ elif pdf[i].find("C12")!=-1 or pdf[i].find("12_6")!=-1:
1592+ PDFs[i]='C'
1593+ PDFL[i]='A'
1594+ elif pdf[i].find("O16")!=-1 or pdf[i].find("16_8")!=-1:
1595+ PDFs[i]='O'
1596+ PDFL[i]='A'
1597+ elif pdf[i].find("Xe131")!=-1 or pdf[i].find("131_54")!=-1:
1598+ PDFs[i]='Xe'
1599+ PDFL[i]='A'
1600+ elif pdf[i].find("Au197")!=-1 or pdf[i].find("197_79")!=-1:
1601+ PDFs[i]='Au'
1602+ PDFL[i]='A'
1603+ elif pdf[i].find("Pb208")!=-1 or pdf[i].find("208_82")!=-1:
1604+ PDFs[i]='Pb'
1605+ PDFL[i]='A'
1606+ else:
1607+ PDFs[i]='X'
1608+ PDFL[i]='A'
1609+
1610+
1611+
1612+ import numpy as np
1613+ sub_string = 'ebeam1'
1614+ with open(str(path_to_card)) as file:
1615+ lines = file.readlines()
1616+ for line in lines:
1617+ if sub_string in line:
1618+ s=line
1619+ word_list = s.split()
1620+ Energy1=float(word_list[0])
1621+
1622+ sub_string = 'ebeam2'
1623+ with open(str(path_to_card)) as file:
1624+ lines = file.readlines()
1625+ for line in lines:
1626+ if sub_string in line:
1627+ s=line
1628+ word_list = s.split()
1629+ Energy2=float(word_list[0])
1630+
1631+
1632+ if (Energy1<=0.938 or Energy2<=0.938):
1633+ sqrtS=float(np.sqrt(2*Energy2*Energy1)/1000)
1634+ else:
1635+ sqrtS=float(np.sqrt(4*Energy2*Energy1)/1000)
1636+
1637+
1638+
1639+
1640+ sqrtS=float(np.sqrt(4*Energy2*Energy1)/1000)
1641+ # Add a margin on upper and lower bound.
1642+ #ymax = ymax# + 0.1# * (ymax - ymin)2
1643+ #ymin = ymin# - 0.2# * (ymax - ymin)0.5
1644+ replacement_dic['unset label'] = 'unset label'
1645+
1646+ gnuplot_out.extend(['unset multiplot'])
1647+ gnuplot_out.extend(['set multiplot'])
1648+ gnuplot_out.extend(['set key spacing 0'])
1649+ replacement_dic['origin_x'] = 0.0000e+00
1650+ replacement_dic['origin_y'] = 0.85
1651+ replacement_dic['size_y'] = 1.75000e-01
1652+ replacement_dic['size_x'] = 1.0000e+00
1653+ replacement_dic['mytics'] = 5
1654+ replacement_dic['set_ytics'] = 'set ytics auto'
1655+ replacement_dic['set_format_x'] = "set format x"
1656+ replacement_dic['set_yscale'] = "unset logscale y"
1657+ replacement_dic['set_format_y'] = 'unset format'
1658+ replacement_dic['set_ylabel'] = 'set ylabel "R_{%s%s}"'%(PDFL[0],PDFL[1])
1659+
1660+
1661+ s=self[0].get_HwU_histogram_name(format='human-no_type')
1662+ ylab='y'
1663+ ptlab='p_{T}, [GeV/c^{2}]'
1664+
1665+ if s.find("rap")!=-1 or s.find("y")!=-1:
1666+ gnuplot_out.extend(['set xlabel "%s"'%(ylab)])
1667+ if s.find("pt")!=-1:
1668+ gnuplot_out.extend(['set xlabel "%s"'%(ptlab)])
1669+
1670+ replacement_dic['set_histo_label'] = 'set label "%s + %s, {/Symbol=\\\%d}s = %.2f TeV, PDFs=%s, %s" font ",9" at graph 0, graph 1.07'%(PDFs[0],PDFs[1],326,sqrtS, pdf[0].replace('_','\\\_'), pdf[1].replace('_','\\\\_'))
1671+
1672+ replacement_dic['subhistogram_type'] = 'R_{%s%s}'%(PDFL[0],PDFL[1])
1673+ gnuplot_out.append(subhistogram_header1%replacement_dic)
1674+
1675+
1676+ uncertainty_plot_lines = []
1677+ plot_lines = []
1678+
1679+ # Some crap to get the colors right I suppose...
1680+ n=-1
1681+ n=n+1
1682+ copy_swap_re = r"perl -pe 's/^\s*(?<x1>[\+|-]?\d+(\.\d*)?([EeDd][\+|-]?\d+)?)\s*(?<x2>[\+|-]?\d+(\.\d*)?([EeDd][\+|-]?\d+)?)(?<rest>.*)\n/ $+{x1} $+{x2} $+{rest}\n$+{x2} $+{x1} $+{rest}\n/g'"
1683+ # Gnuplot escapes the antislash, so we must esacape then once more O_o.
1684+ # Gnuplot doesn't have raw strings, what a shame...
1685+ copy_swap_re = copy_swap_re.replace('\\','\\\\')
1686+
1687+ for i_histo_ratio, histo_ration in enumerate(self[n_histograms:]):
1688+ n=n+1
1689+ k=n
1690+ block_ratio_pos = block_position+n_histograms+i_histo_ratio
1691+ color_index = n%self.number_line_colors_defined+1
1692+ # Now add the subhistograms
1693+###############A
1694+
1695+ if not RPA_var_pos is None:
1696+ for j,RPA_var in enumerate(RPA_var_pos):
1697+ uncertainty_plot_lines.append({})
1698+ if j==0:
1699+ color_index = k%self.number_line_colors_defined+1
1700+ else:
1701+ n=n+1
1702+ color_index = n%self.number_line_colors_defined+1
1703+ # Add the central line only if advanced pdf variation
1704+ if j>0 or rpa[j]!='none':
1705+ plot_lines.append(
1706+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title ''"\
1707+ %(HwU_name,block_position+1,RPA_var+4,2))
1708+ plot_lines.append(
1709+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title ''"\
1710+ %(HwU_name,block_position+1,RPA_var+5,2))
1711+
1712+ plot_lines.append(
1713+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title 'R_{%s%s} central value + PDF uncertainties, LO'"\
1714+ %(HwU_name,block_position+1,RPA_var+3,1,PDFL[0],PDFL[1]))
1715+ plot_lines.append(
1716+' "<%s %s" index %d using ($1):(safe($%d,$3,1.0)):(safe($%d,$3,1.0)) with filledcurve ls %d fs transparent pattern 2 title " " '\
1717+ %(copy_swap_re,HwU_name,block_position+1,RPA_var+4,RPA_var+5,2))
1718+###############ANton
1719+ plot_lines.append("1 ls 999 title '', 0.8 ls 999 title '', 1.2 ls 999 title ''")
1720+
1721+ # Now add the uncertainty lines, those not using a band so that they
1722+ # are not covered by those using a band after we reverse plo_lines
1723+ for one_plot in uncertainty_plot_lines:
1724+ for uncertainty_type, lines in one_plot.items():
1725+ if not uncertainty_type in use_band:
1726+ plot_lines.extend(lines)
1727+ # then those using a band
1728+ for one_plot in uncertainty_plot_lines:
1729+ for uncertainty_type, lines in one_plot.items():
1730+ if uncertainty_type in use_band:
1731+ plot_lines.extend(lines)
1732+
1733+ # Reverse so that bands appear first
1734+ plot_lines.reverse()
1735+ # Add the plot lines
1736+ if not no_uncertainties:
1737+ gnuplot_out.append(',\\\n'.join(plot_lines))
1738+ # We finish here when no ratio plot are asked for.
1739+ if len(self)-n_histograms==0:
1740+ # Now add the tail for this group
1741+ gnuplot_out.extend(['','unset label','',
1742+'################################################################################'])
1743+ # Return the starting data_block position for the next histogram group
1744+ return block_position+len(self)
1745+####################################################################################################################22222
1746+#################################################################################################################22222
1747+ # We can finally add the last subhistograms for the ratios.
1748+ for i, histo in enumerate(self[:n_histograms]):
1749+ if i==0: continue
1750+
1751+ replacement_dic['unset label'] = 'unset label'
1752+
1753+ replacement_dic['origin_x'] = 0.0000e+00
1754+ replacement_dic['origin_y'] = 0.60
1755+ replacement_dic['size_y'] = 1.75000e-01
1756+ replacement_dic['size_x'] = 1.0000e+00
1757+ replacement_dic['mytics'] = 5
1758+
1759+ replacement_dic['set_ytics'] = 'set ytics auto'
1760+ replacement_dic['set_format_x'] = "set format x"
1761+ replacement_dic['set_yscale'] = "unset logscale y"
1762+ replacement_dic['set_format_y'] = 'unset format'
1763+ replacement_dic['set_ylabel'] = 'set ylabel "R_{%s%s}"'%(PDFL[0],PDFL[1])
1764+
1765+ replacement_dic['set_histo_label'] = \
1766+ 'set label "%s + %s, {/Symbol=\\\%d}s = %.2f TeV, PDFs=%s, %s" font ",9" at graph 0, graph 1.07'%(PDFs[0],PDFs[1],326,sqrtS, pdf[0].replace('_','\\\_'), pdf[1].replace('_','\\\\_'))#replace('_','\_')
1767+
1768+ replacement_dic['subhistogram_type'] = 'R_{%s%s}'%(PDFL[0],PDFL[1])
1769+ gnuplot_out.append(subhistogram_header1%replacement_dic)
1770+
1771+
1772+ uncertainty_plot_lines = []
1773+ plot_lines = []
1774+
1775+ # Some crap to get the colors right I suppose...
1776+ n=-1
1777+ n=n+1
1778+
1779+ for i_histo_ratio, histo_ration in enumerate(self[n_histograms:]):
1780+ n=n+1
1781+ k=n
1782+ block_ratio_pos = block_position+n_histograms+i_histo_ratio
1783+ color_index = n%self.number_line_colors_defined+1
1784+ # Now add the subhistograms
1785+###############A
1786+
1787+ if not RPA_var_pos is None:
1788+ for j,RPA_var in enumerate(RPA_var_pos):
1789+ uncertainty_plot_lines.append({})
1790+ if j==0:
1791+ color_index = k%self.number_line_colors_defined+1
1792+ else:
1793+ n=n+1
1794+ color_index = n%self.number_line_colors_defined+1
1795+ # Add the central line only if advanced pdf variation
1796+ if j>0 or rpa[j]!='none':
1797+ plot_lines.append(
1798+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title ''"\
1799+ %(HwU_name,block_position,RPA_var+4,2))
1800+ plot_lines.append(
1801+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title ''"\
1802+ %(HwU_name,block_position,RPA_var+5,2))
1803+
1804+ plot_lines.append(
1805+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title 'R_{%s%s} central value + PDF uncertainties, NLO'"\
1806+ %(HwU_name,block_position,RPA_var+3,1,PDFL[0],PDFL[1]))
1807+ plot_lines.append(
1808+' "<%s %s" index %d using ($1):(safe($%d,$3,1.0)):(safe($%d,$3,1.0)) with filledcurve ls %d fs transparent pattern 2 title " " '\
1809+ %(copy_swap_re,HwU_name,block_position,RPA_var+4,RPA_var+5,2))
1810+
1811+###############ANton
1812+ plot_lines.append("1 ls 999 title '', 0.8 ls 999 title '', 1.2 ls 999 title ''")
1813+
1814+ # Now add the uncertainty lines, those not using a band so that they
1815+ # are not covered by those using a band after we reverse plo_lines
1816+ for one_plot in uncertainty_plot_lines:
1817+ for uncertainty_type, lines in one_plot.items():
1818+ if not uncertainty_type in use_band:
1819+ plot_lines.extend(lines)
1820+ # then those using a band
1821+ for one_plot in uncertainty_plot_lines:
1822+ for uncertainty_type, lines in one_plot.items():
1823+ if uncertainty_type in use_band:
1824+ plot_lines.extend(lines)
1825+
1826+ # Reverse so that bands appear first
1827+ plot_lines.reverse()
1828+ # Add the plot lines
1829+ if not no_uncertainties:
1830+ gnuplot_out.append(',\\\n'.join(plot_lines))
1831+ # We finish here when no ratio plot are asked for.
1832+ if len(self)-n_histograms==0:
1833+ # Now add the tail for this group
1834+ gnuplot_out.extend(['','unset label','',
1835+'################################################################################'])
1836+ # Return the starting data_block position for the next histogram group
1837+ return block_position+len(self)
1838+####################################################################################################################22222
1839+#################################################################################################################22222
1840+ # We can finally add the last subhistograms for the ratios.
1841+ for i, histo in enumerate(self[:n_histograms]):
1842+ if i==0: continue
1843+
1844+ replacement_dic['unset label'] = 'unset label'
1845+
1846+ replacement_dic['ymin'] = ymin
1847+ replacement_dic['ymax'] = ymax
1848+ replacement_dic['origin_x'] = 0.0000e+00
1849+ replacement_dic['origin_y'] = 0.35
1850+ replacement_dic['size_y'] = 1.75000e-01
1851+ replacement_dic['size_x'] = 1.0000e+00
1852+ replacement_dic['mytics'] = 5
1853+ replacement_dic['set_ytics'] = 'set ytics auto'
1854+ replacement_dic['set_format_x'] = "set format x"
1855+ replacement_dic['set_yscale'] = "unset logscale y"
1856+ replacement_dic['set_format_y'] = 'unset format'
1857+ replacement_dic['set_ylabel'] = 'set ylabel "R_{%s%s}"'%(PDFL[1],PDFL[0])
1858+ replacement_dic['set_histo_label'] = \
1859+ 'set label "%s + %s, {/Symbol=\\\%d}s = %.2f TeV, PDFs=%s, %s" font ",9" at graph 0, graph 1.07'%(PDFs[1],PDFs[0],326,sqrtS, pdf[1].replace('_','\\\_'), pdf[0].replace('_','\\\\_'))#replace('_','\_')
1860+ replacement_dic['subhistogram_type'] = 'R_{%s%s}'%(PDFL[1],PDFL[0])
1861+ gnuplot_out.append(subhistogram_header1%replacement_dic)
1862+
1863+
1864+ uncertainty_plot_lines = []
1865+ plot_lines = []
1866+
1867+ # Some crap to get the colors right I suppose...
1868+ n=-1
1869+ n=n+1
1870+
1871+ for i_histo_ratio, histo_ration in enumerate(self[n_histograms:]):
1872+ n=n+1
1873+ k=n
1874+ block_ratio_pos = block_position+n_histograms+i_histo_ratio
1875+ color_index = n%self.number_line_colors_defined+1
1876+ # Now add the subhistograms
1877+###############ANton
1878+ if not RAP_var_pos is None:
1879+ for j,RAP_var in enumerate(RAP_var_pos):
1880+ uncertainty_plot_lines.append({})
1881+ if j==0:
1882+ color_index = k%self.number_line_colors_defined+1
1883+ else:
1884+ n=n+1
1885+ color_index = n%self.number_line_colors_defined+1
1886+ # Add the central line only if advanced pdf variation
1887+ if j>0 or rap[j]!='none':
1888+ plot_lines.append(
1889+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title ''"\
1890+ %(HwU_name,block_position+1,RAP_var+4,4))
1891+ plot_lines.append(
1892+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title ''"\
1893+ %(HwU_name,block_position+1,RAP_var+5,4))
1894+ plot_lines.append(
1895+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title 'R_{%s%s} central value + PDF uncertainties, LO'"\
1896+ %(HwU_name,block_position+1,RAP_var+3,3,PDFL[1],PDFL[0]))
1897+ plot_lines.append(
1898+' "<%s %s" index %d using ($1):(safe($%d,$3,1.0)):(safe($%d,$3,1.0)) with filledcurve ls %d fs transparent pattern 2 title " " '\
1899+ %(copy_swap_re,HwU_name,block_position+1,RAP_var+4,RAP_var+5,4))
1900+###############ANton
1901+ plot_lines.append("1 ls 999 title '', 0.8 ls 999 title '', 1.2 ls 999 title ''")
1902+
1903+ # Now add the uncertainty lines, those not using a band so that they
1904+ # are not covered by those using a band after we reverse plo_lines
1905+ for one_plot in uncertainty_plot_lines:
1906+ for uncertainty_type, lines in one_plot.items():
1907+ if not uncertainty_type in use_band:
1908+ plot_lines.extend(lines)
1909+ # then those using a band
1910+ for one_plot in uncertainty_plot_lines:
1911+ for uncertainty_type, lines in one_plot.items():
1912+ if uncertainty_type in use_band:
1913+ plot_lines.extend(lines)
1914+
1915+ # Reverse so that bands appear first
1916+ plot_lines.reverse()
1917+ # Add the plot lines
1918+ if not no_uncertainties:
1919+ gnuplot_out.append(',\\\n'.join(plot_lines))
1920+ # We finish here when no ratio plot are asked for.
1921+ if len(self)-n_histograms==0:
1922+ # Now add the tail for this group
1923+ gnuplot_out.extend(['','unset label','',
1924+'################################################################################'])
1925+ # Return the starting data_block position for the next histogram group
1926+ return block_position+len(self)
1927+####################################################################################################################22222
1928+####################################################################################################################22222
1929+
1930+ # We can finally add the last subhistograms for the ratios.
1931+ for i, histo in enumerate(self[:n_histograms]):
1932+ if i==0: continue
1933+
1934+ replacement_dic['unset label'] = 'unset label'
1935+ replacement_dic['ymin'] = ymin
1936+ replacement_dic['ymax'] = ymax
1937+ replacement_dic['origin_x'] = 0.0000e+00
1938+ replacement_dic['origin_y'] = 0.1
1939+ replacement_dic['size_y'] = 1.75000e-01
1940+ replacement_dic['size_x'] = 1.0000e+00
1941+ replacement_dic['mytics'] = 5
1942+ replacement_dic['set_ytics'] = 'set ytics auto'
1943+ replacement_dic['set_format_x'] = "set format x"
1944+ replacement_dic['set_yscale'] = "unset logscale y"
1945+ replacement_dic['set_format_y'] = 'unset format'
1946+ replacement_dic['set_ylabel'] = 'set ylabel "R_{%s%s}"'%(PDFL[1],PDFL[0])
1947+ replacement_dic['subhistogram_type'] = 'R_{%s%s}'%(PDFL[1],PDFL[0])
1948+
1949+ replacement_dic['set_histo_label'] = \
1950+ 'set label "%s + %s, {/Symbol=\\\%d}s = %.2f TeV, PDFs=%s, %s" font ",9" at graph 0, graph 1.07'%(PDFs[1],PDFs[0],326,sqrtS, pdf[1].replace('_','\\\_'), pdf[0].replace('_','\\\\_'))#replace('_','\_')
1951+
1952+ gnuplot_out.append(subhistogram_header1%replacement_dic)
1953+
1954+ uncertainty_plot_lines = []
1955+ plot_lines = []
1956+
1957+ # Some crap to get the colors right I suppose...
1958+ n=-1
1959+ n=n+1
1960+
1961+ for i_histo_ratio, histo_ration in enumerate(self[n_histograms:]):
1962+ n=n+1
1963+ k=n
1964+ block_ratio_pos = block_position+n_histograms+i_histo_ratio
1965+ color_index = n%self.number_line_colors_defined+1
1966+ # Now add the subhistograms
1967+###############A
1968+ if not RAP_var_pos is None:
1969+ for j,RAP_var in enumerate(RAP_var_pos):
1970+ uncertainty_plot_lines.append({})
1971+ if j==0:
1972+ color_index = k%self.number_line_colors_defined+1
1973+ else:
1974+ n=n+1
1975+ color_index = n%self.number_line_colors_defined+1
1976+ # Add the central line only if advanced pdf variation
1977+ if j>0 or rap[j]!='none':
1978+ plot_lines.append(
1979+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title ''"\
1980+ %(HwU_name,block_position,RAP_var+4,4))
1981+ plot_lines.append(
1982+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title ''"\
1983+ %(HwU_name,block_position,RAP_var+5,4))
1984+ #uncertainty_plot_lines[-1]['rap'] = get_uncertainty_lines(
1985+ #HwU_name, block_position+i, RAP_var+4, color_index+20,'',
1986+ #ratio=True, band='rap' in use_band)
1987+ plot_lines.append(
1988+"'%s' index %d using (($1+$2)/2):(safe($%d,$3,1.0)) ls %d title 'R_{%s%s} central value + PDF uncertainties, NLO'"\
1989+ %(HwU_name,block_position,RAP_var+3,3,PDFL[1],PDFL[0]))
1990+ plot_lines.append(
1991+' "<%s %s" index %d using ($1):(safe($%d,$3,1.0)):(safe($%d,$3,1.0)) with filledcurve ls %d fs transparent pattern 2 title " " '\
1992+ %(copy_swap_re,HwU_name,block_position,RAP_var+4,RAP_var+5,4))
1993+###############A
1994+ plot_lines.append("1 ls 999 title '', 0.8 ls 999 title '', 1.2 ls 999 title ''")
1995+
1996+ # Now add the uncertainty lines, those not using a band so that they
1997+ # are not covered by those using a band after we reverse plo_lines
1998+ for one_plot in uncertainty_plot_lines:
1999+ for uncertainty_type, lines in one_plot.items():
2000+ if not uncertainty_type in use_band:
2001+ plot_lines.extend(lines)
2002+ # then those using a band
2003+ for one_plot in uncertainty_plot_lines:
2004+ for uncertainty_type, lines in one_plot.items():
2005+ if uncertainty_type in use_band:
2006+ plot_lines.extend(lines)
2007+
2008+ #plot_lines.append("1.0 ls 999 title ''")
2009+
2010+ # Reverse so that bands appear first
2011+ plot_lines.reverse()
2012+ # Add the plot lines
2013 gnuplot_out.append(',\\\n'.join(plot_lines))
2014
2015 # Now add the tail for this group
2016- gnuplot_out.extend(['','unset label','',
2017+ gnuplot_out.extend(['','unset label','unset multiplot','unset xlabel','set key spacing 1','',
2018 '################################################################################'])
2019
2020 # Return the starting data_block position for the next histogram group
2021 return block_position+len(self)
2022
2023+################################################################################################################
2024+################################################################################################################
2025+################################################################################################################
2026+################################################################################################################
2027+
2028+
2029 ################################################################################
2030 ## matplotlib related function
2031 ################################################################################
2032@@ -3472,13 +4394,13 @@
2033 '--only_scale','--only_pdf','--only_stat','--only_merging','--only_alpsfact',
2034 '--variations','--band','--central_only', '--lhapdf-config','--titles',
2035 '--keep_all_weights','--colours']
2036- n_ratios = -1
2037- uncertainties = ['scale','pdf','statistical','merging_scale','alpsfact']
2038+ n_ratios = -1
2039+ uncertainties = ['scale','pdf','statistical','merging_scale','alpsfact','rpa','rap']#,'scale_pa']#################S.A
2040 # The list of type of uncertainties for which to use bands. None is a 'smart' default
2041 use_band = None
2042 auto_open = True
2043- ratio_correlations = True
2044- consider_reweights = ['pdf','scale','murmuf_scales','merging_scale','alpsfact']
2045+ ratio_correlations = True
2046+ consider_reweights = ['pdf','scale','murmuf_scales','merging_scale','alpsfact','rpa','rap']#,'scale_pa']#################S.A
2047
2048 def log(msg):
2049 print("histograms.py :: %s"%str(msg))

Subscribers

People subscribed via source and target branches

to all changes: