lutest.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lutest.f

    
C*********************************************************************  
    
      SUBROUTINE lutest(MTEST)  
    
C...Purpose: to provide a simple program (disguised as subroutine) to   
C...run at installation as a check that the program works as intended.  
      common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
      SAVE /lujets/ 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      dimension psum(5),pini(6),pfin(6) 
    
C...Loop over events to be generated.   
      IF(mtest.GE.1) CALL lutabu(20)    
      nerr=0    
      DO 170 iev=1,600  
    
C...Reset parameter values. Switch on some nonstandard features.    
      mstj(1)=1 
      mstj(3)=0 
      mstj(11)=1    
      mstj(42)=2    
      mstj(43)=4    
      mstj(44)=2    
      parj(17)=0.1  
      parj(22)=1.5  
      parj(43)=1.   
      parj(54)=-0.05    
      mstj(101)=5   
      mstj(104)=5   
      mstj(105)=0   
      mstj(107)=1   
      IF(iev.EQ.301.OR.iev.EQ.351.OR.iev.EQ.401) mstj(116)=3    
    
C...Ten events each for some single jets configurations.    
      IF(iev.LE.50) THEN    
        ity=(iev+9)/10  
        mstj(3)=-1  
        IF(ity.EQ.3.OR.ity.EQ.4) mstj(11)=2 
        IF(ity.EQ.1) CALL lu1ent(1,1,15.,0.,0.) 
        IF(ity.EQ.2) CALL lu1ent(1,3101,15.,0.,0.)  
        IF(ity.EQ.3) CALL lu1ent(1,-2203,15.,0.,0.) 
        IF(ity.EQ.4) CALL lu1ent(1,-4,30.,0.,0.)    
        IF(ity.EQ.5) CALL lu1ent(1,21,15.,0.,0.)    
    
C...Ten events each for some simple jet systems; string fragmentation.  
      ELSEIF(iev.LE.130) THEN   
        ity=(iev-41)/10 
        IF(ity.EQ.1) CALL lu2ent(1,1,-1,40.)    
        IF(ity.EQ.2) CALL lu2ent(1,4,-4,30.)    
        IF(ity.EQ.3) CALL lu2ent(1,2,2103,100.) 
        IF(ity.EQ.4) CALL lu2ent(1,21,21,40.)   
        IF(ity.EQ.5) CALL lu3ent(1,2101,21,-3203,30.,0.6,0.8)   
        IF(ity.EQ.6) CALL lu3ent(1,5,21,-5,40.,0.9,0.8) 
        IF(ity.EQ.7) CALL lu3ent(1,21,21,21,60.,0.7,0.5)    
        IF(ity.EQ.8) CALL lu4ent(1,2,21,21,-2,40.,0.4,0.64,0.6,0.12,0.2)    
    
C...Seventy events with independent fragmentation and momentum cons.    
      ELSEIF(iev.LE.200) THEN   
        ity=1+(iev-131)/16  
        mstj(2)=1+mod(iev-131,4)    
        mstj(3)=1+mod((iev-131)/4,4)    
        IF(ity.EQ.1) CALL lu2ent(1,4,-5,40.)    
        IF(ity.EQ.2) CALL lu3ent(1,3,21,-3,40.,0.9,0.4) 
        IF(ity.EQ.3) CALL lu4ent(1,2,21,21,-2,40.,0.4,0.64,0.6,0.12,0.2)    
        IF(ity.GE.4) CALL lu4ent(1,2,-3,3,-2,40.,0.4,0.64,0.6,0.12,0.2) 
    
C...A hundred events with random jets (check invariant mass).   
      ELSEIF(iev.LE.300) THEN   
  100   DO 110 j=1,5    
  110   psum(j)=0.  
        njet=2.+6.*rlu(0)   
        DO 120 i=1,njet 
        kfl=21  
        IF(i.EQ.1) kfl=int(1.+4.*rlu(0))    
        IF(i.EQ.njet) kfl=-int(1.+4.*rlu(0))    
        ejet=5.+20.*rlu(0)  
        theta=acos(2.*rlu(0)-1.)    
        phi=6.2832*rlu(0)   
        IF(i.LT.njet) CALL lu1ent(-i,kfl,ejet,theta,phi)    
        IF(i.EQ.njet) CALL lu1ent(i,kfl,ejet,theta,phi) 
        IF(i.EQ.1.OR.i.EQ.njet) psum(5)=psum(5)+ulmass(kfl) 
        DO 120 j=1,4    
  120   psum(j)=psum(j)+p(i,j)  
        IF(psum(4)**2-psum(1)**2-psum(2)**2-psum(3)**2.LT.  
     &  (psum(5)+parj(32))**2) goto 100 
    
C...Fifty e+e- continuum events with matrix elements.   
      ELSEIF(iev.LE.350) THEN   
        mstj(101)=2 
        CALL lueevt(0,40.)  
    
C...Fifty e+e- continuum event with varying shower options. 
      ELSEIF(iev.LE.400) THEN   
        mstj(42)=1+mod(iev,2)   
        mstj(43)=1+mod(iev/2,4) 
        mstj(44)=mod(iev/8,3)   
        CALL lueevt(0,90.)  
    
C...Fifty e+e- continuum events with coherent shower, including top.    
      ELSEIF(iev.LE.450) THEN   
        mstj(104)=6 
        CALL lueevt(0,500.) 
    
C...Fifty Upsilon decays to ggg or gammagg with coherent shower.    
      ELSEIF(iev.LE.500) THEN   
        CALL luonia(5,9.46) 
    
C...One decay each for some heavy mesons.   
      ELSEIF(iev.LE.560) THEN   
        ity=iev-501 
        kfls=2*(ity/20)+1   
        kflb=8-mod(ity/5,4) 
        kflc=kflb-mod(ity,5)    
        CALL lu1ent(1,100*kflb+10*kflc+kfls,0.,0.,0.)   
    
C...One decay each for some heavy baryons.  
      ELSEIF(iev.LE.600) THEN   
        ity=iev-561 
        kfls=2*(ity/20)+2   
        kfla=8-mod(ity/5,4) 
        kflb=kfla-mod(ity,5)    
        kflc=max(1,kflb-1)  
        CALL lu1ent(1,1000*kfla+100*kflb+10*kflc+kfls,0.,0.,0.) 
      ENDIF 
    
C...Generate event. Find total momentum, energy and charge. 
      DO 130 j=1,4  
  130 pini(j)=plu(0,j)  
      pini(6)=plu(0,6)  
      CALL luexec   
      DO 140 j=1,4  
  140 pfin(j)=plu(0,j)  
      pfin(6)=plu(0,6)  
    
C...Check conservation of energy, momentum and charge;  
C...usually exact, but only approximate for single jets.    
      merr=0    
      IF(iev.LE.50) THEN    
        IF((pfin(1)-pini(1))**2+(pfin(2)-pini(2))**2.GE.4.) merr=merr+1 
        epzrem=pini(4)+pini(3)-pfin(4)-pfin(3)  
        IF(epzrem.LT.0..OR.epzrem.GT.2.*parj(31)) merr=merr+1   
        IF(abs(pfin(6)-pini(6)).GT.2.1) merr=merr+1 
      ELSE  
        DO 150 j=1,4    
  150   IF(abs(pfin(j)-pini(j)).GT.0001*pini(4)) merr=merr+1    
        IF(abs(pfin(6)-pini(6)).GT.0.1) merr=merr+1 
      ENDIF 
      IF(merr.NE.0) WRITE(mstu(11),1000) (pini(j),j=1,4),pini(6),   
     &(pfin(j),j=1,4),pfin(6)   
    
C...Check that all KF codes are known ones, and that partons/particles  
C...satisfy energy-momentum-mass relation. Store particle statistics.   
      DO 160 i=1,n  
      IF(k(i,1).GT.20) goto 160 
      IF(lucomp(k(i,2)).EQ.0) THEN  
        WRITE(mstu(11),1100) i  
        merr=merr+1 
      ENDIF 
      pd=p(i,4)**2-p(i,1)**2-p(i,2)**2-p(i,3)**2-p(i,5)**2  
      IF(abs(pd).GT.max(0.1,0.001*p(i,4)**2).OR.p(i,4).LT.0.) THEN  
        WRITE(mstu(11),1200) i  
        merr=merr+1 
      ENDIF 
  160 CONTINUE  
      IF(mtest.GE.1) CALL lutabu(21)    
    
C...List all erroneous events and some normal ones. 
      IF(merr.NE.0.OR.mstu(24).NE.0.OR.mstu(28).NE.0) THEN  
        CALL lulist(2)  
      ELSEIF(mtest.GE.1.AND.mod(iev-5,100).EQ.0) THEN   
        CALL lulist(1)  
      ENDIF 
    
C...Stop execution if too many errors. Endresult of run.    
      IF(merr.NE.0) nerr=nerr+1 
      IF(nerr.GE.10) THEN   
        WRITE(mstu(11),1300) iev    
        stop    
      ENDIF 
  170 CONTINUE  
      IF(mtest.GE.1) CALL lutabu(22)    
      WRITE(mstu(11),1400) nerr 
    
C...Reset commonblock variables changed during run. 
      mstj(2)=3 
      parj(17)=0.   
      parj(22)=1.   
      parj(43)=0.5  
      parj(54)=0.   
      mstj(105)=1   
      mstj(107)=0   
    
C...Format statements for output.   
 1000 FORMAT(/' Momentum, energy and/or charge were not conserved ',    
     &'in following event'/' sum of',9x,'px',11x,'py',11x,'pz',11x, 
     &'E',8x,'charge'/' before',2x,4(1x,f12.5),1x,f8.2/' after',3x, 
     &4(1x,f12.5),1x,f8.2)  
 1100 FORMAT(/5x,'Entry no.',i4,' in following event not known code')   
 1200 FORMAT(/5x,'Entry no.',i4,' in following event has faulty ',  
     &'kinematics') 
 1300 FORMAT(/5x,'Ten errors experienced by event ',i3/ 
     &5x,'Something is seriously wrong! Execution stopped now!')    
 1400 FORMAT(/5x,'Number of erroneous or suspect events in run:',i3/    
     &5x,'(0 fine, 1 acceptable if a single jet, ', 
     &'>=2 something is wrong)')    
    
      RETURN    
      END