All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
syjetf2.f
Go to the documentation of this file.
1 CDECK ID>, SYJETF2.
2 ************************************************************************
3  SUBROUTINE syjetf2(IEXAM,N0JET,NPART,QIN,IJ,FRAC,IJOUT,QJET)
4 * Particle Re-association Routine (Mode=2)
5 * Jet core is always defined by the original jet
6 *
7 * Authour : Satoru Yamashita
8 * Creat : 1-Oct-97
9 * Mod : 1-Feb-98 S.Yamashita buffer size increased
10 *
11 * Inputs:
12 * IEXAM ; Jet-particle Re-association method 1-6
13 * EXAM(1) = THETA
14 * EXAM(2) = PART(4)*PJSTMP(4,J)-PART(1)*PJSTMP(1,J)
15 * + -PART(2)*PJSTMP(2,J)-PART(3)*PJSTMP(3,J)
16 * EXAM(3) = DURHAM
17 * EXAM(4) = E0JADE
18 * EXAM(5) = Geneva Scheme
19 * EXAM(6) = EJADE ! it makes too bad cross talk
20 * N0JET ; Number of Jets
21 * NPART ; Number of particles
22 * QIN ; particles 5-momenta
23 * IJ ; Original Jet assignment for "particle"s
24 * FRAC ; Energy Fraction for the core
25 * Outputs:
26 * IJOUT ; New Jet assignment for "particle"s
27 * QJET ; New Jet 6-momenta
28 *
29  implicit none
30  INTEGER iexam,n,i,j,k,ii
31  REAL emin2,rinv2
32  INTEGER n0jet,npart
33  REAL qin(5,*)
34  INTEGER ij(*)
35  REAL frac
36  INTEGER ijout(*)
37  REAL qjet(6,*)
38  REAL e0jade,ejade,durham
39 * ----
40  INTEGER mxjet
41  parameter(mxjet = 30)
42  INTEGER nj(mxjet),iii(500),isortj(500,mxjet),
43  + jetass(500),icutj(mxjet)
44  REAL pjstmp(6,mxjet)
45  INTEGER iseri(500)
46  REAL pj(5,500,mxjet)
47  INTEGER iserij(50,mxjet)
48  REAL q(5,500),qsave(5,500)
49  REAL rsum,rsum0
50  INTEGER nres,iastmp(500)
51  REAL part(6)
52  REAL a,costh,theta
53  INTEGER mxexam
54  parameter(mxexam = 6)
55  REAL exam(mxexam)
56  REAL examin
57  INTEGER iexamn
58 
59 ** INPUTS
60  n=npart
61  DO i=1,n
62  DO j=1,5
63  q(j,i)=qin(j,i)
64  qsave(j,i)=qin(j,i)
65  ENDDO
66  iseri(i)=i
67  ENDDO
68 ** INITIALYZE NUMBER OF PARTICLES IN JET
69  DO i=1,mxjet
70  nj(i)= 0
71  ENDDO
72 ** PARTICLES ARRANGED IN JET
73  DO i=1,n
74 *--- added for cone jet
75  IF(ij(i).GT.0) THEN
76 *--- add end
77  nj(ij(i))=nj(ij(i))+1
78  DO j=1,5
79  pj(j,nj(ij(i)),ij(i)) = q(j,i)
80  ENDDO
81  iserij(nj(ij(i)),ij(i)) = iseri(i)
82  ENDIF
83  ENDDO
84 
85 C WRITE(6,*)NJ(1),NJ(2),NJ(3),NJ(4)
86 ** LEADING PARTICLES SEARCH IN EACH JETS (ORIGINAL NJET) AND RE-DEFINE JETS
87  DO i=1,n0jet
88  DO j=1,6
89  pjstmp(j,i) = 0.
90  ENDDO
91  rsum=0.
92  icutj(i) = 0
93  DO j=1,nj(i)
94  DO k=1,5
95  q(k,j)=pj(k,j,i)
96  ENDDO
97  rsum = rsum + q(4,j)
98  ENDDO
99  CALL syisrt(nj(i),q,iii)
100  DO j=1,nj(i)
101  isortj(j,i) = iii(j)
102  ENDDO
103  rsum0 = 0.
104  DO j=1,nj(i)
105  rsum0 = rsum0 + q(4,iii(j))
106 * --- RE-CALCULATE JET MOMENTUM
107  DO k=1,4
108  pjstmp(k,i) = pjstmp(k,i) + q(k,iii(j))
109  ENDDO
110  pjstmp(6,i)=
111  + sqrt(pjstmp(1,i)**2+pjstmp(2,i)**2+pjstmp(3,i)**2)
112  pjstmp(5,i)=sqrt(max(0.,pjstmp(4,i)**2-pjstmp(6,i)**2))
113 
114  icutj(i) = j
115 * --- CHECK FRACTION OF ENERGY
116  IF(rsum0.GT.rsum*frac) THEN
117  goto 7
118  ENDIF
119  ENDDO
120  7 CONTINUE
121  ENDDO
122 
123 C---re-calculate reference jet momentum
124  DO i=1,n0jet
125  DO j=1,6
126  pjstmp(j,i) = 0.
127  ENDDO
128  DO j=1,nj(i)
129  DO k=1,4
130  pjstmp(k,i) = pjstmp(k,i) + pj(k,j,i)
131  ENDDO
132  ENDDO
133  pjstmp(6,i)=
134  + sqrt(pjstmp(1,i)**2+pjstmp(2,i)**2+pjstmp(3,i)**2)
135  pjstmp(5,i)=sqrt(max(0.,pjstmp(4,i)**2-pjstmp(6,i)**2))
136  ENDDO
137 
138 ** JET ASSIGN CLEAR
139  DO i=1,n
140  jetass(i)=0
141  ENDDO
142 ** JET ASSIGN AGAIN FOR LEADING PARTICLES
143  DO i=1,n0jet
144  DO j=1,icutj(i)
145  ii=isortj(j,i)
146  ii=iserij(ii,i)
147  jetass(ii) = i
148  ENDDO
149  ENDDO
150 ** PUT REMAINING PARTICLES INTO TEMP BUFFER
151  nres = 0
152  DO i=1,n
153  IF(jetass(i).EQ.0) THEN
154  nres = nres + 1
155  iastmp(nres) = iseri(i)
156  DO j=1,5
157  q(j,nres) = qsave(j,i)
158  ENDDO
159  ENDIF
160  ENDDO
161 
162 ** SORT REMAINING PARTICLES IN ENERGY ORDERING
163  IF(nres.EQ.0) THEN
164 C WRITE(*,*)' NRES 0'
165  goto 99
166  ENDIF
167  CALL syisrt(nres,q,iii)
168  99 CONTINUE
169 ** ASSIGNMENT OF THE PARTICLES INTO JETS
170  DO i=1,nres
171  ii=iii(i) ! ID IN RES-BUFFER
172  ii=iastmp(ii) ! SERIAL NUMBER
173  DO j=1,5
174  part(j)=qsave(j,ii)
175  ENDDO
176  part(6)=sqrt(part(1)**2+part(2)**2+part(3)**2)
177 
178  iexamn = 0
179  examin = 1.e+10
180  DO j=1,n0jet
181  a = pjstmp(1,j)*part(1)+
182  + pjstmp(2,j)*part(2)+pjstmp(3,j)*part(3)
183  a = a/part(6)/pjstmp(6,j)
184  a = min(a, 1.)
185  costh = max(a,-1.) ! COSTH
186  theta = acos(costh) ! THERA (ANGLE BETWEEN JET & PART) (1)
187 *--
188  emin2 = min(part(4)**2,pjstmp(4,j)**2)
189  rinv2 = (part(4)+pjstmp(4,j))**2-(part(1)+pjstmp(1,j))**2
190  + -(part(2)+pjstmp(2,j))**2-(part(3)+pjstmp(3,j))**2
191  ejade = max(0.,rinv2)
192  durham= 2.*emin2*max(0.,1.-costh)
193  e0jade= 2.*part(4)*pjstmp(4,j)*max(0.,1.-costh)
194  exam(1) = theta
195  exam(2) = part(4)*pjstmp(4,j)-part(1)*pjstmp(1,j)
196  + -part(2)*pjstmp(2,j)-part(3)*pjstmp(3,j)
197  exam(3) = durham
198  exam(4) = e0jade
199  exam(5) = 8.*part(4)*pjstmp(4,j)*max(0.,(1.-
200  + (part(1)*pjstmp(1,j)+part(2)*pjstmp(2,j)+part(3)*pjstmp(3,j))/
201  + (part(6)*pjstmp(6,j))))/(9.*(part(4)+pjstmp(4,j))**2) ! geneve
202  exam(6) = ejade ! it makes too bad cross talk
203 C EXAM(2) = - SQRT(EMIN2) * COSTH
204 
205  IF(exam(iexam).LT.examin) THEN
206  examin = exam(iexam)
207  iexamn = j
208  ENDIF
209 
210  ENDDO
211 
212 * NEW ASSIGN AND NEW JETS
213  jetass(ii) = iexamn
214 C==== MODE=2 ; Do NOT recalculate Jets (always use original jets for reference)
215 C DO K=1,4
216 C PJSTMP(K,IEXAMN) = PJSTMP(K,IEXAMN) + PART(K)
217 C ENDDO
218 C PJSTMP(6,IEXAMN) = PJSTMP(1,IEXAMN)**2+
219 C + PJSTMP(2,IEXAMN)**2+PJSTMP(3,IEXAMN)**2
220 C PJSTMP(5,IEXAMN) = PJSTMP(4,IEXAMN)**2-PJSTMP(6,IEXAMN)
221 C PJSTMP(6,IEXAMN) = SQRT(MAX(0.,PJSTMP(6,IEXAMN)))
222 C PJSTMP(5,IEXAMN) = SQRT(MAX(0.,PJSTMP(5,IEXAMN)))
223  ENDDO
224 
225  DO i=1,npart
226  ijout(i) = jetass(i)
227  ENDDO
228  DO i=1,n0jet
229  DO j=1,6
230  qjet(j,i)=pjstmp(j,i)
231  ENDDO
232  ENDDO
233 C---re-calculate reference jet momentum
234  DO i=1,n0jet
235  DO j=1,6
236  qjet(j,i) = 0.
237  ENDDO
238  ENDDO
239  DO i=1,npart
240  IF(ijout(i).GT.0.AND.ijout(i).LE.n0jet) THEN
241  DO j=1,4
242  qjet(j,ijout(i)) = qjet(j,ijout(i)) + qin(j,i)
243  ENDDO
244  ENDIF
245  ENDDO
246  DO i=1,n0jet
247  qjet(6,i)=
248  + sqrt(qjet(1,i)**2+qjet(2,i)**2+qjet(3,i)**2)
249  qjet(5,i)=sqrt(max(0.,qjet(4,i)**2-qjet(6,i)**2))
250  ENDDO
251 
252  END
syjetf2
subroutine syjetf2(IEXAM, N0JET, NPART, QIN, IJ, FRAC, IJOUT, QJET)
Definition: syjetf2.f:3
ThrustReconstruction::min
double min(double a, double b)
Definition: ThrustReconstruction.cc:432
syisrt
subroutine syisrt(N, Q, III)
Definition: syisrt.f:3
FPCCDClustering::theta
double theta[MAX_LINK]
Definition: FPCCDClustering.h:240