| 0795afa3 |
1 | #include "isajet/pilot.h" |
| 2 | SUBROUTINE FRGMNT |
| 3 | C |
| 4 | C Control jet fragmentation. Boost to frames defined in |
| 5 | C EVOLVE and call JETGEN. |
| 6 | C |
| 7 | C EVOLVE initializes /PJSET/ as follows-- |
| 8 | C 1 - 2 = PINITS (except for E+E-) |
| 9 | C N0W - N0W = QWJET (for DRELLYAN, NJET=3) |
| 10 | C N0JETS - N0JETS+NJET = PJETS |
| 11 | C N0PAIR - N0PAIR+NPAIR = PPAIR (for WPAIR) |
| 12 | C |
| 13 | #if defined(CERNLIB_IMPNONE) |
| 14 | IMPLICIT NONE |
| 15 | #endif |
| 16 | #include "isajet/itapes.inc" |
| 17 | #include "isajet/primar.inc" |
| 18 | #include "isajet/jetpar.inc" |
| 19 | #include "isajet/pjets.inc" |
| 20 | #include "isajet/pinits.inc" |
| 21 | #include "isajet/partcl.inc" |
| 22 | #include "isajet/const.inc" |
| 23 | #include "isajet/jetset.inc" |
| 24 | #include "isajet/jwork.inc" |
| 25 | #include "isajet/keys.inc" |
| 26 | #include "isajet/q1q2.inc" |
| 27 | #include "isajet/frame.inc" |
| 28 | C |
| 29 | REAL PSUM(5),PALLJ(5),P12(5),PIN(5,2),PWREST(5),PADD(5) |
| 30 | REAL POLD(5),PNEW(5) |
| 31 | REAL PINPL,PINMN,BP,PT2AVE,PTADD,RANF,PHIADD,PALLPL,PALLMN |
| 32 | REAL PALLX,PALLY |
| 33 | INTEGER K,J,JJET,NZERO,IB,NPTCL1,NPTCL2,IFAIL,JET,NPJET1,NPLV1 |
| 34 | INTEGER NPJET3,IP,NP1,NP2,NFIRST,IP1,IFR,NLJ |
| 35 | DOUBLE PRECISION DSUM(5),DPASS(5) |
| 36 | C |
| 37 | C Initialize |
| 38 | DO 100 K=1,5 |
| 39 | 100 DSUM(K)=0. |
| 40 | NLJ=NJET |
| 41 | IF(KEYS(3)) NLJ=NJET+1 |
| 42 | DO 101 J=1,NLJ |
| 43 | JJET=N0JETS+J-1 |
| 44 | IF(JJET.EQ.N0W) GOTO 101 |
| 45 | CALL DBLVEC(PJSET(1,JJET),DPASS) |
| 46 | DO 102 K=1,4 |
| 47 | 102 DSUM(K)=DSUM(K)+DPASS(K) |
| 48 | 101 CONTINUE |
| 49 | DSUM(5)=DSQRT(DSUM(4)**2-DSUM(1)**2-DSUM(2)**2-DSUM(3)**2) |
| 50 | DO 103 K=1,5 |
| 51 | 103 PALLJ(K)=DSUM(K) |
| 52 | C |
| 53 | NZERO=NPTCL+1 |
| 54 | C |
| 55 | C Fragment partons from initial state shower |
| 56 | C |
| 57 | IF(.NOT.KEYS(2)) THEN |
| 58 | DO 110 J=1,NJSET |
| 59 | IF(JDCAY(J).EQ.JPACK*J+J) THEN |
| 60 | IB=JORIG(J)/JPACK-10 |
| 61 | DO 120 K=1,5 |
| 62 | 120 PIN(K,IB)=PJSET(K,J) |
| 63 | ENDIF |
| 64 | 110 CONTINUE |
| 65 | C |
| 66 | CALL FRGJET(11) |
| 67 | CALL FRGJET(12) |
| 68 | C |
| 69 | NPTCL1=NPTCL+1 |
| 70 | NPTCL2=NPTCL1+1 |
| 71 | IF(NPTCL1.GT.MXPTCL) GO TO 9999 |
| 72 | PINPL=.5*(PIN(4,1)+PIN(3,1)+PIN(4,2)+PIN(3,2)) |
| 73 | PINMN=.5*(PIN(4,1)-PIN(3,1)+PIN(4,2)-PIN(3,2)) |
| 74 | PPTCL(1,NPTCL1)=0. |
| 75 | PPTCL(2,NPTCL1)=0. |
| 76 | PPTCL(3,NPTCL1)=HALFE-PINPL |
| 77 | PPTCL(4,NPTCL1)=HALFE-PINPL |
| 78 | PPTCL(5,NPTCL1)=0. |
| 79 | PPTCL(1,NPTCL2)=0. |
| 80 | PPTCL(2,NPTCL2)=0. |
| 81 | PPTCL(3,NPTCL2)=-(HALFE-PINMN) |
| 82 | PPTCL(4,NPTCL2)=HALFE-PINMN |
| 83 | PPTCL(5,NPTCL2)=0. |
| 84 | DO 130 K=1,4 |
| 85 | 130 PSUM(K)=-PALLJ(K) |
| 86 | PSUM(4)=PSUM(4)+ECM |
| 87 | PSUM(5)=PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2 |
| 88 | IF(PSUM(5).GE.0.) THEN |
| 89 | PSUM(5)=SQRT(PSUM(5)) |
| 90 | CALL RESCAL(NZERO,NPTCL2,PSUM,IFAIL) |
| 91 | ENDIF |
| 92 | C |
| 93 | DO 140 K=1,4 |
| 94 | 140 PBEAMS(K)=PPTCL(K,NPTCL1)+PPTCL(K,NPTCL2) |
| 95 | PBEAMS(5)=SQRT(PBEAMS(4)**2-PBEAMS(1)**2-PBEAMS(2)**2 |
| 96 | $ -PBEAMS(3)**2) |
| 97 | ENDIF |
| 98 | C |
| 99 | C Boost partons from final jets with -FRAME |
| 100 | C |
| 101 | 200 DO 210 J=1,NJSET |
| 102 | JET=JORIG(J)/JPACK |
| 103 | IF ( JET.EQ.0 ) THEN |
| 104 | IFR=1 |
| 105 | ELSE |
| 106 | IF(JET.GT.10) GO TO 210 |
| 107 | IF(KEYS(6)) THEN |
| 108 | IF(IDJETS(JET).EQ.10) GO TO 210 |
| 109 | ENDIF |
| 110 | IFR=IFRAME(JET) |
| 111 | ENDIF |
| 112 | C |
| 113 | C Do this boost in double precision for 32-bit machines |
| 114 | CALL DBOOST(-1,FRAME(1,IFR),PJSET(1,J)) |
| 115 | 210 CONTINUE |
| 116 | C |
| 117 | C Fragment partons from final jets |
| 118 | C |
| 119 | NPJET1=NPTCL+1 |
| 120 | DO 220 K=1,4 |
| 121 | 220 PSUM(K)=0 |
| 122 | C |
| 123 | C Conserve mass of 1+2 for DRELLYAN (automatic for WPAIR) |
| 124 | C |
| 125 | IF(KEYS(3)) THEN |
| 126 | CALL FRGJET(1) |
| 127 | CALL FRGJET(2) |
| 128 | IF(STDDY) THEN |
| 129 | DO 230 K=1,4 |
| 130 | PSUM(K)=PJSET(K,3)+PJSET(K,4) |
| 131 | 230 CONTINUE |
| 132 | ELSE |
| 133 | DO 240 K=1,4 |
| 134 | PSUM(K)=PJSET(K,N0W+1)+PJSET(K,N0W+2) |
| 135 | 240 CONTINUE |
| 136 | ENDIF |
| 137 | PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2) |
| 138 | NPLV1=NPTCL |
| 139 | CALL RESCAL(NPJET1,NPLV1,PSUM,IFAIL) |
| 140 | C EXTRADIM has only jet3 + graviton |
| 141 | ELSEIF(KEYS(11)) THEN |
| 142 | CALL FRGJET(3) |
| 143 | CALL FRGJET(0) |
| 144 | NPLV1=NPTCL |
| 145 | DO 241 K=1,4 |
| 146 | PSUM(K)=PJSET(K,3)+PJSET(K,4) |
| 147 | 241 CONTINUE |
| 148 | PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2) |
| 149 | CALL RESCAL (NPJET1,NPLV1,PSUM,IFAIL) |
| 150 | ELSE |
| 151 | C All other processes |
| 152 | DO 242 J=1,NJET |
| 153 | JJET=N0JETS+J-1 |
| 154 | CALL FRGJET(J) |
| 155 | DO 243 K=1,4 |
| 156 | 243 PSUM(K)=PSUM(K)+PJSET(K,JJET) |
| 157 | 242 CONTINUE |
| 158 | PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2) |
| 159 | NPLV1=NPTCL |
| 160 | CALL RESCAL(NPJET1,NPLV1,PSUM,IFAIL) |
| 161 | ENDIF |
| 162 | C |
| 163 | C Add extra jets for DRELLYAN |
| 164 | IF(KEYS(3).AND..NOT.STDDY) THEN |
| 165 | NPJET3=NPTCL+1 |
| 166 | DO 245 J=3,NJET |
| 167 | 245 CALL FRGJET(J) |
| 168 | NPTCL1=NPTCL+1 |
| 169 | IF(NPTCL1.GT.MXPTCL) GO TO 9999 |
| 170 | DO 250 K=1,4 |
| 171 | PPTCL(K,NPTCL1)=PJSET(K,N0W) |
| 172 | 250 PSUM(K)=PJSET(K,N0W) |
| 173 | PPTCL(5,NPTCL1)=PJSET(5,N0W) |
| 174 | DO 246 J=3,NJET |
| 175 | JJET=N0JETS+J-3 |
| 176 | DO 246 K=1,4 |
| 177 | PSUM(K)=PSUM(K)+PJSET(K,JJET) |
| 178 | 246 CONTINUE |
| 179 | PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2) |
| 180 | CALL RESCAL(NPJET3,NPTCL1,PSUM,IFAIL) |
| 181 | DO 260 K=1,5 |
| 182 | 260 PWREST(K)=PPTCL(K,NPTCL1) |
| 183 | ENDIF |
| 184 | C |
| 185 | C Boost partons back to PP COM |
| 186 | C |
| 187 | DO 300 J=1,NJSET |
| 188 | JET=JORIG(J)/JPACK |
| 189 | IF ( JET.EQ.0 ) THEN |
| 190 | IFR=1 |
| 191 | ELSE |
| 192 | IF(JET.GT.10) GO TO 300 |
| 193 | IF(KEYS(6)) THEN |
| 194 | IF(IDJETS(JET).EQ.10) GO TO 300 |
| 195 | ENDIF |
| 196 | IFR=IFRAME(JET) |
| 197 | ENDIF |
| 198 | BP=0. |
| 199 | DO 305 K=1,3 |
| 200 | 305 BP=BP+FRAME(K,IFR)*PJSET(K,J) |
| 201 | BP=BP/FRAME(5,IFR) |
| 202 | DO 310 K=1,3 |
| 203 | 310 PJSET(K,J)=PJSET(K,J)+FRAME(K,IFR)*PJSET(4,J)/FRAME(5,IFR) |
| 204 | $ +FRAME(K,IFR)*BP/(FRAME(4,IFR)+FRAME(5,IFR)) |
| 205 | PJSET(4,J)=FRAME(4,IFR)*PJSET(4,J)/FRAME(5,IFR)+BP |
| 206 | 300 CONTINUE |
| 207 | C |
| 208 | C Reset FRAME to boost hadrons to PP COM |
| 209 | C |
| 210 | IF(KEYS(1).OR.KEYS(2).OR.(KEYS(3).AND.NJET.EQ.2).OR.KEYS(5) |
| 211 | $.OR.(KEYS(7).AND.NPAIR.EQ.0).OR.KEYS(8)) THEN |
| 212 | DO 410 K=1,5 |
| 213 | FRAME(K,1)=PALLJ(K) |
| 214 | 410 CONTINUE |
| 215 | ELSEIF(KEYS(3).AND.NJET.GT.2) THEN |
| 216 | DO 420 K=1,5 |
| 217 | 420 FRAME(K,1)=PALLJ(K) |
| 218 | BP=0. |
| 219 | DO 430 K=1,3 |
| 220 | 430 BP=BP+FRAME(K,1)*PWREST(K) |
| 221 | BP=BP/FRAME(5,1) |
| 222 | DO 440 K=1,3 |
| 223 | FRAME(K,2)=PWREST(K)+FRAME(K,1)*PWREST(4)/FRAME(5,1) |
| 224 | $ +FRAME(K,1)*BP/(FRAME(4,1)+FRAME(5,1)) |
| 225 | 440 CONTINUE |
| 226 | FRAME(4,2)=FRAME(4,1)*PWREST(4)/FRAME(5,1)+BP |
| 227 | ENDIF |
| 228 | C |
| 229 | C Boost hadrons back to PP COM |
| 230 | C |
| 231 | DO 500 IP=NZERO,NPTCL |
| 232 | JET=IABS(IORIG(IP))/IPACK |
| 233 | IF(JET.GT.10) GO TO 500 |
| 234 | IF(KEYS(6)) THEN |
| 235 | IF(IDJETS(JET).EQ.10) GO TO 500 |
| 236 | ENDIF |
| 237 | IF(JET.EQ.0) THEN |
| 238 | IFR=1 |
| 239 | ELSE |
| 240 | IFR=IFRAME(JET) |
| 241 | ENDIF |
| 242 | BP=0. |
| 243 | DO 510 K=1,3 |
| 244 | 510 BP=BP+FRAME(K,IFR)*PPTCL(K,IP) |
| 245 | BP=BP/FRAME(5,IFR) |
| 246 | DO 520 K=1,3 |
| 247 | 520 PPTCL(K,IP)=PPTCL(K,IP)+FRAME(K,IFR)*PPTCL(4,IP)/FRAME(5,IFR) |
| 248 | $ +FRAME(K,IFR)*BP/(FRAME(4,IFR)+FRAME(5,IFR)) |
| 249 | PPTCL(4,IP)=FRAME(4,IFR)*PPTCL(4,IP)/FRAME(5,IFR)+BP |
| 250 | 500 CONTINUE |
| 251 | C |
| 252 | C Add intrinsic PT |
| 253 | C |
| 254 | IF(.NOT.KEYS(2)) THEN |
| 255 | PT2AVE=.1*SQRT(QSQ) |
| 256 | PTADD=SQRT(-PT2AVE*ALOG(RANF())) |
| 257 | PHIADD=2.*PI*RANF() |
| 258 | PADD(1)=2.*PTADD*COS(PHIADD) |
| 259 | PADD(2)=2.*PTADD*SIN(PHIADD) |
| 260 | C Must use large and small components carefully to calculate |
| 261 | C mass on 32-bit machines. |
| 262 | PALLPL=0. |
| 263 | PALLMN=0. |
| 264 | PALLX=0. |
| 265 | PALLY=0. |
| 266 | DO 525 IP=NZERO,NPTCL |
| 267 | PALLX=PALLX+PPTCL(1,IP) |
| 268 | PALLY=PALLY+PPTCL(2,IP) |
| 269 | IF(PPTCL(3,IP).GT.0.) THEN |
| 270 | PALLPL=PALLPL+(PPTCL(4,IP)+PPTCL(3,IP)) |
| 271 | PALLMN=PALLMN+(PPTCL(1,IP)**2+PPTCL(2,IP)**2+PPTCL(5,IP)**2) |
| 272 | $ /(PPTCL(4,IP)+PPTCL(3,IP)) |
| 273 | ELSE |
| 274 | PALLMN=PALLMN+(PPTCL(4,IP)-PPTCL(3,IP)) |
| 275 | PALLPL=PALLPL+(PPTCL(1,IP)**2+PPTCL(2,IP)**2+PPTCL(5,IP)**2) |
| 276 | $ /(PPTCL(4,IP)-PPTCL(3,IP)) |
| 277 | ENDIF |
| 278 | 525 CONTINUE |
| 279 | POLD(1)=PALLX |
| 280 | POLD(2)=PALLY |
| 281 | POLD(3)=.5*(PALLPL-PALLMN) |
| 282 | POLD(4)=.5*(PALLPL+PALLMN) |
| 283 | POLD(5)=SQRT(PALLPL*PALLMN-PALLX**2-PALLY**2) |
| 284 | PNEW(1)=PADD(1)+POLD(1) |
| 285 | PNEW(2)=PADD(2)+POLD(2) |
| 286 | PNEW(3)=POLD(3) |
| 287 | PNEW(4)=SQRT(PNEW(1)**2+PNEW(2)**2+PNEW(3)**2+POLD(5)**2) |
| 288 | PNEW(5)=POLD(5) |
| 289 | C |
| 290 | DO 530 IP=NZERO,NPTCL |
| 291 | BP=0. |
| 292 | DO 531 K=1,3 |
| 293 | 531 BP=BP+POLD(K)*PPTCL(K,IP) |
| 294 | BP=BP/POLD(5) |
| 295 | DO 532 K=1,3 |
| 296 | 532 PPTCL(K,IP)=PPTCL(K,IP)-POLD(K)*PPTCL(4,IP)/POLD(5) |
| 297 | $ +POLD(K)*BP/(POLD(4)+POLD(5)) |
| 298 | PPTCL(4,IP)=PPTCL(4,IP)*POLD(4)/POLD(5)-BP |
| 299 | C |
| 300 | BP=0. |
| 301 | DO 533 K=1,3 |
| 302 | 533 BP=BP+PNEW(K)*PPTCL(K,IP) |
| 303 | BP=BP/PNEW(5) |
| 304 | DO 534 K=1,3 |
| 305 | 534 PPTCL(K,IP)=PPTCL(K,IP)+PNEW(K)*PPTCL(4,IP)/PNEW(5) |
| 306 | $ +PNEW(K)*BP/(PNEW(4)+PNEW(5)) |
| 307 | PPTCL(4,IP)=PPTCL(4,IP)*PNEW(4)/PNEW(5)+BP |
| 308 | 530 CONTINUE |
| 309 | C |
| 310 | C Add opposite PT to beam jets |
| 311 | DO 541 K=1,4 |
| 312 | 541 PBEAMS(K)=-PNEW(K) |
| 313 | PBEAMS(4)=PBEAMS(4)+ECM |
| 314 | PBEAMS(5)=PBEAMS(4)**2-PBEAMS(1)**2-PBEAMS(2)**2 -PBEAMS(3)**2 |
| 315 | IF ( PBEAMS(5).GT.0 ) THEN |
| 316 | PBEAMS(5)=SQRT(PBEAMS(5)) |
| 317 | ELSE |
| 318 | PBEAMS(4)=SQRT(PBEAMS(4)**2-PBEAMS(5)+4.) |
| 319 | PBEAMS(5)=2. |
| 320 | ENDIF |
| 321 | ENDIF |
| 322 | C |
| 323 | C Decay hadrons |
| 324 | C |
| 325 | NP1=NZERO |
| 326 | 600 NP2=NPTCL |
| 327 | DO 610 IP=NP1,NP2 |
| 328 | NFIRST=NPTCL+1 |
| 329 | JET=IABS(IORIG(IP))/IPACK |
| 330 | CALL DECAY(IP) |
| 331 | DO 620 IP1=NFIRST,NPTCL |
| 332 | 620 IORIG(IP1)=ISIGN(IABS(IORIG(IP1))+IPACK*JET,IORIG(IP1)) |
| 333 | 610 CONTINUE |
| 334 | NP1=NP2+1 |
| 335 | IF(NP1.LE.NPTCL) GO TO 600 |
| 336 | RETURN |
| 337 | C |
| 338 | C Error |
| 339 | C |
| 340 | 9999 CALL PRTEVT(0) |
| 341 | WRITE(ITLIS,9998) NPTCL |
| 342 | 9998 FORMAT(//' ERROR IN FRGMNT ... NPTCL > ',I6) |
| 343 | RETURN |
| 344 | END |
git://git.uio.no / u / mrichter / AliRoot.git / blame