]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA/jetset/lucell.F
This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / PYTHIA / jetset / lucell.F
diff --git a/PYTHIA/jetset/lucell.F b/PYTHIA/jetset/lucell.F
new file mode 100644 (file)
index 0000000..fdb5b81
--- /dev/null
@@ -0,0 +1,205 @@
+C********************************************************************* 
+      SUBROUTINE LUCELL(NJET) 
+C...Purpose: to provide a simple way of jet finding in an eta-phi-ET 
+C...coordinate frame, as used for calorimeters at hadron colliders. 
+      COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) 
+      COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
+      COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) 
+      SAVE /LUJETS/,/LUDAT1/,/LUDAT2/ 
+C...Loop over all particles. Find cell that was hit by given particle. 
+      PTLRAT=1./SINH(PARU(51))**2 
+      NP=0 
+      NC=N 
+      DO 110 I=1,N 
+      IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 110 
+      IF(P(I,1)**2+P(I,2)**2.LE.PTLRAT*P(I,3)**2) GOTO 110 
+      IF(MSTU(41).GE.2) THEN 
+        KC=LUCOMP(K(I,2)) 
+        IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR. 
+     &  KC.EQ.18) GOTO 110 
+        IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0) 
+     &  GOTO 110 
+      ENDIF 
+      NP=NP+1 
+      PT=SQRT(P(I,1)**2+P(I,2)**2) 
+      ETA=SIGN(LOG((SQRT(PT**2+P(I,3)**2)+ABS(P(I,3)))/PT),P(I,3)) 
+      IETA=MAX(1,MIN(MSTU(51),1+INT(MSTU(51)*0.5*(ETA/PARU(51)+1.)))) 
+      PHI=ULANGL(P(I,1),P(I,2)) 
+      IPHI=MAX(1,MIN(MSTU(52),1+INT(MSTU(52)*0.5*(PHI/PARU(1)+1.)))) 
+      IETPH=MSTU(52)*IETA+IPHI 
+C...Add to cell already hit, or book new cell. 
+      DO 100 IC=N+1,NC 
+      IF(IETPH.EQ.K(IC,3)) THEN 
+        K(IC,4)=K(IC,4)+1 
+        P(IC,5)=P(IC,5)+PT 
+        GOTO 110 
+      ENDIF 
+  100 CONTINUE 
+      IF(NC.GE.MSTU(4)-MSTU(32)-5) THEN 
+        CALL LUERRM(11,'(LUCELL:) no more memory left in LUJETS') 
+        NJET=-2 
+        RETURN 
+      ENDIF 
+      NC=NC+1 
+      K(NC,3)=IETPH 
+      K(NC,4)=1 
+      K(NC,5)=2 
+      P(NC,1)=(PARU(51)/MSTU(51))*(2*IETA-1-MSTU(51)) 
+      P(NC,2)=(PARU(1)/MSTU(52))*(2*IPHI-1-MSTU(52)) 
+      P(NC,5)=PT 
+  110 CONTINUE 
+C...Smear true bin content by calorimeter resolution. 
+      IF(MSTU(53).GE.1) THEN 
+        DO 130 IC=N+1,NC 
+        PEI=P(IC,5) 
+        IF(MSTU(53).EQ.2) PEI=P(IC,5)*COSH(P(IC,1)) 
+  120   PEF=PEI+PARU(55)*SQRT(-2.*LOG(MAX(1E-10,RLU(0)))*PEI)* 
+     &  COS(PARU(2)*RLU(0)) 
+        IF(PEF.LT.0..OR.PEF.GT.PARU(56)*PEI) GOTO 120 
+        P(IC,5)=PEF 
+        IF(MSTU(53).EQ.2) P(IC,5)=PEF/COSH(P(IC,1)) 
+  130   CONTINUE 
+      ENDIF 
+C...Remove cells below threshold. 
+      IF(PARU(58).GT.0.) THEN 
+        NCC=NC 
+        NC=N 
+        DO 140 IC=N+1,NCC 
+        IF(P(IC,5).GT.PARU(58)) THEN 
+          NC=NC+1 
+          K(NC,3)=K(IC,3) 
+          K(NC,4)=K(IC,4) 
+          K(NC,5)=K(IC,5) 
+          P(NC,1)=P(IC,1) 
+          P(NC,2)=P(IC,2) 
+          P(NC,5)=P(IC,5) 
+        ENDIF 
+  140   CONTINUE 
+      ENDIF 
+C...Find initiator cell: the one with highest pT of not yet used ones. 
+      NJ=NC 
+  150 ETMAX=0. 
+      DO 160 IC=N+1,NC 
+      IF(K(IC,5).NE.2) GOTO 160 
+      IF(P(IC,5).LE.ETMAX) GOTO 160 
+      ICMAX=IC 
+      ETA=P(IC,1) 
+      PHI=P(IC,2) 
+      ETMAX=P(IC,5) 
+  160 CONTINUE 
+      IF(ETMAX.LT.PARU(52)) GOTO 220 
+      IF(NJ.GE.MSTU(4)-MSTU(32)-5) THEN 
+        CALL LUERRM(11,'(LUCELL:) no more memory left in LUJETS') 
+        NJET=-2 
+        RETURN 
+      ENDIF 
+      K(ICMAX,5)=1 
+      NJ=NJ+1 
+      K(NJ,4)=0 
+      K(NJ,5)=1 
+      P(NJ,1)=ETA 
+      P(NJ,2)=PHI 
+      P(NJ,3)=0. 
+      P(NJ,4)=0. 
+      P(NJ,5)=0. 
+C...Sum up unused cells within required distance of initiator. 
+      DO 170 IC=N+1,NC 
+      IF(K(IC,5).EQ.0) GOTO 170 
+      IF(ABS(P(IC,1)-ETA).GT.PARU(54)) GOTO 170 
+      DPHIA=ABS(P(IC,2)-PHI) 
+      IF(DPHIA.GT.PARU(54).AND.DPHIA.LT.PARU(2)-PARU(54)) GOTO 170 
+      PHIC=P(IC,2) 
+      IF(DPHIA.GT.PARU(1)) PHIC=PHIC+SIGN(PARU(2),PHI) 
+      IF((P(IC,1)-ETA)**2+(PHIC-PHI)**2.GT.PARU(54)**2) GOTO 170 
+      K(IC,5)=-K(IC,5) 
+      K(NJ,4)=K(NJ,4)+K(IC,4) 
+      P(NJ,3)=P(NJ,3)+P(IC,5)*P(IC,1) 
+      P(NJ,4)=P(NJ,4)+P(IC,5)*PHIC 
+      P(NJ,5)=P(NJ,5)+P(IC,5) 
+  170 CONTINUE 
+C...Reject cluster below minimum ET, else accept. 
+      IF(P(NJ,5).LT.PARU(53)) THEN 
+        NJ=NJ-1 
+        DO 180 IC=N+1,NC 
+        IF(K(IC,5).LT.0) K(IC,5)=-K(IC,5) 
+  180   CONTINUE 
+      ELSEIF(MSTU(54).LE.2) THEN 
+        P(NJ,3)=P(NJ,3)/P(NJ,5) 
+        P(NJ,4)=P(NJ,4)/P(NJ,5) 
+        IF(ABS(P(NJ,4)).GT.PARU(1)) P(NJ,4)=P(NJ,4)-SIGN(PARU(2), 
+     &  P(NJ,4)) 
+        DO 190 IC=N+1,NC 
+        IF(K(IC,5).LT.0) K(IC,5)=0 
+  190   CONTINUE 
+      ELSE 
+        DO 200 J=1,4 
+        P(NJ,J)=0. 
+  200   CONTINUE 
+        DO 210 IC=N+1,NC 
+        IF(K(IC,5).GE.0) GOTO 210 
+        P(NJ,1)=P(NJ,1)+P(IC,5)*COS(P(IC,2)) 
+        P(NJ,2)=P(NJ,2)+P(IC,5)*SIN(P(IC,2)) 
+        P(NJ,3)=P(NJ,3)+P(IC,5)*SINH(P(IC,1)) 
+        P(NJ,4)=P(NJ,4)+P(IC,5)*COSH(P(IC,1)) 
+        K(IC,5)=0 
+  210   CONTINUE 
+      ENDIF 
+      GOTO 150 
+C...Arrange clusters in falling ET sequence. 
+  220 DO 250 I=1,NJ-NC 
+      ETMAX=0. 
+      DO 230 IJ=NC+1,NJ 
+      IF(K(IJ,5).EQ.0) GOTO 230 
+      IF(P(IJ,5).LT.ETMAX) GOTO 230 
+      IJMAX=IJ 
+      ETMAX=P(IJ,5) 
+  230 CONTINUE 
+      K(IJMAX,5)=0 
+      K(N+I,1)=31 
+      K(N+I,2)=98 
+      K(N+I,3)=I 
+      K(N+I,4)=K(IJMAX,4) 
+      K(N+I,5)=0 
+      DO 240 J=1,5 
+      P(N+I,J)=P(IJMAX,J) 
+      V(N+I,J)=0. 
+  240 CONTINUE 
+  250 CONTINUE 
+      NJET=NJ-NC 
+C...Convert to massless or massive four-vectors. 
+      IF(MSTU(54).EQ.2) THEN 
+        DO 260 I=N+1,N+NJET 
+        ETA=P(I,3) 
+        P(I,1)=P(I,5)*COS(P(I,4)) 
+        P(I,2)=P(I,5)*SIN(P(I,4)) 
+        P(I,3)=P(I,5)*SINH(ETA) 
+        P(I,4)=P(I,5)*COSH(ETA) 
+        P(I,5)=0. 
+  260   CONTINUE 
+      ELSEIF(MSTU(54).GE.3) THEN 
+        DO 270 I=N+1,N+NJET 
+        P(I,5)=SQRT(MAX(0.,P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2)) 
+  270   CONTINUE 
+      ENDIF 
+C...Information about storage. 
+      MSTU(61)=N+1 
+      MSTU(62)=NP 
+      MSTU(63)=NC-N 
+      IF(MSTU(43).LE.1) MSTU(3)=NJET 
+      IF(MSTU(43).GE.2) N=N+NJET 
+      RETURN 
+      END