All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
syjetf1.f
Go to the documentation of this file.
1 CDECK ID>, SYJETF1.
2 ************************************************************************
3  SUBROUTINE syjetf1(IEXAM,N0JET,NPART,QIN,IJ,FRAC,IJOUT,QJET)
4 * Particle Re-association Routine (Mode=1)
5 * Iteratively, core is re-defined when a particle is added to the 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 ** JET ASSIGN CLEAR
123  DO i=1,n
124  jetass(i)=0
125  ENDDO
126 ** JET ASSIGN AGAIN FOR LEADING PARTICLES
127  DO i=1,n0jet
128  DO j=1,icutj(i)
129  ii=isortj(j,i)
130  ii=iserij(ii,i)
131  jetass(ii) = i
132  ENDDO
133  ENDDO
134 ** PUT REMAINING PARTICLES INTO TEMP BUFFER
135  nres = 0
136  DO i=1,n
137  IF(jetass(i).EQ.0) THEN
138  nres = nres + 1
139  iastmp(nres) = iseri(i)
140  DO j=1,5
141  q(j,nres) = qsave(j,i)
142  ENDDO
143  ENDIF
144  ENDDO
145 
146 ** SORT REMAINING PARTICLES IN ENERGY ORDERING
147  IF(nres.EQ.0) THEN
148 C WRITE(*,*)' NRES 0'
149  goto 99
150  ENDIF
151  CALL syisrt(nres,q,iii)
152  99 CONTINUE
153 ** ASSIGNMENT OF THE PARTICLES INTO JETS
154  DO i=1,nres
155  ii=iii(i) ! ID IN RES-BUFFER
156  ii=iastmp(ii) ! SERIAL NUMBER
157  DO j=1,5
158  part(j)=qsave(j,ii)
159  ENDDO
160  part(6)=sqrt(part(1)**2+part(2)**2+part(3)**2)
161 
162  iexamn = 0
163  examin = 1.e+10
164  DO j=1,n0jet
165  a = pjstmp(1,j)*part(1)+
166  + pjstmp(2,j)*part(2)+pjstmp(3,j)*part(3)
167  a = a/part(6)/pjstmp(6,j)
168  a = min(a, 1.)
169  costh = max(a,-1.) ! COSTH
170  theta = acos(costh) ! THERA (ANGLE BETWEEN JET & PART) (1)
171 *--
172  emin2 = min(part(4)**2,pjstmp(4,j)**2)
173  rinv2 = (part(4)+pjstmp(4,j))**2-(part(1)+pjstmp(1,j))**2
174  + -(part(2)+pjstmp(2,j))**2-(part(3)+pjstmp(3,j))**2
175  ejade = max(0.,rinv2)
176  durham= 2.*emin2*max(0.,1.-costh)
177  e0jade= 2.*part(4)*pjstmp(4,j)*max(0.,1.-costh)
178 
179  exam(1) = theta
180  exam(2) = part(4)*pjstmp(4,j)-part(1)*pjstmp(1,j)
181  + -part(2)*pjstmp(2,j)-part(3)*pjstmp(3,j)
182  exam(3) = durham
183  exam(4) = e0jade
184  exam(5) = 8.*part(4)*pjstmp(4,j)*max(0.,(1.-
185  + (part(1)*pjstmp(1,j)+part(2)*pjstmp(2,j)+part(3)*pjstmp(3,j))/
186  + (part(6)*pjstmp(6,j))))/(9.*(part(4)+pjstmp(4,j))**2) ! Geneva
187  exam(6) = ejade ! it makes too bad cross talk
188 C EXAM(2) = - SQRT(EMIN2) * COSTH
189 
190  IF(exam(iexam).LT.examin) THEN
191  examin = exam(iexam)
192  iexamn = j
193  ENDIF
194 
195  ENDDO
196 
197 * NEW ASSIGN AND NEW JETS
198  jetass(ii) = iexamn
199  DO k=1,4
200  pjstmp(k,iexamn) = pjstmp(k,iexamn) + part(k)
201  ENDDO
202  pjstmp(6,iexamn) = pjstmp(1,iexamn)**2+
203  + pjstmp(2,iexamn)**2+pjstmp(3,iexamn)**2
204  pjstmp(5,iexamn) = pjstmp(4,iexamn)**2-pjstmp(6,iexamn)
205  pjstmp(6,iexamn) = sqrt(max(0.,pjstmp(6,iexamn)))
206  pjstmp(5,iexamn) = sqrt(max(0.,pjstmp(5,iexamn)))
207  ENDDO
208 
209  DO i=1,npart
210  ijout(i) = jetass(i)
211  ENDDO
212  DO i=1,n0jet
213  DO j=1,6
214  qjet(j,i)=pjstmp(j,i)
215  ENDDO
216  ENDDO
217  END
syjetf1
subroutine syjetf1(IEXAM, N0JET, NPART, QIN, IJ, FRAC, IJOUT, QJET)
Definition: syjetf1.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