← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to Template/NLO/SubProcesses/reweight_xsec_events.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/reweight_xsec_events.f
 
1
      program reweight_xsec_events
 
2
c Given a LH file that contains an <rwgt> part, computes the scale 
 
3
c and/or PDF dependence through reweighting. A new file is created,
 
4
c which does not contain the <rwgt> part, but retains only the 
 
5
c information on the maximum and minimum weights due to scale
 
6
c and PDF variations
 
7
c Compile with makefile_rwgt
 
8
      implicit none
 
9
      include "nexternal.inc"
 
10
      include "genps.inc"
 
11
      include "nFKSconfigs.inc"
 
12
      include "reweight_all.inc"
 
13
      include "run.inc"
 
14
      character*7 pdlabel,epa_label
 
15
      integer lhaid
 
16
      common/to_pdf/lhaid,pdlabel,epa_label
 
17
      integer maxevt,ifile,ofile,i,jj,isave,ii
 
18
      double precision saved_weight
 
19
      logical unweighted
 
20
      integer IDBMUP(2),PDFGUP(2),PDFSUP(2),IDWTUP,NPRUP,LPRUP
 
21
      double precision EBMUP(2),XSECUP,XERRUP,XMAXUP
 
22
      INTEGER MAXNUP
 
23
      PARAMETER (MAXNUP=500)
 
24
      INTEGER NUP,IDPRUP,IDUP(MAXNUP),ISTUP(MAXNUP),
 
25
     # MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP)
 
26
      DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
27
     # PUP(5,MAXNUP),VTIMUP(MAXNUP),SPINUP(MAXNUP)
 
28
      integer isorh_lhe,ifks_lhe,jfks_lhe,fksfather_lhe,ipartner_lhe
 
29
      double precision scale1_lhe,scale2_lhe,percentage
 
30
      integer kwgtinfo,kexternal,jwgtnumpartn
 
31
      double precision wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
32
      double precision xmuR_over_ref,xmuF1_over_ref,xmuF2_over_ref,
 
33
     # xQES_over_ref,pr_muR_over_ref,pr_muF1_over_ref,pr_muF2_over_ref,
 
34
     # tmp,yfactR(maxscales),yfactF(maxscales),xsecPDFr(0:maxPDFs)
 
35
      double precision xsecPDFr_acc(0:maxPDFs),xsecScale_acc(maxscales
 
36
     $     ,maxscales)
 
37
      double precision compute_rwgt_wgt_Sev,compute_rwgt_wgt_Sev_nbody
 
38
     &     ,compute_rwgt_wgt_Hev
 
39
      integer kr,kf,n,nng,nps,npairs,nsets,izero,itmp,idpdf(0:maxPDFs)
 
40
      parameter (izero=0)
 
41
      integer lef
 
42
      character*80 event_file,fname1
 
43
      character*140 buff
 
44
      character*10 MonteCarlo
 
45
      character*9 ch1
 
46
      character*20 parm(20)
 
47
      double precision value(20)
 
48
      logical AddInfoLHE
 
49
      external compute_rwgt_wgt_Sev,compute_rwgt_wgt_Sev_nbody
 
50
     &     ,compute_rwgt_wgt_Hev
 
51
      integer i_process
 
52
      common/c_addwrite/i_process
 
53
c
 
54
      call setrun                !Sets up run parameters
 
55
 
 
56
 
 
57
      write(*,*) 'Enter event file name'
 
58
      read(*,*) event_file
 
59
 
 
60
      write(*,*)'Enter 1 to save all cross sections on tape'
 
61
      write(*,*)'      0 otherwise'
 
62
      read(*,*)isave
 
63
      if(isave.eq.1)then
 
64
        isave=9
 
65
      else
 
66
        isave=0
 
67
      endif
 
68
 
 
69
 
 
70
      xQES_over_ref=QES_over_ref
 
71
      xmuR_over_ref=muR_over_ref
 
72
      xmuF1_over_ref=muF1_over_ref
 
73
      xmuF2_over_ref=muF2_over_ref
 
74
      write(*,*) 'Using:'  
 
75
      write(*,*) 'QES_over_ref: ', xQES_over_ref
 
76
      write(*,*) 'muR_over_ref: ', xmuR_over_ref
 
77
      write(*,*) 'muF1_over_ref: ', xmuF1_over_ref
 
78
      write(*,*) 'muF2_over_ref: ', xmuF2_over_ref
 
79
      if (xmuF1_over_ref .ne. xmuF2_over_ref) then
 
80
          write(*,*) "The variables muF1_over_ref and muF2_over_ref" //
 
81
     1     " have to be set equal in the run_card.dat." //
 
82
     1     " Run cannot continue, quitting..."
 
83
          stop
 
84
      endif
 
85
 
 
86
      if(do_rwgt_scale)then
 
87
        yfactR(1)=1.d0
 
88
        yfactR(2)=rw_Rscale_up
 
89
        yfactR(3)=rw_Rscale_down
 
90
        yfactF(1)=1.d0
 
91
        yfactF(2)=rw_Fscale_up
 
92
        yfactF(3)=rw_Fscale_down
 
93
        write(*,*) 'Doing scale reweight:'
 
94
        write(*,*) rw_Fscale_down, ' < mu_F < ', rw_Fscale_up
 
95
        write(*,*) rw_Rscale_down, ' < mu_R < ', rw_Rscale_up
 
96
        numscales=3
 
97
      else
 
98
        numscales=0
 
99
      endif
 
100
 
 
101
c Note: when ipdf#0, the central PDF set will be used also as a reference
 
102
c for the scale uncertainty
 
103
      if(do_rwgt_pdf)then
 
104
 
 
105
        idpdf(0)=lhaid
 
106
        idpdf(1)=pdf_set_min
 
107
        itmp=pdf_set_max
 
108
        nsets=itmp-idpdf(1)+1
 
109
        write(*,*) 'Doing PDF reweight:'
 
110
        write(*,*) 'Central set id: ', idpdf(0)
 
111
        write(*,*) 'Min error set id: ', idpdf(1)
 
112
        write(*,*) 'Max error set id: ', itmp
 
113
        if(mod(nsets,2).ne.0)then
 
114
          write(*,*)'The number of error sets must be even',nsets
 
115
          stop
 
116
        else
 
117
          npairs=nsets/2
 
118
        endif
 
119
        do i=2,nsets
 
120
          idpdf(i)=idpdf(1)+i-1
 
121
        enddo
 
122
        if(nsets.gt.maxPDFs)then
 
123
          write(*,*)'Too many PDFs: increase maxPDFs in reweight0.inc'
 
124
          stop
 
125
        endif
 
126
c
 
127
        value(1)=idpdf(0)
 
128
        parm(1)='DEFAULT'
 
129
        call pdfset(parm,value)
 
130
c
 
131
        numPDFpairs=npairs
 
132
      else
 
133
        numPDFpairs=0
 
134
      endif
 
135
 
 
136
c$$$      call fk88strcat(event_file,'.rwgt',fname1)
 
137
      lef=index(event_file,' ')-1
 
138
      fname1=event_file(1:lef)//'.rwgt'
 
139
 
 
140
      ifile=34
 
141
      open (unit=ifile,file=event_file,status='old')
 
142
      AddInfoLHE=.true.
 
143
      unweighted=.true.
 
144
      call read_lhef_header(ifile,maxevt,MonteCarlo)
 
145
      call read_lhef_init(ifile,
 
146
     &     IDBMUP,EBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,
 
147
     &     XSECUP,XERRUP,XMAXUP,LPRUP)
 
148
 
 
149
      do i=1,min(10,maxevt)
 
150
        call read_lhef_event(ifile,
 
151
     &       NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
152
     &       IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
 
153
 
 
154
        if(buff(1:1).ne.'#')then
 
155
          write(*,*)'This event file cannot be reweighted [1]',i
 
156
          stop
 
157
        endif
 
158
        read(buff,*)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,
 
159
     #                    fksfather_lhe,ipartner_lhe,
 
160
     #                    scale1_lhe,scale2_lhe,
 
161
     #                    kwgtinfo,kexternal,jwgtnumpartn,
 
162
     #         wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
163
        if(kwgtinfo.lt.1.or.kwgtinfo.gt.5)then
 
164
          write(*,*)'This event file cannot be reweighted [2]',i
 
165
          write(*,*)kwgtinfo
 
166
          stop
 
167
        endif
 
168
        if(i.eq.1)then
 
169
          saved_weight=abs(XWGTUP)
 
170
        else
 
171
          unweighted=unweighted.and.
 
172
     #               abs(1.d0-abs(XWGTUP)/saved_weight).lt.1.d-5
 
173
        endif
 
174
      enddo
 
175
 
 
176
      write(*,*)'  '
 
177
      if(unweighted)then
 
178
        write(*,*)'The events appear to be unweighted'
 
179
        write(*,*)' Will store the ratios of recomputed weights'
 
180
        write(*,*)' over reference weights'
 
181
      else
 
182
        write(*,*)'The events appear to be weighted'
 
183
        write(*,*)' Will store recomputed weights'
 
184
      endif
 
185
 
 
186
      rewind(34)
 
187
 
 
188
      ofile=35
 
189
      open(unit=ofile,file=fname1,status='unknown')
 
190
 
 
191
      call read_lhef_header(ifile,maxevt,MonteCarlo)
 
192
      call write_lhef_header(ofile,maxevt,MonteCarlo)
 
193
      call read_lhef_init(ifile,
 
194
     &     IDBMUP,EBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,
 
195
     &     XSECUP,XERRUP,XMAXUP,LPRUP)
 
196
      call write_lhef_init(ofile,
 
197
     &     IDBMUP,EBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,
 
198
     &     XSECUP,XERRUP,XMAXUP,LPRUP)
 
199
 
 
200
      nScontributions=1
 
201
 
 
202
c To keep track of the accumulated results:
 
203
        do ii=1,numscales
 
204
           do jj=1,numscales
 
205
              xsecScale_acc(jj,ii)=0d0
 
206
           enddo
 
207
        enddo
 
208
        do n=0,nsets
 
209
           xsecPDFr_acc(n)=0d0
 
210
        enddo
 
211
 
 
212
c Determine the flavor map between the NLO and Born
 
213
      call find_iproc_map()
 
214
 
 
215
      do i=1,maxevt
 
216
        call read_lhef_event(ifile,
 
217
     &       NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
218
     &       IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
 
219
        call reweight_fill_extra_inverse()
 
220
 
 
221
        if(buff(1:1).ne.'#')then
 
222
          write(*,*)'This event file cannot be reweighted [3]',i
 
223
          stop
 
224
        endif
 
225
        read(buff,*)ch1,iSorH_lhe,ifks_lhe,jfks_lhe,
 
226
     #                    fksfather_lhe,ipartner_lhe,
 
227
     #                    scale1_lhe,scale2_lhe,
 
228
     #                    kwgtinfo,kexternal,jwgtnumpartn,
 
229
     #         wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
230
        if(kwgtinfo.lt.1.or.kwgtinfo.gt.5)then
 
231
          write(*,*)'This event file cannot be reweighted [4]',i
 
232
          write(*,*)kwgtinfo
 
233
          stop
 
234
        endif
 
235
        if(wgtcentral.ne.0.d0.or.wgtmumin.ne.0.d0.or.
 
236
     #     wgtmumax.ne.0.d0.or.wgtpdfmin.ne.0.d0.or.
 
237
     #     wgtpdfmax.ne.0.d0)then
 
238
          write(*,*)'This event file was already reweighted',i
 
239
          write(*,*)wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
240
          stop
 
241
        endif
 
242
 
 
243
        if (kwgtinfo.eq.5) call reweight_settozero()
 
244
 
 
245
        if(do_rwgt_scale)then
 
246
 
 
247
          wgtmumin=1.d40
 
248
          wgtmumax=-1.d40
 
249
 
 
250
          do kr=1,3
 
251
            do kf=1,3
 
252
              wgtref=0d0
 
253
              pr_muR_over_ref=xmuR_over_ref*yfactR(kr)
 
254
              pr_muF1_over_ref=xmuF1_over_ref*yfactF(kf)
 
255
              pr_muF2_over_ref=pr_muF1_over_ref
 
256
              wgtxsecmu(kr,kf)=0d0
 
257
              if(iSorH_lhe.eq.1)then
 
258
c The nbody contributions
 
259
                 if (kwgtinfo.eq.5) then
 
260
                    call fill_reweight0inc_nbody(i_process)
 
261
                    wgtxsecmu(kr,kf)=wgtxsecmu(kr,kf)
 
262
     &                   +compute_rwgt_wgt_Sev_nbody(pr_muR_over_ref
 
263
     &                   ,pr_muF1_over_ref, pr_muF2_over_ref
 
264
     &                   ,xQES_over_ref, kwgtinfo)
 
265
                    call reweight_settozero()
 
266
                 endif
 
267
                 do ii=1,nScontributions
 
268
                    nFKSprocess_used=nFKSprocess_reweight(ii)
 
269
                    if (kwgtinfo.eq.5)
 
270
     &                   call fill_reweight0inc(nFKSprocess_used*2-1
 
271
     $                   ,i_process)
 
272
                    wgtxsecmu(kr,kf)=wgtxsecmu(kr,kf)+
 
273
     &                   compute_rwgt_wgt_Sev(pr_muR_over_ref
 
274
     &                   ,pr_muF1_over_ref, pr_muF2_over_ref
 
275
     &                   ,xQES_over_ref, kwgtinfo)
 
276
                    if (kwgtinfo.eq.5) call reweight_settozero()
 
277
                 enddo
 
278
              elseif(iSorH_lhe.eq.2)then
 
279
                 if (kwgtinfo.eq.5)
 
280
     &                call fill_reweight0inc(nFKSprocess_used*2,
 
281
     &                i_process)
 
282
                 wgtxsecmu(kr,kf)=wgtxsecmu(kr,kf)+
 
283
     &                compute_rwgt_wgt_Hev(pr_muR_over_ref
 
284
     &                ,pr_muF1_over_ref, pr_muF2_over_ref
 
285
     &                ,xQES_over_ref, kwgtinfo)
 
286
                 if (kwgtinfo.eq.5) call reweight_settozero()
 
287
              else
 
288
                 write(*,*)'Invalid value of iSorH_lhe',iSorH_lhe
 
289
                 stop
 
290
              endif
 
291
c
 
292
              tmp=wgtxsecmu(kr,kf)
 
293
              if(tmp.lt.wgtmumin)wgtmumin=tmp
 
294
              if(tmp.gt.wgtmumax)wgtmumax=tmp
 
295
            enddo
 
296
          enddo
 
297
 
 
298
          if (kwgtinfo.eq.5) then
 
299
             if (iSorH_lhe.eq.1) then
 
300
                wgtref=wgtref_nbody_all(i_process)
 
301
                do ii=1,nScontributions
 
302
                   wgtref=wgtref+ wgtref_all(nFKSprocess_reweight(ii)*2
 
303
     $                  -1,i_process)
 
304
                enddo
 
305
             else
 
306
                wgtref=wgtref_all(nFKSprocess_used*2,i_process)
 
307
             endif
 
308
          endif
 
309
 
 
310
          if(unweighted)then
 
311
            wgtcentral=wgtxsecmu(1,1)/wgtref
 
312
            wgtmumin=wgtmumin/wgtref
 
313
            wgtmumax=wgtmumax/wgtref
 
314
          else
 
315
            wgtcentral=wgtxsecmu(1,1)
 
316
          endif
 
317
 
 
318
        endif
 
319
 
 
320
        if(do_rwgt_pdf)then
 
321
 
 
322
          do n=0,nsets
 
323
             wgtref=0d0
 
324
             call InitPDF(n)
 
325
             wgtxsecPDF(n)=0d0
 
326
 
 
327
             if(iSorH_lhe.eq.1)then
 
328
c The nbody contributions
 
329
                if (kwgtinfo.eq.5) then
 
330
                   call fill_reweight0inc_nbody(i_process)
 
331
                   wgtxsecPDF(n)=wgtxsecPDF(n)
 
332
     $                  +compute_rwgt_wgt_Sev_nbody(xmuR_over_ref
 
333
     $                  ,xmuF1_over_ref, xmuF2_over_ref ,xQES_over_ref,
 
334
     $                  kwgtinfo)
 
335
                   call reweight_settozero()
 
336
                endif
 
337
                do ii=1,nScontributions
 
338
                   nFKSprocess_used=nFKSprocess_reweight(ii)
 
339
                   if (kwgtinfo.eq.5)
 
340
     &                  call fill_reweight0inc(nFKSprocess_used*2-1
 
341
     $                   ,i_process)
 
342
                   wgtxsecPDF(n)=wgtxsecPDF(n)+
 
343
     &                  compute_rwgt_wgt_Sev(xmuR_over_ref
 
344
     &                  ,xmuF1_over_ref, xmuF2_over_ref ,xQES_over_ref,
 
345
     &                  kwgtinfo)
 
346
                   if (kwgtinfo.eq.5) call reweight_settozero()
 
347
                enddo
 
348
             elseif(iSorH_lhe.eq.2)then
 
349
                if (kwgtinfo.eq.5)
 
350
     &               call fill_reweight0inc(nFKSprocess_used*2,
 
351
     &               i_process)
 
352
                wgtxsecPDF(n)=wgtxsecPDF(n)+
 
353
     &               compute_rwgt_wgt_Hev(xmuR_over_ref ,xmuF1_over_ref,
 
354
     &               xmuF2_over_ref ,xQES_over_ref, kwgtinfo)
 
355
                if (kwgtinfo.eq.5) call reweight_settozero()
 
356
             else
 
357
                write(*,*)'Invalid value of iSorH_lhe',iSorH_lhe
 
358
                stop
 
359
             endif
 
360
c
 
361
             if (kwgtinfo.eq.5) then
 
362
                if (iSorH_lhe.eq.1) then
 
363
                   wgtref=wgtref_nbody_all(i_process)
 
364
                   do ii=1,nScontributions
 
365
                      wgtref=wgtref+wgtref_all(nFKSprocess_reweight(ii)
 
366
     $                     *2-1,i_process)
 
367
                   enddo
 
368
                else
 
369
                   wgtref=wgtref_all(nFKSprocess_used*2,i_process)
 
370
                endif
 
371
             endif
 
372
             
 
373
             if(unweighted)then
 
374
                xsecPDFr(n)=wgtxsecPDF(n)/wgtref
 
375
             else
 
376
                xsecPDFr(n)=wgtxsecPDF(n)
 
377
             endif
 
378
          enddo
 
379
 
 
380
          if(do_rwgt_scale)then
 
381
            if(abs(xsecPDFr(0)/wgtcentral-1.d0).gt.1.d-6)then
 
382
              write(*,*)'Central valued computed with mu and PDF differ'
 
383
              write(*,*)xsecPDFr(0),wgtcentral
 
384
              stop
 
385
            endif
 
386
          else
 
387
            wgtcentral=xsecPDFr(0)
 
388
c The following serves to write on tape the reference cross section
 
389
c computed with the new parameters
 
390
            wgtxsecmu(1,1)=wgtxsecPDF(0)
 
391
          endif
 
392
 
 
393
          wgtpdfmin=0.d0
 
394
          wgtpdfmax=0.d0
 
395
 
 
396
          do n=1,npairs
 
397
            nps=2*n-1
 
398
            nng=2*n
 
399
 
 
400
            wgtpdfmin=wgtpdfmin+
 
401
     #                ( max(0.d0,
 
402
     #                      xsecPDFr(0)-xsecPDFr(nps),
 
403
     #                      xsecPDFr(0)-xsecPDFr(nng)) )**2
 
404
            wgtpdfmax=wgtpdfmax+
 
405
     #                ( max(0.d0,
 
406
     #                      xsecPDFr(nps)-xsecPDFr(0),
 
407
     #                      xsecPDFr(nng)-xsecPDFr(0)) )**2
 
408
          enddo
 
409
          wgtpdfmin=wgtcentral-sqrt(wgtpdfmin)
 
410
          wgtpdfmax=wgtcentral+sqrt(wgtpdfmax)
 
411
 
 
412
c Restore default PDFs
 
413
          call InitPDF(izero)
 
414
 
 
415
        endif
 
416
 
 
417
c Keep track of the accumulated results:
 
418
        if (numscales.gt.0) then
 
419
           do ii=1,numscales
 
420
              do jj=1,numscales
 
421
                 xsecScale_acc(ii,jj)=xsecScale_acc(ii,jj)+wgtxsecmu(ii
 
422
     $                ,jj)/wgtref*XWGTUP
 
423
              enddo
 
424
           enddo
 
425
        endif
 
426
        if (nsets.gt.0) then
 
427
           do n=0,nsets
 
428
              xsecPDFr_acc(n)=xsecPDFr_acc(n)+wgtxsecPDF(n)/wgtref
 
429
     $             *XWGTUP
 
430
           enddo
 
431
        endif
 
432
 
 
433
c renormalize all the scale & PDF weights to have the same normalization
 
434
c as XWGTUP
 
435
        if(do_rwgt_scale)then
 
436
           do kr=1,3
 
437
              do kf=1,3
 
438
                 wgtxsecmu(kr,kf)=wgtxsecmu(kr,kf)/wgtref*XWGTUP
 
439
              enddo
 
440
           enddo
 
441
        endif
 
442
        if(do_rwgt_pdf)then
 
443
           do n=0,nsets
 
444
              wgtxsecPDF(n)=wgtxsecPDF(n)/wgtref*XWGTUP
 
445
           enddo
 
446
        endif
 
447
 
 
448
c Write event to disk:
 
449
        write(buff,201)'#aMCatNLO',iSorH_lhe,ifks_lhe,jfks_lhe,
 
450
     #                     fksfather_lhe,ipartner_lhe,
 
451
     #                     scale1_lhe,scale2_lhe,
 
452
     #                     isave,izero,izero,
 
453
     #          wgtcentral,wgtmumin,wgtmumax,wgtpdfmin,wgtpdfmax
 
454
 
 
455
        call write_lhef_event(ofile,
 
456
     &       NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
 
457
     &       IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
 
458
 
 
459
      enddo
 
460
 
 
461
      write(ofile,'(a)')'</LesHouchesEvents>'
 
462
      close(34)
 
463
      close(35)
 
464
 
 
465
c Write the accumulated results to a file
 
466
      open (unit=34,file='scale_pdf_dependence.dat',status='unknown')
 
467
      write (34,*) numscales**2
 
468
      if (numscales.gt.0) then
 
469
         write (34,*) ((xsecScale_acc(ii,jj),ii=1,numscales),jj=1
 
470
     $        ,numscales)
 
471
      else
 
472
         write (34,*) ''
 
473
      endif
 
474
      if (nsets.gt.0) then
 
475
         write (34,*) nsets + 1
 
476
         write (34,*) (xsecPDFr_acc(n),n=0,nsets)
 
477
      else
 
478
         write(34,*) nsets
 
479
         write (34,*) ''
 
480
      endif
 
481
      close(34)
 
482
 
 
483
 201  format(a9,1x,i1,4(1x,i2),2(1x,e14.8),1x,i1,2(1x,i2),5(1x,e14.8))
 
484
 
 
485
      end
 
486
 
 
487
 
 
488
 
 
489
      subroutine set_cms_stuff(icountevts)
 
490
      implicit none
 
491
      include "run.inc"
 
492
 
 
493
      integer icountevts
 
494
 
 
495
      double precision ybst_til_tolab,ybst_til_tocm,sqrtshat,shat
 
496
      common/parton_cms_stuff/ybst_til_tolab,ybst_til_tocm,
 
497
     #                        sqrtshat,shat
 
498
 
 
499
      double precision sqrtshat_ev,shat_ev
 
500
      common/parton_cms_ev/sqrtshat_ev,shat_ev
 
501
 
 
502
      double precision sqrtshat_cnt(-2:2),shat_cnt(-2:2)
 
503
      common/parton_cms_cnt/sqrtshat_cnt,shat_cnt
 
504
 
 
505
      double precision tau_ev,ycm_ev
 
506
      common/cbjrk12_ev/tau_ev,ycm_ev
 
507
 
 
508
      double precision tau_cnt(-2:2),ycm_cnt(-2:2)
 
509
      common/cbjrk12_cnt/tau_cnt,ycm_cnt
 
510
 
 
511
      double precision xbjrk_ev(2),xbjrk_cnt(2,-2:2)
 
512
      common/cbjorkenx/xbjrk_ev,xbjrk_cnt
 
513
 
 
514
c rapidity of boost from \tilde{k}_1+\tilde{k}_2 c.m. frame to lab frame --
 
515
c same for event and counterevents
 
516
c This is the rapidity that enters in the arguments of the sinh() and
 
517
c cosh() of the boost, in such a way that
 
518
c       y(k)_lab = y(k)_tilde - ybst_til_tolab
 
519
c where y(k)_lab and y(k)_tilde are the rapidities computed with a generic
 
520
c four-momentum k, in the lab frame and in the \tilde{k}_1+\tilde{k}_2 
 
521
c c.m. frame respectively
 
522
      ybst_til_tolab=-ycm_cnt(0)
 
523
      if(icountevts.eq.-100)then
 
524
c set Bjorken x's in run.inc for the computation of PDFs in auto_dsig
 
525
        xbk(1)=xbjrk_ev(1)
 
526
        xbk(2)=xbjrk_ev(2)
 
527
c shat=2*k1.k2 -- consistency of this assignment with momenta checked
 
528
c in phspncheck_nocms
 
529
        shat=shat_ev
 
530
        sqrtshat=sqrtshat_ev
 
531
c rapidity of boost from \tilde{k}_1+\tilde{k}_2 c.m. frame to 
 
532
c k_1+k_2 c.m. frame
 
533
        ybst_til_tocm=ycm_ev-ycm_cnt(0)
 
534
      else
 
535
c do the same as above for the counterevents
 
536
        xbk(1)=xbjrk_cnt(1,icountevts)
 
537
        xbk(2)=xbjrk_cnt(2,icountevts)
 
538
        shat=shat_cnt(icountevts)
 
539
        sqrtshat=sqrtshat_cnt(icountevts)
 
540
        ybst_til_tocm=ycm_cnt(icountevts)-ycm_cnt(0)
 
541
      endif
 
542
      return
 
543
      end
 
544
 
 
545