Updating info for ACORDE and TRD
[u/mrichter/AliRoot.git] / PYTHIA6 / pythia6214.f
index e8e6c35..28e5d1e 100644 (file)
@@ -286,7 +286,7 @@ C...and particle, decay and process data.
 C...Double precision and integer declarations.
       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
       IMPLICIT INTEGER(I-N)
-      INTEGER PYK,PYCHGE,PYCOMP
+C      INTEGER PYK,PYCHGE,PYCOMP
 C...Commonblocks.
       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
@@ -2563,8 +2563,8 @@ C...Local arrays and character variables.
       CHARACTER CHFRAM*12,CHBEAM*12,CHTARG*12,CHLH(2)*6
  
 C...Interface to PDFLIB.
-      COMMON/W50512/QCDL4,QCDL5
-      SAVE /W50512/
+      COMMON/LW50512/QCDL4,QCDL5
+      SAVE /LW50512/
       DOUBLE PRECISION VALUE(20),QCDL4,QCDL5
       CHARACTER*20 PARM(20)
       DATA VALUE/20*0D0/,PARM/20*' '/
@@ -3979,32 +3979,37 @@ C...Find mass and evaluate width.
         MINT(51)=0
  
 C...Evaluate suppression factors due to non-simulated channels.
-        IF(KCHG(KC,3).EQ.0) THEN
-          WIDS(KC,1)=((WDTE(0,1)+WDTE(0,2))**2+
-     &    2D0*(WDTE(0,1)+WDTE(0,2))*(WDTE(0,4)+WDTE(0,5))+
-     &    2D0*WDTE(0,4)*WDTE(0,5))/WDTP(0)**2
-          WIDS(KC,2)=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/WDTP(0)
-          WIDS(KC,3)=0D0
-          WIDS(KC,4)=0D0
-          WIDS(KC,5)=0D0
-        ELSE
-          IF(MWID(KC).EQ.3) MINT(63)=1
-          CALL PYWIDT(-KF,PMR**2,WDTPM,WDTEM)
-          MINT(51)=0
-          WIDS(KC,1)=((WDTE(0,1)+WDTE(0,2))*(WDTEM(0,1)+WDTEM(0,3))+
-     &    (WDTE(0,1)+WDTE(0,2))*(WDTEM(0,4)+WDTEM(0,5))+
-     &    (WDTE(0,4)+WDTE(0,5))*(WDTEM(0,1)+WDTEM(0,3))+
-     &    WDTE(0,4)*WDTEM(0,5)+WDTE(0,5)*WDTEM(0,4))/WDTP(0)**2
-          WIDS(KC,2)=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/WDTP(0)
-          WIDS(KC,3)=(WDTEM(0,1)+WDTEM(0,3)+WDTEM(0,4))/WDTP(0)
-          WIDS(KC,4)=((WDTE(0,1)+WDTE(0,2))**2+
-     &    2D0*(WDTE(0,1)+WDTE(0,2))*(WDTE(0,4)+WDTE(0,5))+
-     &    2D0*WDTE(0,4)*WDTE(0,5))/WDTP(0)**2
-          WIDS(KC,5)=((WDTEM(0,1)+WDTEM(0,3))**2+
-     &    2D0*(WDTEM(0,1)+WDTEM(0,3))*(WDTEM(0,4)+WDTEM(0,5))+
-     &    2D0*WDTEM(0,4)*WDTEM(0,5))/WDTP(0)**2
+C...AM
+C...Protection against division by 0 since rho_21_tc is causing problem here
+        IF (WDTP(0) .GT. 0.) THEN
+           
+           IF(KCHG(KC,3).EQ.0) THEN
+              WIDS(KC,1)=((WDTE(0,1)+WDTE(0,2))**2+
+     &             2D0*(WDTE(0,1)+WDTE(0,2))*(WDTE(0,4)+WDTE(0,5))+
+     &             2D0*WDTE(0,4)*WDTE(0,5))/WDTP(0)**2
+              WIDS(KC,2)=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/WDTP(0)
+              WIDS(KC,3)=0D0
+              WIDS(KC,4)=0D0
+              WIDS(KC,5)=0D0
+           ELSE
+              IF(MWID(KC).EQ.3) MINT(63)=1
+              CALL PYWIDT(-KF,PMR**2,WDTPM,WDTEM)
+              MINT(51)=0
+              WIDS(KC,1)=((WDTE(0,1)+WDTE(0,2))*(WDTEM(0,1)+WDTEM(0,3))+
+     &             (WDTE(0,1)+WDTE(0,2))*(WDTEM(0,4)+WDTEM(0,5))+
+     &             (WDTE(0,4)+WDTE(0,5))*(WDTEM(0,1)+WDTEM(0,3))+
+     &             WDTE(0,4)*WDTEM(0,5)+WDTE(0,5)*WDTEM(0,4))/WDTP(0)**2
+              WIDS(KC,2)=(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/WDTP(0)
+              WIDS(KC,3)=(WDTEM(0,1)+WDTEM(0,3)+WDTEM(0,4))/WDTP(0)
+              WIDS(KC,4)=((WDTE(0,1)+WDTE(0,2))**2+
+     &             2D0*(WDTE(0,1)+WDTE(0,2))*(WDTE(0,4)+WDTE(0,5))+
+     &             2D0*WDTE(0,4)*WDTE(0,5))/WDTP(0)**2
+              WIDS(KC,5)=((WDTEM(0,1)+WDTEM(0,3))**2+
+     &             2D0*(WDTEM(0,1)+WDTEM(0,3))*(WDTEM(0,4)+WDTEM(0,5))+
+     &             2D0*WDTEM(0,4)*WDTEM(0,5))/WDTP(0)**2
+           ENDIF
+           
         ENDIF
 C...Set resonance widths and branching ratios;
 C...also on/off switch for decays.
         IF(MWID(KC).EQ.1.OR.MWID(KC).EQ.3) THEN
@@ -27687,8 +27692,8 @@ C...Local arrays.
      &XPPI(-6:6),XPPR(-6:6)
  
 C...Interface to PDFLIB.
-      COMMON/W50513/XMIN,XMAX,Q2MIN,Q2MAX
-      SAVE /W50513/
+      COMMON/LW50513/XMIN,XMAX,Q2MIN,Q2MAX
+      SAVE /LW50513/
       DOUBLE PRECISION XX,QQ,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU,
      &VALUE(20),XMIN,XMAX,Q2MIN,Q2MAX
       CHARACTER*20 PARM(20)
@@ -28353,8 +28358,8 @@ C...Local arrays.
       DIMENSION XPEL(-25:25),XPGA(-6:6),SXP(0:6)
  
 C...Interface to PDFLIB.
-      COMMON/W50513/XMIN,XMAX,Q2MIN,Q2MAX
-      SAVE /W50513/
+      COMMON/LW50513/XMIN,XMAX,Q2MIN,Q2MAX
+      SAVE /LW50513/
       DOUBLE PRECISION XX,QQ,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU,
      &VALUE(20),XMIN,XMAX,Q2MIN,Q2MAX
       CHARACTER*20 PARM(20)
@@ -45407,7 +45412,6 @@ C...Commonblocks.
       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYINT4/
 C...Local array.
       DIMENSION PS(2,6),IJOIN(100)
 C...Initialize and reset.
       MSTU(24)=0
       IF(MSTU(12).GE.1) CALL PYLIST(0)
@@ -50094,6 +50098,41 @@ C.. Form baryon. Distinguish Lambda- and Sigmalike baryons.
         KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+KFLS,KFL1)
         IF(KSIG.EQ.0) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+KFLS,KFL1)
       ENDIF
+C -------------------------------------------------------------------------
+C Extracted from a private e-mail exchange with Torbjorn Sjostrand
+C 
+C No, Lambda(1520) is not included and not foreseen.
+C So if you want it in Pythia, it would have to be a hack.
+C What you could do is:
+C 1) In PYKFDI, just before the RETURN above label 140, you could check if
+C a Lambda, Sigma0 or Sigma*0 has been produced, and with some small
+C probability switch such a particle to the Lambda(1520) code. That is,
+C if KF = 3122, 3212, or 3214 and a random number below some number, switch
+C to KF = 3124. (And correspondingly for anticparticles.)
+C 2) Use the PYUPDA routine (see manual) to include particle and decay data
+C for the Lambda(1520).
+C -------------------------------------------------------------------------
+      IF (IABS(KF).EQ.3122) THEN
+C     Converting a fraction (0.20) of Lambda0 to Lambda(1520) + c.c.
+C     This fraction is based on the experimental measurement at ISR
+C     Bobbink 83, NP B217,11 (1983)
+C     The region 0.5 < XF < 1.0 has been extrapolated to XF=0
+         IF(PYR(0).LE.0.20) KF=ISIGN(3124,KF)
+      ENDIF
+
+      IF(IABS(KF).EQ.3212) THEN
+C     Converting a fraction (0.20) of Sigma0 to Lambda(1520) + c.c.
+C     We suppose the same fraction as for Lambda0
+         IF(PYR(0).LE.0.20) KF=ISIGN(3124,KF)
+      ENDIF
+
+      IF (IABS(KF).EQ.3214) THEN
+C     Converting a fraction (0.30) of Sigma0(1385) to Lambda(1520) + c.c.
+C     This is conservative extimate supposing that the ratio
+C     scales as (M_Sigma1385/M_Lambda0)^2 ~ 1.5 
+         IF(PYR(0).LE.0.30) KF=ISIGN(3124,KF)
+      ENDIF
       RETURN
  
 C...Use tabulated probabilities to select new flavour and hadron.
@@ -51375,10 +51414,13 @@ C...Calculate allowed z range.
  
 C...Integral of Altarelli-Parisi z kernel for QCD.
 C...(Includes squark and gluino; with factor N_C/C_F extra for latter).
+      FMED = PARJ(200)
       IF(MSTJ(49).EQ.0.AND.KFL(1).EQ.21) THEN
-        FBR=6D0*LOG((1D0-ZC)/ZC)+MSTJ(45)*0.5D0
+C Nestor
+        FBR=(1.D0+FMED)*6D0*LOG((1D0-ZC)/ZC)+MSTJ(45)*0.5D0
       ELSEIF(MSTJ(49).EQ.0) THEN
-        FBR=(8D0/3D0)*LOG((1D0-ZC)/ZC)
+C Nestor
+        FBR=(1.D0+FMED)*(8D0/3D0)*LOG((1D0-ZC)/ZC)
         IF(IGLUI.EQ.1.AND.IR.GE.31) FBR=FBR*(9D0/4D0)
  
 C...Integral of Altarelli-Parisi z kernel for scalar gluon.
@@ -51390,9 +51432,9 @@ C...Integral of Altarelli-Parisi z kernel for scalar gluon.
  
 C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon.
       ELSEIF(KFL(1).EQ.21) THEN
-        FBR=6D0*MSTJ(45)*(0.5D0-ZC)
+        FBR=(1.D0+FMED)*6D0*MSTJ(45)*(0.5D0-ZC)
       ELSE
-        FBR=2D0*LOG((1D0-ZC)/ZC)
+        FBR=(1.D0+FMED)*2D0*LOG((1D0-ZC)/ZC)
       ENDIF
  
 C...Reset QCD probability for colourless.
@@ -54958,24 +55000,26 @@ C...Identify particle code and whether already defined  (for MUPDA=3).
   150         CONTINUE
             ENDIF
 C...Remove duplicate old decay data.
-            IF(KCREP.NE.0.AND.MDCY(KCREP,3).GT.0) THEN
-              IDCREP=MDCY(KCREP,2)
-              NDCREP=MDCY(KCREP,3)
-              DO 160 I=1,KCC
-                IF(MDCY(I,2).GT.IDCREP) MDCY(I,2)=MDCY(I,2)-NDCREP
-  160         CONTINUE
-              DO 180 I=IDCREP,NDC-NDCREP
-                MDME(I,1)=MDME(I+NDCREP,1)
-                MDME(I,2)=MDME(I+NDCREP,2)
-                BRAT(I)=BRAT(I+NDCREP)
-                DO 170 J=1,5
-                  KFDP(I,J)=KFDP(I+NDCREP,J)
-  170           CONTINUE
-  180         CONTINUE
-              NDC=NDC-NDCREP
-              KC=KCREP
-            ELSEIF(KCREP.NE.0) THEN
-              KC=KCREP
+            IF(KCREP.NE.0) THEN
+               IF(MDCY(KCREP,3).GT.0) THEN
+                  IDCREP=MDCY(KCREP,2)
+                  NDCREP=MDCY(KCREP,3)
+                  DO 160 I=1,KCC
+                     IF(MDCY(I,2).GT.IDCREP) MDCY(I,2)=MDCY(I,2)-NDCREP
+ 160              CONTINUE
+                  DO 180 I=IDCREP,NDC-NDCREP
+                     MDME(I,1)=MDME(I+NDCREP,1)
+                     MDME(I,2)=MDME(I+NDCREP,2)
+                     BRAT(I)=BRAT(I+NDCREP)
+                     DO 170 J=1,5
+                        KFDP(I,J)=KFDP(I+NDCREP,J)
+ 170                 CONTINUE
+ 180              CONTINUE
+                  NDC=NDC-NDCREP
+                  KC=KCREP
+               ELSE
+                  KC=KCREP
+               ENDIF
             ELSE
               KCC=KCC+1
               KC=KCC
@@ -56194,6 +56238,7 @@ C...Loop over all particles. Find cell that was hit by given particle.
      &    KC.EQ.18) GOTO 110
           IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.PYCHGE(K(I,2)).EQ.0)
      &    GOTO 110
+
         ENDIF
         NP=NP+1
         PT=SQRT(P(I,1)**2+P(I,2)**2)
@@ -59654,14 +59699,14 @@ C      IDATI(5)=IMIN
 C      IDATI(6)=ISEC
  
 C...Example 4: GNU LINUX libU77, SunOS.
-      CALL IDATE(IDTEMP)
-      IDATI(1)=IDTEMP(3)
-      IDATI(2)=IDTEMP(2)
-      IDATI(3)=IDTEMP(1)
-      CALL ITIME(IDTEMP)
-      IDATI(4)=IDTEMP(1)
-      IDATI(5)=IDTEMP(2)
-      IDATI(6)=IDTEMP(3)
+c      CALL IDATE(IDTEMP)
+c      IDATI(1)=IDTEMP(3)
+c      IDATI(2)=IDTEMP(2)
+c      IDATI(3)=IDTEMP(1)
+c      CALL ITIME(IDTEMP)
+c      IDATI(4)=IDTEMP(1)
+c      IDATI(5)=IDTEMP(2)
+c      IDATI(6)=IDTEMP(3)
  
 C...Common code to ensure right century.
       IDATI(1)=2000+MOD(IDATI(1),100)