c HWPLOT c Version 2.0 26/1/2000 c See also http://www.hep.phy.cam.ac.uk/~lester/hwplot.html c Author: c Christopher Lester, St John's College, Cambridge, CB2 1TP, UK. c Thanks for useful discussions to: c James Hetherington, King's College, Cambridge, UK. c lester@hep.phy.cam.ac.uk c jamesh@hep.phy.cam.ac.uk cc--------------------------------------------------------------------- c subroutine exampleAnalysisProgram cc--------------------------------------------------------------------- c include 'HERWIG60.INC' c include 'HWPLOT.INC' c c logical iWantTextLabelsToo c character*(hwpltf) myFilePrefix c integer MaxGen c c ! ... lots of analysis ... c c iWantTextLabelsToo=.false. c myFilePrefix='myfile' c maxGen=5 c c CALL hwplot(maxGen, c & iWantTextLabelsToo, c & myFilePrefix, c & NEVHEP) ! Creates a diagram with up to c ! five generations of perturbative c ! particles. The info is saved in c ! a file called myprefix72.in (say) c c ! ... some more analysis ... c c end c--------------------------------------------------------------------- subroutine hwplot(maxdepth, labels, filename, filedigit) c--------------------------------------------------------------------- include 'HERWIG60.INC' include 'HWPLOT.INC' c-----Input params logical labels integer maxdepth integer filedigit character*(hwpltf) filename c-----Local vars character*(hwpltx) tmpFileName character*(hwpltx) fullFileName character*(hwpltx) tmpDescriptor character*(hwpltx) descriptor c-----Functions integer hwplpu integer hard parameter (hard=6) hwpllabels=labels write (tmpDescriptor, '("(A",I3,",I",I3,",A",I3,")")') & hwpltf,hwpltn,hwplte c print *,tmpDescriptor call StringCompactifier(tmpDescriptor, descriptor) c print *,descriptor write (tmpFileName,descriptor) filename, filedigit, hwpltg c print *,tmpFileName call StringCompactifier(tmpFileName, fullFileName) c print *,fullFileName open(unit=hwplun, file=fullFileName, status='unknown') hwplmd=maxdepth hwplde=1 hwpldy(hwplde)=20 hwplx(hwplde)=1 hwply(hwplde)=10 hwplid(hwplde)=hwplpu(hard) write (hwplun,*) 'mark ',hwplde,' ',hwplx(hwplde), & ' ',hwply(hwplde) call hwplre call hwplai(hwplx(hwplde)+1.0,hwply(hwplde),nevhep) write (hwplun,*) '!' close(hwplun) end c--------------------------------------------------------------------- subroutine hwplotdefault c--------------------------------------------------------------------- include 'HERWIG60.INC' include 'HWPLOT.INC' character*(hwpltf) myfile myfile='temp.' call hwplot(6,.true.,myfile,NEVHEP) end c--------------------------------------------------------------------- subroutine hwplre c--------------------------------------------------------------------- include 'HERWIG60.INC' include 'HWPLOT.INC' c-----Input vars c-----Local vars logical woah integer i logical theyalldid integer typ c-----Functions integer hwplpu logical hwplcl integer hwplgreen integer hwplred integer hwplblack integer hwplblue integer hwplyellow integer hwplcyan integer hwplmagenta integer hwplwhite parameter (hwplred=4) parameter (hwplblack=0) parameter (hwplwhite=7) parameter (hwplcyan=3) parameter (hwplyellow=6) parameter (hwplblue=1) parameter (hwplgreen=2) parameter (hwplmagenta=5) integer hwplzero integer hwplpos integer hwplneg integer hwpldef parameter (hwplzero=hwplblack) parameter (hwplpos=hwplred) parameter (hwplneg=hwplblue) parameter (hwpldef=hwplgreen) integer photon parameter (photon=10) integer gluon parameter (gluon=11) integer dfermion parameter (dfermion=12) integer draw parameter (draw=13) integer fermion parameter (fermion=14) integer wplus parameter (wplus=15) integer wminus parameter (wminus=16) integer eplus parameter (eplus=17) integer eminus parameter (eminus=18) integer nutrino parameter (nutrino=19) integer neutralinoA parameter (neutralinoA=20) integer neutralinoB parameter (neutralinoB=21) integer charginoA parameter (charginoA=22) integer charginoB parameter (charginoB=23) integer gravitino parameter (gravitino=24) integer squark parameter (squark=25) integer gluinoA parameter (gluinoA=26) integer gluinoB parameter (gluinoB=27) integer sneutrino parameter (sneutrino=28) integer negslepton parameter (negslepton=29) integer posslepton parameter (posslepton=30) C formvar(10) = '(''photon #'',I1,'' #'',I1)' 10 format ('photon #',I1,' #',I1) 11 format ('gluon #',I1,' #',I1) 12 format ('dfermion #',I1,' #',I1) 13 format ('draw #',I1,' #',I1) 14 format ('fermion #',I1,' #',I1) 15 format ('cphoton #',I1,' #',I1,' 45.0') 16 format ('cphoton #',I1,' #',I1,' -45.0') 17 format ('cfermion #',I1,' #',I1,' -30.0') 18 format ('cfermion #',I1,' #',I1,' -30.0') 19 format ('fermion #',I1,' #',I1) 20 format ('photon #',I1,' #',I1) 21 format ('draw #',I1,' #',I1) 22 format ('cphoton #',I1,' #',I1,' 45.0') 23 format ('cline #',I1,' #',I1,' 45.0') 24 format ('scalar #',I1,' #',I1) 25 format ('dfermion #',I1,' #',I1) 26 format ('gluon #',I1,' #',I1) 27 format ('draw #',I1,' #',I1) 28 format ('dfermion #',I1,' #',I1) 29 format ('cdfermion #',I1,' #',I1,' -30.0') 30 format ('cdfermion #',I1,' #',I1,' -30.0') if (hwplde.le.hwplmd) then c--------Go as far as one can with the variable hwpld1(hwplde)=JDAHEP(1,hwplid(hwplde)) hwpld2(hwplde)=JDAHEP(2,hwplid(hwplde)) if (hwpld1(hwplde).gt.0) then c-----------At least one particle was made hwplnd(hwplde)=hwpld2(hwplde)-hwpld1(hwplde)+1 hwpllo(hwplde)=hwpld1(hwplde)-1 do while (hwpllo(hwplde).lt.hwpld2(hwplde)) hwpllo(hwplde)=hwpllo(hwplde)+1 hwplid(hwplde+1)=hwplpu(hwpllo(hwplde)) c *(1.1**(hwplde)) hwpldy(hwplde+1)=hwpldy(hwplde)/hwplnd(hwplde) hwply(hwplde+1)=hwply(hwplde)+ & ((( real(hwpllo(hwplde)-hwpld1(hwplde) ) +0.5)/ & real(hwplnd(hwplde)))-0.5) & *hwpldy(hwplde) hwplx(hwplde+1)=hwplx(hwplde)+1+0.4*sqrt(5**2-( & hwply(hwplde+1)-hwply(hwplde) & )**2) write (hwplun,*) 'mark ',hwplde+1, & ' ',hwplx(hwplde+1), & ' ',hwply(hwplde+1) typ=idhep(hwpllo(hwplde)) woah=.false. if (hwpllabels) then call hwplal( & 0.0*hwplx(hwplde)+1.0*hwplx(hwplde+1), & 0.0*hwply(hwplde)+1.0*hwply(hwplde+1), & rname(idhw(hwpllo(hwplde)))) endif if (abs(typ).ge.1000000) then c-----------------Susy particle if (typ.eq.1000021) then c--------------------Gluino write (hwplun,'("colour ",I3)') hwplzero write (hwplun,gluinoA) hwplde,hwplde+1 write (hwplun,gluinoB) hwplde,hwplde+1 else if (abs(typ).eq.1000012.or. & abs(typ).eq.1000014.or. & abs(typ).eq.1000016) then c--------------------Sneutrino write (hwplun,'("colour ",I3)') hwplzero write (hwplun,sneutrino) hwplde,hwplde+1 else if (abs(typ).eq.1000022.or. & abs(typ).eq.1000023.or. & abs(typ).eq.1000025.or. & abs(typ).eq.1000035) then c--------------------Neutralino write (hwplun,'("colour ",I3)') hwplzero if (typ.gt.0) then write (hwplun,neutralinoA) hwplde,hwplde+1 write (hwplun,neutralinoB) hwplde,hwplde+1 else write (hwplun,neutralinoA) hwplde+1,hwplde write (hwplun,neutralinoB) hwplde+1,hwplde endif else if (abs(typ).eq.1000011.or. & abs(typ).eq.1000013.or. & abs(typ).eq.1000015.or. & abs(typ).eq.2000011.or. & abs(typ).eq.2000013.or. & abs(typ).eq.2000015) then c--------------------Slepton if (typ.gt.0) then write (hwplun,'("colour ",I3)') hwplneg write (hwplun,negslepton) hwplde,hwplde+1 else write (hwplun,'("colour ",I3)') hwplpos write (hwplun,posslepton) hwplde+1,hwplde endif else if (abs(typ).eq.1000024.or. & abs(typ).eq.1000037) then c--------------------Chargino if (typ.gt.0) then write (hwplun,'("colour ",I3)') hwplpos write (hwplun,charginoA) hwplde,hwplde+1 write (hwplun,charginoB) hwplde,hwplde+1 else write (hwplun,'("colour ",I3)') hwplneg write (hwplun,charginoA) hwplde+1,hwplde write (hwplun,charginoB) hwplde+1,hwplde endif else if (typ.eq.1000039) then c--------------------Gravitino write (hwplun,'("colour ",I3)') hwplzero write (hwplun,gravitino) hwplde,hwplde+1 else if ((abs(typ).ge.1000001.and. & abs(typ).le.1000006) & .or. & (abs(typ).ge.2000001.and. & abs(typ).le.2000006)) then c--------------------Squark c--------------------These are to be coloured by charge if (typ.gt.0) then if ((typ/2)*2.eq.typ) then write (hwplun,'("colour ",I3)') hwplpos else write (hwplun,'("colour ",I3)') hwplneg endif write (hwplun,squark) hwplde,hwplde+1 else if ((typ/2)*2.eq.typ) then write (hwplun,'("colour ",I3)') hwplneg else write (hwplun,'("colour ",I3)') hwplpos endif write (hwplun,squark) hwplde+1,hwplde endif else print *,'Didnt know ',typ write (hwplun,'("colour ",I3)') hwpldef write (hwplun,dfermion) hwplde,hwplde+1 endif else if (typ.ge.1.and.typ.le.6) then c-----------------Quark if ((typ/2)*2.eq.typ) then write (hwplun,'("colour ",I3)') hwplpos else write (hwplun,'("colour ",I3)') hwplneg endif write (hwplun,fermion) hwplde,hwplde+1 else if (typ.le.-1.and.typ.ge.-6) then c-----------------Anti quark if ((typ/2)*2.eq.typ) then write (hwplun,'("colour ",I3)') hwplneg else write (hwplun,'("colour ",I3)') hwplpos endif write (hwplun,fermion) hwplde+1,hwplde else if (typ.eq.21) then c-----------------Gluon write (hwplun,'("colour ",I3)') hwplzero write (hwplun,gluon) hwplde,hwplde+1 else if (typ.eq.22) then c-----------------Photon write (hwplun,'("colour ",I3)') hwplzero write (hwplun,photon) hwplde,hwplde+1 else if (typ.eq.11.or. & typ.eq.13.or. & typ.eq.15) then c-----------------Neg lepton write (hwplun,'("colour ",I3)') hwplneg write (hwplun,eminus) hwplde,hwplde+1 else if (typ.eq.-11.or. & typ.eq.-13.or. & typ.eq.-15) then c-----------------Pos lepton write (hwplun,'("colour ",I3)') hwplpos write (hwplun,eplus) hwplde+1,hwplde else if (typ.eq.12.or. & typ.eq.14.or. & typ.eq.16) then c-----------------nutrino write (hwplun,'("colour ",I3)') hwplzero write (hwplun,nutrino) hwplde,hwplde+1 else if (typ.eq.-12.or. & typ.eq.-14.or. & typ.eq.-16) then c-----------------anti nutrino write (hwplun,'("colour ",I3)') hwplzero write (hwplun,nutrino) hwplde+1,hwplde else if (abs(typ).eq.24) then c-----------------wplus or minus if (typ.eq.24) then write (hwplun,'("colour ",I3)') hwplpos write (hwplun,wplus) hwplde,hwplde+1 else write (hwplun,'("colour ",I3)') hwplneg write (hwplun,wminus) hwplde,hwplde+1 endif else if (typ.eq.91) then c-----------------Cluster - do nothing. Do not go on. write (hwplun,*) '! Ignored cluster' woah=.true. else c-----------------Unknown print *,'Didnt know ',typ write (hwplun,'("colour ",I3)') hwpldef write (hwplun,draw) hwplde,hwplde+1 endif c--------------Also remember to go no further if this is a c--------------gluon decaying non perturbatively if (typ.eq.21) then c print *,'Testing gluon',hwplid(hwplde+1), c & 'for terbivity.' if (jdahep(1,hwplid(hwplde+1)).ne.0.and. & jdahep(2,hwplid(hwplde+1)).ne.0) then c--------------------The gluon decayed into at least two things c--------------------Did all those things cluster immediately? theyalldid=.true. do i=jdahep(1,hwplid(hwplde+1)), & jdahep(2,hwplid(hwplde+1)) if (.not.hwplcl(i)) then theyalldid=.false. endif enddo if (theyalldid) then c-----------------------All products clustered woah=.true. endif endif endif c--------------Here is the start of the recursive bit if (.not.woah) then hwplde=hwplde+1 call hwplre hwplde=hwplde-1 endif c--------------And there the recursive bit ends. enddo endif endif return end c--------------------------------------------------------------------- subroutine hwplai(digx,digy,lab) c--------------------------------------------------------------------- implicit none include 'HWPLOT.INC' c-----Input variables real digx,digy integer lab real x,y call hwplct(digx,digy,x,y) write (hwplun,10) int(x),int(y) write (hwplun,20) lab 10 format('text ',I5,' ',I5) 20 format('Event ',I10) end c--------------------------------------------------------------------- subroutine hwplal(digx,digy,lab) c--------------------------------------------------------------------- implicit none include 'HWPLOT.INC' c-----Input variables real digx,digy character*8 lab real x,y call hwplct(digx,digy,x,y) write (hwplun,10) int(x),int(y) write (hwplun,20) lab 10 format('text ',I5,' ',I5) 20 format(A8) end c--------------------------------------------------------------------- subroutine hwplct(digx,digy,figx,figy) c--------------------------------------------------------------------- implicit none c-----Input vars real digx,digy c-----Output vars real figx,figy c-----Local stuff real ox,oy,sx,sy parameter (ox=0.0,oy=67.5,sx=450.0,sy=450.0) figx=int(digx*sx+ox) figy=int(digy*sy+oy) end c--------------------------------------------------------------------- function hwplpu(id) c--------------------------------------------------------------------- include 'HERWIG60.INC' integer hwplpu integer id integer d1,d2 logical whoa whoa=.false. hwplpu=id d1=JDAHEP(1,hwplpu) d2=JDAHEP(2,hwplpu) do while (d1.ge.d2.and.d1.gt.0.and.(.not.whoa)) if (idhep(d1).ne.91) then hwplpu=d1 d1=JDAHEP(1,hwplpu) d2=JDAHEP(2,hwplpu) else whoa=.true. endif enddo return end c--------------------------------------------------------------------- function hwplcl(id) c--------------------------------------------------------------------- include 'HERWIG60.INC' logical hwplcl integer hid integer id integer d1,d2 logical whoa integer i logical theyalldid whoa=.false. hid=id d1=JDAHEP(1,hid) d2=JDAHEP(2,hid) hwplcl=.false. do while (d1.ge.d2.and.d1.gt.0.and.(.not.whoa)) if (idhep(d1).ne.91) then hid=d1 d1=JDAHEP(1,hid) d2=JDAHEP(2,hid) else c--------------Possibility of full clustering theyalldid=.true. do i=d1+1,d2 if (idhep(i).ne.91) then theyalldid=.false. endif enddo if (theyalldid) then c-----------------Full clustering! hwplcl=.true. return endif c--------------Non total clustering at best c--------------Do nothing whoa=.true. endif enddo return end c--------------------------------------------------------------------- subroutine StringCompactifier(in, out) c--------------------------------------------------------------------- implicit none include 'HWPLOT.INC' character*(hwpltx) in character*(hwpltx) out character*(hwpltx) buffer integer i buffer=in out =in do i=hwpltx,1,-1 c--------Look for a blank if (buffer(i:i).eq.' ') then c-----------We have a blank, so copy to end out(i:)=buffer(i+1:) buffer=out endif enddo end c---------------------------------------------------------------------