pyhitest.f

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

    
*********************************************************************** 
    
      SUBROUTINE pyhitest(MTEST)  
    
C...Purpose: to provide a simple program (disguised as a 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/ 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /ludat2/ 
      common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)    
      SAVE /ludat3/ 
      common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200) 
      SAVE /pyhisubs/ 
      common/pyhipars/mstp(200),parp(200),msti(200),pari(200) 
      SAVE /pyhipars/ 
    
C...Common initial values. Loop over initiating conditions. 
      mstp(122)=1   
      IF(mtest.LE.0) mstp(122)=0    
      mdcy(lucomp(111),1)=0 
      nerr=0    
      DO 130 iproc=1,7  
    
C...Reset process type, kinematics cuts, and the flags used.    
      msel=0    
      DO 100 isub=1,200 
  100 msub(isub)=0  
      ckin(1)=2.    
      ckin(3)=0.    
      mstp(2)=1 
      mstp(33)=0    
      mstp(81)=1    
      mstp(82)=1    
      mstp(111)=1   
      mstp(131)=0   
      mstp(133)=0   
      parp(131)=0.01    
    
C...Prompt photon production at fixed target.   
      IF(iproc.EQ.1) THEN   
        pzsum=300.  
        pesum=sqrt(pzsum**2+ulmass(211)**2)+ulmass(2212)    
        pqsum=2.    
        msel=10 
        ckin(3)=5.  
        CALL pyhiinit('FIXT','pi+','p',pzsum) 
    
C...QCD processes at ISR energies.  
      ELSEIF(iproc.EQ.2) THEN   
        pesum=63.   
        pzsum=0.    
        pqsum=2.    
        msel=1  
        ckin(3)=5.  
        CALL pyhiinit('CMS','p','p',pesum)    
    
C...W production + multiple interactions at CERN Collider.  
      ELSEIF(iproc.EQ.3) THEN   
        pesum=630.  
        pzsum=0.    
        pqsum=0.    
        msel=12 
        ckin(1)=20. 
        mstp(82)=4  
        mstp(2)=2   
        mstp(33)=3  
        CALL pyhiinit('CMS','p','pbar',pesum) 
    
C...W/Z gauge boson pairs + overlayed events at the Tevatron.   
      ELSEIF(iproc.EQ.4) THEN   
        pesum=1800. 
        pzsum=0.    
        pqsum=0.    
        msub(22)=1  
        msub(23)=1  
        msub(25)=1  
        ckin(1)=200.    
        mstp(111)=0 
        mstp(131)=1 
        mstp(133)=2 
        parp(131)=0.04  
        CALL pyhiinit('CMS','p','pbar',pesum) 
    
C...Higgs production at LHC.    
      ELSEIF(iproc.EQ.5) THEN   
        pesum=17000.    
        pzsum=0.    
        pqsum=0.    
        msel=16 
        pmas(25,1)=300. 
        ckin(1)=200.    
        mstp(81)=0  
        mstp(111)=0 
        CALL pyhiinit('CMS','p','pbar',pesum) 
    
C...Z' production at SSC.   
      ELSEIF(iproc.EQ.6) THEN   
        pesum=40000.    
        pzsum=0.    
        pqsum=0.    
        msel=21 
        pmas(32,1)=600. 
        ckin(1)=400.    
        mstp(81)=0  
        mstp(111)=0 
        CALL pyhiinit('CMS','p','pbar',pesum) 
    
C...W pair production at 1 TeV e+e- collider.   
      ELSEIF(iproc.EQ.7) THEN   
        pesum=1000. 
        pzsum=0.    
        pqsum=0.    
        msub(25)=1  
        CALL pyhiinit('CMS','e+','e-',pesum)  
      ENDIF 
    
C...Generate 20 events of each required type.   
      DO 120 iev=1,20   
      CALL pyhithia   
      pesumm=pesum  
      IF(iproc.EQ.4) pesumm=msti(41)*pesum  
    
C...Check conservation of energy/momentum/flavour.  
      merr=0    
      deve=abs(plu(0,4)-pesumm)+abs(plu(0,3)-pzsum) 
      devt=abs(plu(0,1))+abs(plu(0,2))  
      devq=abs(plu(0,6)-pqsum)  
      IF(deve.GT.1e-3*pesum.OR.devt.GT.max(0.01,1e-5*pesum).OR. 
     &devq.GT.0.1) merr=1   
      IF(merr.NE.0) WRITE(mstu(11),1000) iproc,iev  
    
C...Check that all KF codes are known ones, and that partons/particles  
C...satisfy energy-momentum-mass relation.  
      DO 110 i=1,n  
      IF(k(i,1).GT.20) goto 110 
      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* 
     &sign(1.,p(i,5))   
      IF(abs(pd).GT.max(0.1,0.002*p(i,4)**2,0.002*p(i,5)**2).OR.    
     &(p(i,5).GE.0..AND.p(i,4).LT.0.)) THEN 
        WRITE(mstu(11),1200) i  
        merr=merr+1 
      ENDIF 
  110 CONTINUE  
    
C...Listing of erronoeus events, and first event of each type.  
      IF(merr.GE.1) nerr=nerr+1 
      IF(nerr.GE.10) THEN   
        WRITE(mstu(11),1300)    
        CALL lulist(1)  
        stop    
      ENDIF 
      IF(mtest.GE.1.AND.(merr.GE.1.OR.iev.EQ.1)) THEN   
        IF(merr.GE.1) WRITE(mstu(11),1400)  
        CALL lulist(1)  
      ENDIF 
  120 CONTINUE  
    
C...List statistics for each process type.  
      IF(mtest.GE.1) CALL pyhistat(1) 
  130 CONTINUE  
    
C...Summarize result of run.    
      IF(nerr.EQ.0) WRITE(mstu(11),1500)    
      IF(nerr.GT.0) WRITE(mstu(11),1600) nerr   
      RETURN    
    
C...Formats for information.    
 1000 FORMAT(/5x,'Energy/momentum/flavour nonconservation for process', 
     &i2,', event',i4)  
 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,'This is the tenth error experienced! Something is ',  
     &'wrong.'/5x,'Execution will be stopped after listing of event.')  
 1400 FORMAT(5x,'Faulty event follows:')    
 1500 FORMAT(//5x,'End result of run: no errors detected.') 
 1600 FORMAT(//5x,'End result of run:',i2,' errors detected.'/  
     &5x,'This should not have happened!')  
      END