ISAJET/code/rescal.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE RESCAL(N1,N2,PSUM,IFAIL)
   3 C          RESCALE MOMENTA OF PARTICLES N1...N2 TO GIVE TOTAL
   4 C          FOUR-MOMENTUM PSUM.
   5 C          RETURN IFAIL=0 IF OK, IFAIL=1 IF NO GOOD.
   6 #include "isajet/itapes.inc"
   7 #include "isajet/partcl.inc"
   8       DIMENSION PSUM(5),PADD(5),BETA(3)
   9       DATA ERRLIM/.0001/
  10 C          ORIGIONAL MOMENTUM IS PADD.
  11       IFAIL=1
  12       IF(N1.GE.N2) RETURN
  13       DO 100 K=1,5
  14 100   PADD(K)=0.
  15       DO 110 IP=N1,N2
  16       DO 110 K=1,5
  17       PADD(K)=PADD(K)+PPTCL(K,IP)
  18 110   CONTINUE
  19       IF(PADD(5).GE.PSUM(5)) RETURN
  20       PADD(5)=PADD(4)**2-PADD(1)**2-PADD(2)**2-PADD(3)**2
  21       IF(PADD(5).LE.0) RETURN
  22       PADD(5)=SQRT(PADD(5))
  23       DO 120 K=1,3
  24 120   BETA(K)=-PADD(K)/PADD(5)
  25       GAMMA=PADD(4)/PADD(5)
  26 C          BOOST PARTICLES TO REST.
  27 200   CONTINUE
  28       DO 210 IP=N1,N2
  29       BP=0.
  30       DO 220 K=1,3
  31 220   BP=BP+PPTCL(K,IP)*BETA(K)
  32       DO 230 K=1,3
  33 230   PPTCL(K,IP)=PPTCL(K,IP)+BETA(K)*PPTCL(4,IP)
  34      $+BETA(K)*BP/(GAMMA+1.)
  35       PPTCL(4,IP)=GAMMA*PPTCL(4,IP)+BP
  36 210   CONTINUE
  37       IF(IFAIL.EQ.0) RETURN
  38 C          RESCALE MOMENTA IN REST FRAME.
  39       SCAL=1.
  40       DO 301 IPASS=1,200
  41       SUM=0.
  42       DO 310 IP=N1,N2
  43       DO 320 K=1,3
  44 320   PPTCL(K,IP)=SCAL*PPTCL(K,IP)
  45       PPTCL(4,IP)=SQRT(PPTCL(1,IP)**2+PPTCL(2,IP)**2+PPTCL(3,IP)**2
  46      $+PPTCL(5,IP)**2)
  47       SUM=SUM+PPTCL(4,IP)
  48 310   CONTINUE
  49       SCAL=PSUM(5)/SUM
  50 301   IF(ABS(SCAL-1.).LE.ERRLIM) GO TO 300
  51 300   CONTINUE
  52 C          BOOST BACK WITH PSUM.
  53       BMAG=0.
  54       DO 400 K=1,3
  55       BETA(K)=PSUM(K)/PSUM(5)
  56       BMAG=BMAG+ABS(BETA(K))
  57 400   CONTINUE
  58       GAMMA=PSUM(4)/PSUM(5)
  59       IFAIL=0
  60       IF(BMAG.EQ.0.) RETURN
  61       GO TO 200
  62       END