| 1 | subroutine histo_init(num_channel)
|
|---|
| 2 | implicit none
|
|---|
| 3 | c
|
|---|
| 4 | c argument
|
|---|
| 5 | c
|
|---|
| 6 | integer num_channel
|
|---|
| 7 | c
|
|---|
| 8 | c local
|
|---|
| 9 | c
|
|---|
| 10 | character*40 t1, t2, t3, t4, t0
|
|---|
| 11 | integer i
|
|---|
| 12 | c
|
|---|
| 13 | c gobal
|
|---|
| 14 | c
|
|---|
| 15 | double precision GeVbin, etabin, xbin, PTbin, THETAbin
|
|---|
| 16 | common/to_bin/ GeVbin, etabin, xbin
|
|---|
| 17 | c
|
|---|
| 18 | c
|
|---|
| 19 | call inihist
|
|---|
| 20 | c
|
|---|
| 21 | GeVbin=24d0
|
|---|
| 22 | etabin=0.1d0
|
|---|
| 23 | PTbin=10d0
|
|---|
| 24 | THETAbin=4d-2
|
|---|
| 25 | c
|
|---|
| 26 | t0='histo for one vegas iteration (prov)'
|
|---|
| 27 | t1='final state invariant mass '
|
|---|
| 28 | t2='fct error'
|
|---|
| 29 | CALL MBOOK(1,200,t0,GeVbin,100d0,1300d0)
|
|---|
| 30 | CALL MBOOK(1,-200,t0,GeVbin,100d0,1300d0)
|
|---|
| 31 | c
|
|---|
| 32 | t3='Top quark PT '
|
|---|
| 33 | CALL MBOOK(2,200,t0,PTbin,0d0,500d0)
|
|---|
| 34 | CALL MBOOK(2,-200,t0,PTbin,0d0,500d0)
|
|---|
| 35 | c
|
|---|
| 36 | t4='THETA angle '
|
|---|
| 37 | CALL MBOOK(3,200,t0,THETAbin,-1d0,1d0)
|
|---|
| 38 | CALL MBOOK(3,-200,t0,THETAbin,-1d0,1d0)
|
|---|
| 39 |
|
|---|
| 40 | do i=1,num_channel
|
|---|
| 41 | CALL MBOOK(1,i,t1,GeVbin,100d0,1300d0)
|
|---|
| 42 | CALL MBOOK(1,-1*i,t2,GeVbin,100d0,1300d0)
|
|---|
| 43 | CALL MBOOK(2,i,t3,PTbin,0d0,500d0)
|
|---|
| 44 | CALL MBOOK(2,-1*i,t2,PTbin,0d0,500d0)
|
|---|
| 45 | CALL MBOOK(3,i,t4,THETAbin,-1d0,1d0)
|
|---|
| 46 | CALL MBOOK(3,-1*i,t2,THETAbin,-1d0,1d0)
|
|---|
| 47 | enddo
|
|---|
| 48 | return
|
|---|
| 49 |
|
|---|
| 50 | end
|
|---|
| 51 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|---|
| 52 |
|
|---|
| 53 | subroutine FILL_plot(fct,weight)
|
|---|
| 54 | implicit none
|
|---|
| 55 | include '../../nexternal.inc'
|
|---|
| 56 | include '../../genps.inc'
|
|---|
| 57 | c
|
|---|
| 58 | c arguments
|
|---|
| 59 | c
|
|---|
| 60 | double precision weight,fct
|
|---|
| 61 | c
|
|---|
| 62 | c local
|
|---|
| 63 | c
|
|---|
| 64 | double precision Ptot(0:3), minv
|
|---|
| 65 | integer i,nu
|
|---|
| 66 |
|
|---|
| 67 | double precision Ptot1(0:3),Ptot2(0:3),PT1,PT2
|
|---|
| 68 | double precision BETA,GAM,PZ_CM,THETA
|
|---|
| 69 | c
|
|---|
| 70 | c common
|
|---|
| 71 | c
|
|---|
| 72 | double precision momenta(0:3,-max_branch:2*nexternal) ! records the momenta of external/intermediate legs (MG order)
|
|---|
| 73 | double precision mvir2(-max_branch:2*nexternal) ! records the sq invariant masses of intermediate particles (MG order)
|
|---|
| 74 | common /to_diagram_kin/ momenta, mvir2
|
|---|
| 75 |
|
|---|
| 76 | double precision PT
|
|---|
| 77 | external PT
|
|---|
| 78 |
|
|---|
| 79 |
|
|---|
| 80 | c---
|
|---|
| 81 | c Begin code for variable 1 (invariant mass of the top quark pair)
|
|---|
| 82 | c---
|
|---|
| 83 |
|
|---|
| 84 | do nu=0,3
|
|---|
| 85 | Ptot(nu)=0d0
|
|---|
| 86 | do i=3,nexternal
|
|---|
| 87 | Ptot(nu)=Ptot(nu)+momenta(nu,i)
|
|---|
| 88 | enddo
|
|---|
| 89 | enddo
|
|---|
| 90 | minv=dsqrt(Ptot(0)**2-Ptot(1)**2-Ptot(2)**2-Ptot(3)**2)
|
|---|
| 91 |
|
|---|
| 92 | CALL MFILL(1,200, minv, fct*weight)
|
|---|
| 93 | CALL MFILL(1,-200, minv, fct**2*weight)
|
|---|
| 94 |
|
|---|
| 95 | c---
|
|---|
| 96 | c Begin code for variable 2 (PT of the top quark)
|
|---|
| 97 | c---
|
|---|
| 98 |
|
|---|
| 99 | do nu=0,3
|
|---|
| 100 | Ptot1(nu)=0d0
|
|---|
| 101 | Ptot2(nu)=0d0
|
|---|
| 102 | do i=3,5
|
|---|
| 103 | Ptot1(nu)=Ptot1(nu)+momenta(nu,i)
|
|---|
| 104 | Ptot2(nu)=Ptot2(nu)+momenta(nu,i+3)
|
|---|
| 105 | enddo
|
|---|
| 106 | enddo
|
|---|
| 107 | PT1=PT(Ptot1)
|
|---|
| 108 | PT2=PT(Ptot2)
|
|---|
| 109 |
|
|---|
| 110 | CALL MFILL(2,200,PT1,fct*weight)
|
|---|
| 111 | CALL MFILL(2,-200,PT1,fct**2*weight)
|
|---|
| 112 |
|
|---|
| 113 | c CALL MFILL()
|
|---|
| 114 | c CALL MFILL()
|
|---|
| 115 |
|
|---|
| 116 | c---
|
|---|
| 117 | c Begin code for variable 3 (angle of the top quark in the center of mass referentiel)
|
|---|
| 118 | c---
|
|---|
| 119 |
|
|---|
| 120 |
|
|---|
| 121 | BETA = (momenta(0,1)-momenta(0,2))/(momenta(0,1)+momenta(0,2))
|
|---|
| 122 | GAM = dsqrt(1-BETA**2)
|
|---|
| 123 |
|
|---|
| 124 | PZ_CM = GAM*Ptot1(3)-GAM*BETA*Ptot1(0)
|
|---|
| 125 | THETA = PZ_CM/(Ptot1(1)**2+Ptot1(2)**2+PZ_CM**2)
|
|---|
| 126 |
|
|---|
| 127 | CALL MFILL(3,200,THETA,fct*weight)
|
|---|
| 128 | CALL MFILL(3,-200,THETA,fct**2*weight)
|
|---|
| 129 |
|
|---|
| 130 | return
|
|---|
| 131 | end
|
|---|
| 132 |
|
|---|
| 133 |
|
|---|
| 134 |
|
|---|
| 135 | cCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcccccccccccccccccccccccccccc
|
|---|
| 136 | subroutine histo_final(num_channel)
|
|---|
| 137 | c
|
|---|
| 138 | c argument
|
|---|
| 139 | c
|
|---|
| 140 | integer num_per
|
|---|
| 141 | c
|
|---|
| 142 | c local
|
|---|
| 143 | c
|
|---|
| 144 | integer i,num_iter,j,k
|
|---|
| 145 | character*30 BUFFER
|
|---|
| 146 | character istring
|
|---|
| 147 | c
|
|---|
| 148 | c gobal
|
|---|
| 149 | c
|
|---|
| 150 | double precision GeVbin, etabin, xbin
|
|---|
| 151 | common/to_bin/ GeVbin, etabin, xbin
|
|---|
| 152 | c
|
|---|
| 153 | c include histo globlal
|
|---|
| 154 | c
|
|---|
| 155 | include 'dbook.inc'
|
|---|
| 156 |
|
|---|
| 157 | do k=1,NHistovar
|
|---|
| 158 | do i=1,num_channel
|
|---|
| 159 | c if (i.eq.1) istring='1'
|
|---|
| 160 | c if (i.eq.2) istring='2'
|
|---|
| 161 |
|
|---|
| 162 | call MFINAL(k, i)
|
|---|
| 163 | call MFINAL(k, -i)
|
|---|
| 164 | call MOPERA(k, -1*i, 'I',-1*i,-1*i,1d0,1d0) !contains final error square
|
|---|
| 165 | call MOPERA(k, i, '*',-i,i,1d0,1d0) !contains final estimation of integral
|
|---|
| 166 |
|
|---|
| 167 | if(i.ne.1) then
|
|---|
| 168 | call MOPERA(k,i, '+',1,1,1d0,1d0)
|
|---|
| 169 | call MOPERA(k, -i, '+',-1,-1,1d0,1d0)
|
|---|
| 170 | endif
|
|---|
| 171 |
|
|---|
| 172 | call MFINAL(k, 1)
|
|---|
| 173 | call MOPERA(k, -1,'R',-1,-1,1d0,1d0) ! contains the error for the sum of all permutations
|
|---|
| 174 | call MFINAL(k, -1)
|
|---|
| 175 | enddo
|
|---|
| 176 | buffer='plot_00.dat'
|
|---|
| 177 | if (k.ge.10)then
|
|---|
| 178 | write(buffer(6:7),'(I2)') k
|
|---|
| 179 | else
|
|---|
| 180 | write(buffer(7:7),'(I1)') k
|
|---|
| 181 | endif
|
|---|
| 182 |
|
|---|
| 183 | OPEN(UNIT=98,file=buffer,STATUS='UNKNOWN')
|
|---|
| 184 | buffer=buffer(1:7)//'.top'
|
|---|
| 185 | OPEN(UNIT=99,file=buffer,STATUS='UNKNOWN')
|
|---|
| 186 | CALL MPRINT(k, 1)
|
|---|
| 187 | CALL MTOP(k, 1, 50,' M_INV (GeV) '
|
|---|
| 188 | & ,'d P/d M_inv ','log')
|
|---|
| 189 |
|
|---|
| 190 | close(98)
|
|---|
| 191 | close(99)
|
|---|
| 192 |
|
|---|
| 193 |
|
|---|
| 194 | buffer='err_'//buffer(1:7)//'.dat'
|
|---|
| 195 | OPEN(UNIT=98,file=buffer,STATUS='UNKNOWN')
|
|---|
| 196 | buffer=buffer(1:11)//'.top'
|
|---|
| 197 | OPEN(UNIT=99,file=buffer,STATUS='UNKNOWN')
|
|---|
| 198 | CALL MPRINT(k, -1)
|
|---|
| 199 | CALL MTOP(k, -1,50,' M_INV (GeV) '
|
|---|
| 200 | & ,'d P/d M_inv ','log')
|
|---|
| 201 | c
|
|---|
| 202 | close(98)
|
|---|
| 203 | close(99)
|
|---|
| 204 | enddo
|
|---|
| 205 | return
|
|---|
| 206 | end
|
|---|
| 207 |
|
|---|
| 208 | CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
|
|---|
| 209 | subroutine histo_combine_iteration(channel_pos)
|
|---|
| 210 | implicit none
|
|---|
| 211 | c
|
|---|
| 212 | c argument
|
|---|
| 213 | c
|
|---|
| 214 | integer channel_pos
|
|---|
| 215 | c
|
|---|
| 216 | c local
|
|---|
| 217 | c
|
|---|
| 218 | integer i,num_iter,j,k
|
|---|
| 219 | character*30 BUFFER
|
|---|
| 220 | character istring
|
|---|
| 221 | c
|
|---|
| 222 | c gobal
|
|---|
| 223 | c
|
|---|
| 224 | double precision GeVbin, etabin, xbin
|
|---|
| 225 | common/to_bin/ GeVbin, etabin, xbin
|
|---|
| 226 | c
|
|---|
| 227 | c include histo globlal
|
|---|
| 228 | c
|
|---|
| 229 | include 'dbook.inc'
|
|---|
| 230 | c call MFINAL(i)
|
|---|
| 231 | c call MFINAL(-1*i)
|
|---|
| 232 |
|
|---|
| 233 | do k=1,NHistoVar
|
|---|
| 234 |
|
|---|
| 235 | call MFINAL(k, 200)
|
|---|
| 236 | call MFINAL(k, -200)
|
|---|
| 237 |
|
|---|
| 238 | write(*,*) 'combine for histo', channel_pos
|
|---|
| 239 | call MOPERA(k, 200,'V',-200,-200,1d0,1d0)
|
|---|
| 240 | call MOPERA(k, -200,'S',0,-200,1d0,1d0)
|
|---|
| 241 | call MOPERA(k, 200,'/',-200,200,1d0,1d0)
|
|---|
| 242 | call MOPERA(k, channel_pos,'+',200,channel_pos,1d0,1d0)
|
|---|
| 243 |
|
|---|
| 244 | !define a unity histogram in 200 (after simply redifine to zero)
|
|---|
| 245 | call MOPERA(k, -200,'I',-200,-200,1d0,1d0)
|
|---|
| 246 | call MOPERA(k, -1*channel_pos,'+',-200,-1*channel_pos,1d0,1d0)
|
|---|
| 247 |
|
|---|
| 248 | call MOPERA(k, 200,'+',200,200,0d0,0d0)
|
|---|
| 249 | call MOPERA(k, -200,'+',-200,-200,0d0,0d0)
|
|---|
| 250 |
|
|---|
| 251 | enddo
|
|---|
| 252 | return
|
|---|
| 253 | end
|
|---|
| 254 |
|
|---|
