DifferentialWeight: topgraph.f

File topgraph.f, 6.2 KB (added by anonymous, 7 years ago)
Line 
1      subroutine histo_init(num_channel)
2      implicit none
3c
4c     argument
5c
6      integer num_channel
7c
8c    local
9c
10      character*40 t1, t2, t3, t4, t0
11      integer i
12c
13c     gobal
14c
15      double precision GeVbin, etabin, xbin, PTbin, THETAbin
16      common/to_bin/ GeVbin, etabin, xbin
17c
18c
19      call inihist
20c
21      GeVbin=24d0
22      etabin=0.1d0
23      PTbin=10d0
24      THETAbin=4d-2
25c
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)
31c
32      t3='Top quark PT '
33      CALL MBOOK(2,200,t0,PTbin,0d0,500d0)
34      CALL MBOOK(2,-200,t0,PTbin,0d0,500d0)
35c
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
51ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
52
53      subroutine FILL_plot(fct,weight)
54      implicit none
55      include '../../nexternal.inc'
56      include '../../genps.inc'
57c
58c     arguments
59c
60      double precision weight,fct
61c
62c     local
63c
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
69c
70c     common
71c
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
80c---
81c Begin code for variable 1 (invariant mass of the top quark pair)
82c---
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
95c---
96c Begin code for variable 2 (PT of the top quark)
97c---
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
113c      CALL MFILL()
114c      CALL MFILL()
115
116c---
117c Begin code for variable 3 (angle of the top quark in the center of mass referentiel)
118c---
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
135cCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcccccccccccccccccccccccccccc
136      subroutine histo_final(num_channel)
137c
138c     argument
139c   
140      integer num_per
141c
142c     local
143c     
144      integer i,num_iter,j,k
145      character*30 BUFFER
146      character istring
147c
148c     gobal
149c
150      double precision GeVbin, etabin, xbin
151      common/to_bin/ GeVbin, etabin, xbin
152c     
153c     include histo globlal
154c     
155      include 'dbook.inc'
156
157      do k=1,NHistovar
158         do i=1,num_channel
159c     if (i.eq.1) istring='1'
160c     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')
201c     
202         close(98)
203         close(99)
204      enddo
205      return
206      end
207
208CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
209      subroutine histo_combine_iteration(channel_pos)
210      implicit none
211c
212c     argument
213c   
214      integer channel_pos
215c
216c     local
217c     
218      integer i,num_iter,j,k
219      character*30 BUFFER
220      character istring
221c
222c     gobal
223c
224      double precision GeVbin, etabin, xbin
225      common/to_bin/ GeVbin, etabin, xbin
226c     
227c     include histo globlal
228c     
229      include 'dbook.inc'
230c      call MFINAL(i)
231c      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