JHUGen MELA  v2.4.1
Matrix element calculations as used in JHUGen. MELA is an important tool that was used for the Higgs boson discovery and for precise measurements of its structure and interactions. Please see the website https://spin.pha.jhu.edu/ and papers cited there for more details, and kindly cite those papers when using this code.
democache.f90
Go to the documentation of this file.
1 !!
2 !! File democache.f90 is part of COLLIER
3 !! - A Complex One-Loop Library In Extended Regularizations
4 !!
5 !! Copyright (C) 2015, 2016 Ansgar Denner, Stefan Dittmaier, Lars Hofer
6 !!
7 !! COLLIER is licenced under the GNU GPL version 3, see COPYING for details.
8 !!
9 
10 program democache
11 
12  use collier
13 
14  implicit none
15 
16  double complex :: mominv6(15), masses2_6(0:5), masses2_6b(0:5)
17  double complex :: mominv5(10), masses2_5(0:4)
18  double complex :: mominv4(6), masses2_4(0:3)
19  double complex :: mominv3(3), masses2_3(0:2)
20  double complex :: mominv2(1), masses2_2(0:1)
21  double complex, allocatable :: t1coeff(:),t1coeffuv(:)
22  double complex, allocatable :: t2coeff(:),t2coeffuv(:)
23  double complex, allocatable :: t3coeff(:),t3coeffuv(:)
24  double complex, allocatable :: t4coeff(:),t4coeffuv(:)
25  double complex, allocatable :: t5coeff(:),t5coeffuv(:)
26  double complex, allocatable :: t6coeff(:),t6coeffuv(:)
27  double precision, allocatable :: t1err(:),t2err(:),t3err(:),t4err(:)
28  double precision, allocatable :: t5err(:),t6err(:)
29  real :: starttime,endtime
30  integer :: rank2,rank3,rank4,rank5,rank6,nmax,rmax
31  integer :: i,k,nevents,nloops,icache,mycase
32 
33 ! Nmax = maximal degree N of N-point function
34  nmax = 6
35  rmax = 4
36 
37 ! Initialize COLLIER
38 ! the default directory for diagnostic messages is "output_cll"
39  call init_cll(nmax,rmax,"output_cll")
40 
41 ! INITIALIZE CACHE nr. "icache" for up to Nmax-point functions
42  icache = 1
43  call initcachesystem_cll(icache,nmax)
44 
45 
46 ! allocate fields for coefficients
47  rank2=3
48  rank3=4
49  rank4=4
50  rank5=4
51  rank6=4
52  allocate(t2coeff(getnc_cll(2,rank2)))
53  allocate(t2coeffuv(getnc_cll(2,rank2)))
54  allocate(t2err(0:rank2))
55  allocate(t3coeff(getnc_cll(3,rank3)))
56  allocate(t3coeffuv(getnc_cll(3,rank3)))
57  allocate(t3err(0:rank3))
58  allocate(t4coeff(getnc_cll(4,rank4)))
59  allocate(t4coeffuv(getnc_cll(4,rank4)))
60  allocate(t4err(0:rank4))
61  allocate(t5coeff(getnc_cll(5,rank5)))
62  allocate(t5coeffuv(getnc_cll(5,rank5)))
63  allocate(t5err(0:rank5))
64  allocate(t6coeff(getnc_cll(6,rank6)))
65  allocate(t6coeffuv(getnc_cll(6,rank6)))
66  allocate(t6err(0:rank6))
67 
68 ! set some values for masses
69  masses2_6(0) = dcmplx(1d2,-1d0)
70  masses2_6(1) = dcmplx(1d2,-1d0)
71  masses2_6(2) = dcmplx(0d0,0d0)
72  masses2_6(3) = dcmplx(1d1,0d0)
73  masses2_6(4) = dcmplx(1d1,0d0)
74  masses2_6(5) = dcmplx(0d0,0d0)
75 
76 
77 ! call InitMonitoring_cll
78 
79 ! issue a bunch of calls of COLLIER for different tensor integrals
80 ! for several "events"
81 
82  nevents = 2000
83 
84  nloops = 2
85 
86  write(*,'(/(a))') ' *** Demonstration of cache system ***'
87  write(*,'(/(a))') ' NOTE: This example merely illustrates the usage of the cache system.'
88  write(*,'((a))') ' The results are not representative for the performance'
89  write(*,'((a))') ' of Collier, neither wrt to speed nor precision.'
90  write(*,'(/a21,i10/)') ' Number of "Events": ',nevents
91  write(*,'(/a39,i4/)') ' Number of "identical calls/integral": ',nloops
92 
93  do mycase=1,4
94 
95  if (mycase.eq.1) then
96 ! use COLI branch
97  call setmode_cll(1)
98 ! initialize cache nr. icache for tensor integrals up to N=6
99  call initcachesystem_cll(icache,6)
100 
101  write(*,'(/(a))') ' COLI branch with cache: '
102 
103  elseif (mycase.eq.2) then
104 ! use COLI branch
105  call setmode_cll(1)
106 ! switch off cache
107  call switchoffcachesystem_cll
108 
109  write(*,'(/(a))') ' COLI branch without cache: '
110 
111  else if (mycase.eq.3) then
112 ! use DD branch
113  call setmode_cll(2)
114 ! reinitialize cache nr. icache
115  call initcachesystem_cll(icache,6)
116 
117  write(*,'(/(a))') ' DD branch with cache: '
118 
119  else if (mycase.eq.4) then
120 ! use DD branch
121  call setmode_cll(2)
122 ! switch off cache
123  call switchoffcachesystem_cll
124 
125  write(*,'(/(a))') ' DD branch without cache '
126 
127  end if
128 
129  call cpu_time(starttime)
130 
131 ! write(*,*) 'start',starttime
132 
133  eventloop: do i=1,nevents
134 
135 ! get momenta and invariants (typically from Monte Carlo generator)
136  call getinvariants(mominv6,i)
137 
138 ! TELL COLLIER START OF NEW EVENT for cache "icache"
139  call initevent_cll(icache)
140 
141 ! call bunch of tensor integrals or coefficients
142 
143 ! loop to generate identical master calls
144  testloop: do k=1,nloops
145 
146  call tn_cll(t6coeff,t6coeffuv,mominv6,masses2_6,6,rank6,t6err)
147 
148 ! pick arguments for tensor integrals of lower rank
149  mominv5 = submominv(6,0,mominv6)
150  masses2_5 = submasses(6,0,masses2_6)
151  call tn_cll(t5coeff,t5coeffuv,mominv5,masses2_5,5,rank5,t5err)
152 
153  mominv4 = submominv(5,0,mominv5)
154  masses2_4 = submasses(5,0,masses2_5)
155  mominv3 = submominv(4,0,mominv4)
156  masses2_3 = submasses(4,0,masses2_4)
157  call tn_cll(t3coeff,t3coeffuv,mominv3,masses2_3,3,rank3,t3err)
158  call tn_cll(t4coeff,t4coeffuv,mominv4,masses2_4,4,rank4,t4err)
159 
160 ! change one of the masses
161  masses2_6b(0) = (0d0,0d0)
162  masses2_6b(1:5) = masses2_6(1:5)
163  call tn_cll(t6coeff,t6coeffuv,mominv6,masses2_6b,6,rank6,t6err)
164 
165  mominv5 = submominv(6,1,mominv6)
166  masses2_5 = submasses(6,1,masses2_6b)
167  mominv4 = submominv(5,2,mominv5)
168  masses2_4 = submasses(5,2,masses2_5)
169  mominv3 = submominv(4,0,mominv4)
170  masses2_3 = submasses(4,0,masses2_4)
171  mominv2 = submominv(3,0,mominv3)
172  masses2_2 = submasses(3,0,masses2_3)
173  call tn_cll(t2coeff,t2coeffuv,mominv2,masses2_2,2,rank2,t2err)
174  call tn_cll(t5coeff,t5coeffuv,mominv5,masses2_5,5,rank5,t5err)
175  call tn_cll(t4coeff,t4coeffuv,mominv4,masses2_4,4,rank4,t4err)
176 
177  end do testloop
178 
179  end do eventloop
180 
181  call cpu_time(endtime)
182 
183  write(*,'(/a29,f10.1/)') ' Needed CPU time in seconds: ',endtime-starttime
184 
185  end do
186 
187 ! call PrintStatistics_cll
188 
189  deallocate(t2coeff,t2coeffuv,t2err)
190  deallocate(t3coeff,t3coeffuv,t3err)
191  deallocate(t4coeff,t4coeffuv,t4err)
192  deallocate(t5coeff,t5coeffuv,t5err)
193  deallocate(t6coeff,t6coeffuv,t6err)
194 
195 contains
196 
197  subroutine getinvariants(MomInv6,ievent)
198 
199  implicit none
200 
201  double precision p(1:6,0:3)
202  double complex MomInv6(1:15)
203  double precision :: s(0:5,0:5),fact=47d0/49d0
204  integer :: N=6,k,i,newmom,ievent
205  save s,fact,n
206 
207  newmom=0
208 
209  if (ievent.eq.1) then
210  newmom=1
211 ! use momentum set from Racoon
212  p(1,0) = -0.9500000000000000d+02
213  p(1,1) = 0.0000000000000000d+00
214  p(1,2) = 0.0000000000000000d+00
215  p(1,3) = 0.9500000000000000d+02
216 
217  p(2,0) = -0.9500000000000000d+02
218  p(2,1) = 0.0000000000000000d+00
219  p(2,2) = 0.0000000000000000d+00
220  p(2,3) = -0.9500000000000000d+02
221 
222  p(3,0) = 0.2046111001757171d+02
223  p(3,1) = 0.1057734233089455d+02
224  p(3,2) = -0.2324961261504543d+01
225  p(3,3) = 0.1736005205921753d+02
226 
227  p(4,0) = 0.3558305163378869d+01
228  p(4,1) = 0.1436222934374051d+01
229  p(4,2) = -0.2174258125294355d+01
230  p(4,3) = -0.2423097382091398d+01
231 
232  p(5,0) = 0.8154540918019539d+02
233  p(5,1) = -0.5230395944682889d+02
234  p(5,2) = 0.3083642435466509d+02
235  p(5,3) = 0.5443403822581044d+02
236 
237  p(6,0) = 0.8443517563885433d+02
238  p(6,1) = 0.4029039418156027d+02
239  p(6,2) = -0.2633720496786619d+02
240  p(6,3) = -0.6937099290293661d+02
241 
242 
243  else if (ievent.eq.50) then
244  newmom=1
245 ! set with small momentum 1.9 versus 3.0 seconds
246  p(1,0) = -250.0000000000000d0
247  p(1,1) = 0.0000000000000000d+00
248  p(1,2) = 0.0000000000000000d+00
249  p(1,3) = -250.0000000000000d0
250 
251  p(2,0) = -250.0000000000000d0
252  p(2,1) = 0.0000000000000000d+00
253  p(2,2) = 0.0000000000000000d+00
254  p(2,3) = 250.0000000000000d0
255 
256  p(3,0) = 239.91697591218110d0
257  p(3,1) = 31.640107125154238d0
258  p(3,2) = -116.12928664566658d0
259  p(3,3) = -207.54047252312847d0
260 
261  p(4,0) = 10.276659578244022d0
262  p(4,1) = 1.2584629023924521d0
263  p(4,2) = 10.131487190720421d0
264  p(4,3) = -1.1742957526504370d0
265 
266  p(5,0) = 181.18800109198216d0
267  p(5,1) = -35.098898474333737d0
268  p(5,2) = 109.26326872331617d0
269  p(5,3) = 140.20947604742324d0
270 
271  p(6,0) = 68.618363417592718d0
272  p(6,1) = 2.2003284467870472d0
273  p(6,2) = -3.2654692683700119d0
274  p(6,3) = 68.505292228355671d0
275 
276  else if (ievent.eq.100) then
277  newmom=1
278  p(1,0) = -100.0000000000000d0
279  p(1,1) = 0.000000000000000d0
280  p(1,2) = 0.000000000000000d0
281  p(1,3) = -100.0000000000000d0
282 
283  p(2,0) = -100.0000000000000d0
284  p(2,1) = 0.000000000000000d0
285  p(2,2) = 0.000000000000000d0
286  p(2,3) = 100.0000000000000d0
287 
288  p(5,0) = 42.13045937773281d0
289  p(5,1) = 38.71045260055064d0
290  p(5,2) = 16.54764683624944d0
291  p(5,3) = 1.628450497276325d0
292 
293  p(6,0) = 73.31699208541762d0
294  p(6,1) = -66.17895647220249d0
295  p(6,2) = 31.53024134467396d0
296  p(6,3) = -1.253367244715983d0
297 
298  p(4,0) = 69.60509150033207d0
299  p(4,1) = 34.42307953372518d0
300  p(4,2) = -60.31848303205667d0
301  p(4,3) = 4.647683605402023d0
302 
303  p(3,0) = 14.94745703651754d0
304  p(3,1) = -6.954575662073326d0
305  p(3,2) = 12.24059485113326d0
306  p(3,3) = -5.022766857962365d0
307  else if (ievent.eq.150) then
308  newmom=1
309  p(1,0) = -100.0000000000000d0
310  p(1,1) = 0.0000000000000000d+00
311  p(1,2) = 0.0000000000000000d+00
312  p(1,3) = -100.0000000000000d0
313 
314  p(2,0) = -100.0000000000000d0
315  p(2,1) = 0.0000000000000000d+00
316  p(2,2) = 0.0000000000000000d+00
317  p(2,3) = 100.0000000000000d0
318 
319  p(3,0) = 72.02905523601244d0
320  p(3,1) = 32.33193950844802d0
321  p(3,2) = 64.10893732125734d0
322  p(3,3) = 5.733641195059672d0
323 
324  p(4,0) = 25.07519219892066d0
325  p(4,1) = 2.906699172283076d0
326  p(4,2) = -24.65281899347465d0
327  p(4,3) = 3.543286526606821d0
328 
329  p(5,0) = 60.19609558319355d0
330  p(5,1) = 2.282219923475870d0
331  p(5,2) = -59.72688022443283d0
332  p(5,3) = -7.145710205299344d0
333 
334  p(6,0) = 42.69965698187335d0
335  p(6,1) = -37.52085860420697d0
336  p(6,2) = 20.27076189665014d0
337  p(6,3) = -2.131217516367149d0
338  else if (ievent.eq.200) then
339  newmom=1
340  p(1,0) = -100.0000000000000d0
341  p(1,1) = 0.0000000000000000d+00
342  p(1,2) = 0.0000000000000000d+00
343  p(1,3) = -100.0000000000000d0
344 
345  p(2,0) = -100.0000000000000d0
346  p(2,1) = 0.0000000000000000d+00
347  p(2,2) = 0.0000000000000000d+00
348  p(2,3) = 100.0000000000000d0
349 
350  p(3,0) = 69.63925440091566d0
351  p(3,1) = -68.33683461638704d0
352  p(3,2) = -3.803463998624776d0
353  p(3,3) = -12.85443307731553d0
354 
355  p(4,0) = 34.49634785978589d0
356  p(4,1) = 21.54290146454439d0
357  p(4,2) = -15.60187514686951d0
358  p(4,3) = 21.96549348532329d0
359 
360  p(5,0) = 66.20628864972379d0
361  p(5,1) = 65.60907499667245d0
362  p(5,2) = -3.500460884404577d0
363  p(5,3) = -8.152834381263276d0
364 
365  p(6,0) = 29.65810908957467d0
366  p(6,1) = -18.81514184482980d0
367  p(6,2) = 22.90580002989887d0
368  p(6,3) = -0.9582260267444891d0
369 
370  else if (ievent.eq.250) then
371  newmom=1
372  p( 1, 0) = -159.2798637127597488d0
373  p( 1, 1) = 0.0000000000000000d0
374  p( 1, 2) = 0.0000000000000000d0
375  p( 1, 3) = -159.2798637127597488d0
376  p( 2, 0) = -237.6261878491357891d0
377  p( 2, 1) = 0.0000000000000000d0
378  p( 2, 2) = 0.0000000000000000d0
379  p( 2, 3) = 237.6261878491357891d0
380  p( 3, 0) = 133.6018149124976446d0
381  p( 3, 1) = 60.5134748218625731d0
382  p( 3, 2) = 112.6312302087531094d0
383  p( 3, 3) = -38.7526811273741174d0
384  p( 4, 0) = 69.9256648527542950d0
385  p( 4, 1) = 66.6463990867762135d0
386  p( 4, 2) = 7.6538437283820429d0
387  p( 4, 3) = -19.7300473909055754d0
388  p( 5, 0) = 85.7314709474216841d0
389  p( 5, 1) = -38.6625297933730749d0
390  p( 5, 2) = -63.4432116127613241d0
391  p( 5, 3) = -42.7791164126088006d0
392  p( 6, 0) = 107.6471008492219283d0
393  p( 6, 1) = -88.4973441152657188d0
394  p( 6, 2) = -56.8418623243738281d0
395  p( 6, 3) = 22.9155207945124531d0
396 
397 
398  end if
399 
400  if (newmom.eq.1) then
401  ! define invariants
402  n=6
403  do k=0,n-1
404  s(modulo(k+1,n),k) = p(k+1,0)**2
405  do i=1,3
406  s(modulo(k+1,n),k) = s(modulo(k+1,n),k) - p(k+1,i)**2
407  end do
408  if(abs(s(modulo(k+1,n),k)).lt.1d-14* abs(p(k+1,0))**2) s(modulo(k+1,n),k) = 0d0
409 
410  s(k,modulo(k+1,n)) = s(modulo(k+1,n),k)
411  end do
412 
413  do k=0,n-1
414  s(modulo(k+2,n),k) = (p(modulo(k+1,n)+1,0) + p(k+1,0))**2
415  do i=1,3
416  s(modulo(k+2,n),k) = s(modulo(k+2,n),k) &
417  - (p(modulo(k+1,n)+1,i) + p(k+1,i))**2
418  end do
419  s(k,modulo(k+2,n)) = s(modulo(k+2,n),k)
420  end do
421  do k=0,(n-1)/2
422  s(modulo(k+3,n),k) = (p(modulo(k+2,n)+1,0) + p(modulo(k+1,n)+1,0) + p(k+1,0))**2
423  do i=1,3
424  s(modulo(k+3,n),k) = s(modulo(k+3,n),k) &
425  - (p(modulo(k+2,n)+1,i) + p(modulo(k+1,n)+1,i) + p(k+1,i))**2
426  end do
427  s(k,modulo(k+3,n)) = s(modulo(k+3,n),k)
428  end do
429 
430  do k=1,n
431  mominv6(k) = s(modulo(k,n),k-1)
432  end do
433  do k=1,n
434  mominv6(k+n) = s(modulo(k+1,n),k-1)
435  end do
436  do k=1,n/2
437  mominv6(k+2*n) = s(modulo(k+2,n),k-1)
438  end do
439  fact = 47d0/97d0
440 
441  else
442 ! rescale invariants
443  fact = 4*fact*(1-fact)
444 
445  do k=1,n
446  mominv6(k) = s(modulo(k,n),k-1)*fact
447  end do
448  do k=1,n
449  mominv6(k+n) = s(modulo(k+1,n),k-1)*fact
450  end do
451  do k=1,n/2
452  mominv6(k+2*n) = s(modulo(k+2,n),k-1)*fact
453  end do
454  end if
455 
456  end subroutine getinvariants
457 
458 
459 
460 
461 
462  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
463  ! function SubMasses
464  !
465  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
466 
467  function submasses(N,k,masses)
468 
469  integer :: n,k,i
470  double complex :: masses(0:n-1), submasses(0:n-2)
471 
472  do i=0,n-2
473  if (i.lt.k) then
474  submasses(i) = masses(i)
475  else
476  submasses(i) = masses(i+1)
477  end if
478  end do
479 
480  end function
481 
482 
483 
484 
485 
486  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
487  ! function SubMomInv
488  ! ******************************************
489  ! * rewritten by Ansgar Denner 25.03.15 *
490  ! ******************************************
491  !
492  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
493 
494  function submominv(N,k,MomInv)
495 
496  integer, intent(in) :: n,k
497  integer :: i,cnt,moms,limit
498  double complex :: mominv(n*(n-1)/2), submominv((n-1)*(n-2)/2)
499 
500  if (n.eq.2) return
501 
502  cnt = 1
503  do moms=1,(n-3)/2
504  do i=1,n
505  if (i.ne.k+1) then
506  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
507  submominv(cnt) = mominv(n*moms+i)
508  else
509  submominv(cnt) = mominv(n*(moms-1)+i)
510  end if
511 
512  cnt = cnt+1
513  end if
514  end do
515  end do
516 
517  if (mod(n,2).eq.1) then
518  moms = (n-1)/2
519  if(k.le.(n-1)/2) then
520  limit = (n+1)/2
521  else
522  limit = (n-1)/2
523  end if
524  do i=1,limit
525  if (i.ne.k+1) then
526  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
527  submominv(cnt) = mominv(n*moms+i-(n-1)/2)
528  else
529  submominv(cnt) = mominv(n*(moms-1)+i)
530  end if
531 
532  cnt = cnt+1
533  end if
534  end do
535 
536  else
537  do i=1,n
538  if (i.ne.k+1) then
539  moms = (n/2-1)
540  if ((k.ge.i).and.(i.ge.k-moms+1).or.(n+k.ge.i).and.(i.ge.n+k-moms+1)) then
541  if (i.gt.n/2) then
542  submominv(cnt) = mominv(n*moms-n/2+i)
543  else
544  submominv(cnt) = mominv(n*moms+i)
545  endif
546  else
547  submominv(cnt) = mominv(n*(moms-1)+i)
548  end if
549 
550  cnt = cnt+1
551  end if
552  end do
553 
554  end if
555 
556 
557  end function
558 
559 
560 end program democache
submasses
double complex function, dimension(0:n-2) submasses(N, k, masses)
Definition: democache.f90:468
collier
Definition: COLLIER.F90:30
submominv
double complex function, dimension((n-1) *(n-2)/2) submominv(N, k, MomInv)
Definition: democache.f90:495
democache
program democache
Definition: democache.f90:10
getinvariants
subroutine getinvariants(N, MomInv, MomVec)
Definition: demo.f90:3262
Generated by   doxygen 1.8.17