Updates (Leticia)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 23 Apr 2009 08:44:57 +0000 (08:44 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 23 Apr 2009 08:44:57 +0000 (08:44 +0000)
- corrected boosts
- corrected glauber

PYTHIA6/QPYTHIA/q-pyshow.1.0.F

index e2eefe7..e545b7f 100644 (file)
@@ -246,18 +246,14 @@ c     alphas=1, cr=1.
       dimension xkap2(101), xlkap2(101), xome(101), xlome(101)
       dimension xspec(101,101)
       dimension aux1(101), aux2(101)
-      character*1000 filnam
-      character*1000 chroot
       save xkap2, xlkap2, xome, xlome, xspec, npkap, npome
       DATA IFLAG/0/
 c     WE READ THE GRID ONLY THE FIRST TIME.
       IF (IFLAG .EQ. 0) THEN
-         chroot=' '
-         call getenvf('ALICE_ROOT',chroot)
-         lnroot= lnblnk(chroot)
-         filnam=chroot(1:lnroot)//'/PYTHIA6/QPYTHIA/qgrid'
-         write(6,*) "Opening file ", filnam
-         open(11,file=filnam, status='old')
+c         print*, 'reading qgrid'
+         open(11,file=
+     +'/ed22/dfs/work/ALICE/offline/AliRoot/trunk/PYTHIA6/qgrid',
+     +status='old')
          read(11,*) npkap
          read(11,*) npome
          npkap=npkap+1
@@ -571,6 +567,88 @@ c     RI: VALUE OF THE INTEGRAL.
   
       RETURN
       END
+
+
+      SUBROUTINE QPYROBO(XI,YI,ZI,TI,THE,PHI,BEX,BEY,BEZ,XP,YP,ZP,TP)
+C     N. Armesto, 16.04.2009
+C     performs a boost and rotation of (t,x,y,z) to (tp,xp,yp,zp):
+C     cut version of PYROBO, angles and boost parameters identical.
+      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
+C...Local arrays.
+      DIMENSION ROT(3,3),VR(3),DV(4)
+C
+      X=XI
+      Y=YI
+      Z=ZI
+      T=TI
+C...Rotate, typically from z axis to direction (theta,phi).
+      IF(THE**2+PHI**2.GT.1D-20) THEN
+        ROT(1,1)=COS(THE)*COS(PHI)
+        ROT(1,2)=-SIN(PHI)
+        ROT(1,3)=SIN(THE)*COS(PHI)
+        ROT(2,1)=COS(THE)*SIN(PHI)
+        ROT(2,2)=COS(PHI)
+        ROT(2,3)=SIN(THE)*SIN(PHI)
+        ROT(3,1)=-SIN(THE)
+        ROT(3,2)=0D0
+        ROT(3,3)=COS(THE)
+C   Instead of loop 120 in PYROBO.
+        VR(1)=X
+        VR(2)=Y
+        VR(3)=Z
+C   Instead of loop 130 in PYROBO.
+        J=1
+        X=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3)
+        J=2
+        Y=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3)
+        J=3
+        Z=ROT(J,1)*VR(1)+ROT(J,2)*VR(2)+ROT(J,3)*VR(3)
+      ENDIF
+C   If nothing happens...
+      XP=X
+      YP=Y
+      ZP=Z
+      TP=T
+C...Boost, typically from rest to momentum/energy=beta.
+      IF(BEX**2+BEY**2+BEZ**2.GT.1D-20) THEN
+        DBX=BEX
+        DBY=BEY
+        DBZ=BEZ
+        DB=SQRT(DBX**2+DBY**2+DBZ**2)
+        EPS1=1D0-1D-12
+        IF(DB.GT.EPS1) THEN
+C...Rescale boost vector if too close to unity.
+          CALL PYERRM(3,'(PYROBO:) boost vector too large')
+          DBX=DBX*(EPS1/DB)
+          DBY=DBY*(EPS1/DB)
+          DBZ=DBZ*(EPS1/DB)
+          DB=EPS1
+        ENDIF
+        DGA=1D0/SQRT(1D0-DB**2)
+C    Instead of loop 150 in PYROBO.
+        DV(1)=X
+        DV(2)=Y
+        DV(3)=Z
+        DV(4)=T
+        DBV=DBX*DV(1)+DBY*DV(2)+DBZ*DV(3)
+        DGABV=DGA*(DGA*DBV/(1D0+DGA)+DV(4))
+        XP=DV(1)+DGABV*DBX
+        YP=DV(2)+DGABV*DBY
+        ZP=DV(3)+DGABV*DBZ
+        TP=DGA*(DV(4)+DBV)
+      ENDIF
+C 
+      RETURN
+      END
+      
+
+
+
+
+
+
+
+
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C     PYSHOW ROUTINE FOR Q-PYTHIA version 1.0.
@@ -1125,7 +1203,7 @@ c     &-PS(2)/PS(4),-PS(3)/PS(4))
         qpltbx=-PS(1)/PS(4)
         qpltby=-PS(2)/PS(4)
         qpltbz=-PS(3)/PS(4)
-        print*, "caca" 
+        
       ENDIF
 Cacs-
  
@@ -1332,6 +1410,9 @@ Cacs+
          xbx=p(idf,1)/p(idf,4)
          xby=p(idf,2)/p(idf,4)
          xbz=p(idf,3)/p(idf,4)
+
+
+
          ppos(iep(1),1)=ppos(idf,1)+xbx*xlcoh
          ppos(iep(1),2)=ppos(idf,2)+xby*xlcoh
          ppos(iep(1),3)=ppos(idf,3)+xbz*xlcoh
@@ -1340,7 +1421,7 @@ Cacs+
      >               ppos(iep(1),4),xbx,xby,xbz,qhatl,omegac)
       endif
 Cacs-
-  
+     
       DO 390 I=2,NEP
         IEP(I)=IEP(I-1)+1
         IF(IEP(I).GT.N+NEP) IEP(I)=N+1
@@ -2352,11 +2433,21 @@ cforalice+
 c     Here the transverse coordinates of the hard scattering are set by
 c     glauber geometry. 
       call GetRandomXY(xrang,yrang) 
-      x0=xrang ! fm
-      y0=yrang ! fm
+      xin=xrang ! fm
+      yin=yrang ! fm
 cforalice-
-      z0=0.d0 ! fm
-      t0=0.d0 ! fm
+      zin=0.d0 ! fm
+      tin=0.d0 ! fm
+      call qpyrobo(xin,yin,zin,tin,0.d0,0.d0,bbx,bby,bbz,
+     + x1,y1,z1,t1)
+      call qpyrobo(x1,y1,z1,t1,0.d0,aa2,0.d0,0.d0,0.d0,
+     + x2,y2,z2,t2)
+      call qpyrobo(x2,y2,z2,t2,aa1,0.d0,0.d0,0.d0,0.d0,
+     + xout,yout,zout,tout)
+      x0=xout
+      y0=yout
+      z0=zout
+      t0=tout 
       RETURN
       END
 C
@@ -2382,47 +2473,46 @@ C     SCATTERING. THEY ARE THE ENTRIES THREE TO SEVEN IN ROUTINE PYROBO.
       bimp=parj(197)
 c     Here we give five options for the geometry of the medium:
 
-cforalice+
-c     (0) we call routine CalculateI0I1 in AliFastGlauber. Given the position
-c     of the parton in the reaction plane (x,y), the direction in the 
-c     reaction plane phi=atan(py/px) and the impact parameter of the 
-c     collision bimp, it gives back the transverse path length to the 
-c     "end" of the medium and the integrated qhat along that path length. 
-c     See the seter for this option in Configqpythia.C(one can set the 
-c      value of xkscale by doing SetQhat(xkscale), with xkscale in fm). 
-c     The set value is passed here through the pythia free parameter parj(198). 
+c$$$cforalice+
+c$$$     (0) we call routine CalculateI0I1 in AliFastGlauber. Given the position
+c$$$     of the parton in the reaction plane (x,y), the direction in the 
+c$$$     reaction plane phi=atan(py/px) and the impact parameter of the 
+c$$$     collision bimp, it gives back the transverse path length to the 
+c$$$     "end" of the medium and the integrated qhat along that path length. 
+c$$$     See the seter for this option in Configqpythia.C(one can set the 
+c$$$     value of xkscale by doing SetQhat(xkscale), with xkscale in fm). 
+c$$$     The set value is passed here through the pythia free parameter parj(198). 
+
 
       xkscale=parj(198)
       ellcut=20.d0
       xlcero=0.d0
       xlone=0.d0 
+      
+      call qpyrobo(x,y,z,t,-aa1,-aa2,-bbx,-bby,-bbz,xout,
+     + yout,zout,tout)             
+
+      call qpyrobo(bx,by,bz,1.d0,-aa1,-aa2,-bbx,-bby,-bbz,
+     + bx1,by1,bz1,bt1)
+     
 
-      if((bx.eq.0.d0).and.(by.eq.0.d0)) then
-      if(bz.ge.0.d0) then
-         phi=-aa2
-      else
-         phi=dacos(-1.d0)-aa2     
-      endif 
-      else
-      bx1=dcos(aa2)*bx+dsin(-1.d0*aa2)*by
-      by1=dsin(aa2)*bx+dcos(aa2)*by
       phi=datan2(by1,bx1)
-      endif
-    
       phia=phi
-      xa=x
-      ya=y
       bimpa=bimp
       ellcuta=ellcut
+      xa=xout
+      ya=yout
       call CalculateI0I1(xlcero,xlone,bimpa,xa,ya,phia,ellcuta)
       if(xlcero.eq.0.d0) then
            xlp=0.d0
            qhl=0.d0
       else
       xlp=2.d0*xlone/xlcero
-      qhl=0.1973d0*0.1973d0*xlcero*xkscale !GeV**2
-      endif 
-cforalice-
+      qhl=0.1973d0*0.1973d0*xlcero*xkscale 
+      endif
+   
+
+c$$$cforalice-
 c     To use any of these folowing 1,2,3 or 4 options the user should specify
 c     a constant value for the transport coefficient and an initial in medium length. 
 c     This can be done in the user Config file by setting: SetQhat(qhat), with qhat in
@@ -2434,7 +2524,7 @@ c     and comment the other definitions of xlp and qhl above and below.
 c      xlp=xl
 c      if (xlp .lt. 0.d0) xlp=0.d0
 c      qhl=xlp*qhat ! GeV**2
-
+c        print*, xlp,qhl
 c     (2) simplest ansatz: for an initial parton along the z-axis (approximate)
 c      starting in the center of a medium (-xl,+xl) along the z-axis
 c       if (bz .gt. 0.d0) then