ReggeGribovPartonMCInterface/src/epos-ids.f

0001 c-----------------------------------------------------------------------
0002       subroutine iclass(id,icl)
0003 c-----------------------------------------------------------------------
0004 c      determines hadron class
0005 c-----------------------------------------------------------------------
0006       ida=iabs(id)
0007       if(ida.eq.0.or.(ida.ge.17.and.ida.le.19))then
0008        icl=2
0009       elseif(ida.eq.130.or.ida.eq.230.or.ida.eq.20)then
0010        icl=3
0011       elseif(ida.eq.140.or.ida.eq.240.or.ida.eq.340.or.ida/10.eq.44)then
0012        icl=4
0013       elseif(ida.ge.100.and.ida.le.999)then
0014        icl=1
0015       elseif(ida.ge.1000.and.ida.le.9999)then
0016        icl=2
0017       else
0018        stop'iclass: id not known'
0019       endif
0020       end
0021
0022 c-----------------------------------------------------------------------
0023       subroutine idchrg(id,chrg)
0024 c     computes charge of particle with ident code id
0025 c     ichrg must be dimensioned nqlep+12
0026 c-----------------------------------------------------------------------
0027       dimension ichrg(53),ifl(3)
0028       data ichrg/0,2,-1,-1,2,-1,2,-1,2,0,0,0,-3,0,-3,0,-3,1,1,2,2*0
0029      *,2,-1,-1,2,-1,2,-1,2,0,0,0,-3,0,-3,0,-3,0,-3,3,0
0030      *,3,0,0,0,3,3,3,6,6,6,0/
0031       idabs=iabs(id)
0032       call idflav(id,ifl(1),ifl(2),ifl(3),jspin,ind)
0033       if(idabs.lt.100) goto 200
0034       isum=0
0035       do 100 i=1,3
0036       if(abs(ifl(i)).gt.52)goto 100
0037       isum=isum+ichrg(iabs(ifl(i))+1)*isign(1,ifl(i))
0038   100 continue
0039       chrg=isum/3.
0040       return
0041 200   chrg=ichrg(ind+1)*isign(1,id)
0042       chrg=chrg/3.
0043       return
0044       end
0045
0046 c-----------------------------------------------------------------------
0047       subroutine idspin(id,iso,jspin,istra)
0048 c     computes iso (isospin), jspin and istra (strangeness) of particle
0049 c     with ident code id
0050 c-----------------------------------------------------------------------
0051       include 'epos.inc'
0052       dimension ifl(3)
0053       iso=0
0054       jspin=0
0055       istra=0
0056       idabs=abs(id)
0057       if(idabs.le.nflav)then
0058         iso=isospin(idabs)*sign(1,id)
0059         if(idabs.ge.3)istra=sign(1,id)
0060         return
0061       endif
0062       call idflav(id,ifl(1),ifl(2),ifl(3),jspin,ind)
0063       iq1=abs(ifl(1))
0064       iq2=abs(ifl(2))
0065       iq3=abs(ifl(3))
0066       if(iq1.ge.3)istra=istra+sign(1,ifl(1))
0067       if(iq2.ge.3)istra=istra+sign(1,ifl(2))
0068       if(iq3.ge.3)istra=istra+sign(1,ifl(3))
0069       if(iq1.ne.0)then         !baryon
0070         iso=(isospin(iq1)+isospin(iq2)+isospin(iq3))*sign(1,iq1)
0071       else
0072         iso=isospin(iq2)*sign(1,ifl(2))
0073         iso=iso+isospin(iq3)*sign(1,ifl(3))
0074       endif
0075       return
0076       end
0077
0078 c-----------------------------------------------------------------------
0079       subroutine idcomk(ic)
0080 c     compactifies ic
0081 c-----------------------------------------------------------------------
0082       parameter (nflav=6)
0083       integer ic(2),icx(2),jc(nflav,2)
0084       call idcomp(ic,icx,jc,1)
0085       ic(1)=icx(1)
0086       ic(2)=icx(2)
0087       return
0088       end
0089
0090 cc-----------------------------------------------------------------------
0091 c      subroutine idcomi(ic,icx)
0092 cc     compactifies ic
0093 cc-----------------------------------------------------------------------
0094 c      parameter (nflav=6)
0095 c      integer ic(2),icx(2),jc(nflav,2)
0096 c      call idcomp(ic,icx,jc,1)
0097 c      return
0098 c      end
0099 c
0100 c-----------------------------------------------------------------------
0101       subroutine idcomj(jc)
0102 c     compactifies jc
0103 c-----------------------------------------------------------------------
0104       parameter (nflav=6)
0105       integer ic(2),icx(2),jc(nflav,2)
0106       call idcomp(ic,icx,jc,2)
0107       return
0108       end
0109
0110 c-----------------------------------------------------------------------
0111       subroutine idcomp(ic,icx,jc,im)
0112 c-----------------------------------------------------------------------
0113 c     compactifies ic,jc
0114 c     input: im (1 or 2)
0115 c            ic (if im=1)
0116 c            jc (if im=2)
0117 c     output: icx (if im=1)
0118 c             jc
0119 c-----------------------------------------------------------------------
0120       parameter (nflav=6)
0121       integer ic(2),icx(2),jc(nflav,2)
0122       if(im.eq.1)call iddeco(ic,jc)
0123       icx(1)=0
0124       icx(2)=0
0125            do n=1,nflav
0126            do j=1,2
0127       if(jc(n,j).ne.0)goto1
0128            enddo
0129            enddo
0130       return
0131 1     continue
0132       nq=0
0133       na=0
0134            do n=1,nflav
0135       nq=nq+jc(n,1)
0136       na=na+jc(n,2)
0137            enddo
0138       l=0
0139            do n=1,nflav
0140       k=min0(jc(n,1),jc(n,2))
0141       if(nq.eq.1.and.na.eq.1)k=0
0142       jc(n,1)=jc(n,1)-k
0143       jc(n,2)=jc(n,2)-k
0144       if(jc(n,1).lt.0.or.jc(n,2).lt.0)
0145      *call utstop('idcomp: jc negative&')
0146       l=l+jc(n,1)+jc(n,2)
0147            enddo
0148            if(l.eq.0)then
0149       jc(1,1)=1
0150       jc(1,2)=1
0151            endif
0152            if(im.eq.1)then
0153       call idenco(jc,icx,ireten)
0154       if(ireten.eq.1)call utstop('idcomp: idenco ret code = 1&')
0155            endif
0156       return
0157       end
0158
0159 c-----------------------------------------------------------------------
0160       subroutine iddeco(ic,jc)
0161 c     decode particle id
0162 c-----------------------------------------------------------------------
0163       parameter (nflav=6)
0164       integer jc(nflav,2),ic(2)
0165       ici=ic(1)
0166       jc(6,1)=mod(ici,10)
0167       jc(5,1)=mod(ici/10,10)
0168       jc(4,1)=mod(ici/100,10)
0169       jc(3,1)=mod(ici/1000,10)
0170       jc(2,1)=mod(ici/10000,10)
0171       jc(1,1)=mod(ici/100000,10)
0172       ici=ic(2)
0173       jc(6,2)=mod(ici,10)
0174       jc(5,2)=mod(ici/10,10)
0175       jc(4,2)=mod(ici/100,10)
0176       jc(3,2)=mod(ici/1000,10)
0177       jc(2,2)=mod(ici/10000,10)
0178       jc(1,2)=mod(ici/100000,10)
0179       return
0180       end
0181
0182 c-----------------------------------------------------------------------
0183       subroutine idenco(jc,ic,ireten)
0184 c     encode particle id
0185 c-----------------------------------------------------------------------
0186       parameter (nflav=6)
0187       integer jc(nflav,2),ic(2)
0188       ireten=0
0189       ic(1)=0
0190       do 20 i=1,nflav
0191       if(jc(i,1).ge.10)goto22
0192 20    ic(1)=ic(1)+jc(i,1)*10**(nflav-i)
0193       ic(2)=0
0194       do 21 i=1,nflav
0195       if(jc(i,2).ge.10)goto22
0196 21    ic(2)=ic(2)+jc(i,2)*10**(nflav-i)
0197       return
0198 22    ireten=1
0199       ic(1)=0
0200       ic(2)=0
0201       return
0202       end
0203
0204 c-----------------------------------------------------------------------
0205       subroutine idenct(jc,id,ib1,ib2,ib3,ib4)
0206 c     encode particle id
0207 c-----------------------------------------------------------------------
0208       parameter (nflav=6)
0209       integer jc(nflav,2),ic(2)
0210
0211       do 40 nf=1,nflav
0212       do 40 ij=1,2
0213       if(jc(nf,ij).ge.10)id=7*10**8
0214 40    continue
0215            if(id/10**8.ne.7)then
0216       call idenco(jc,ic,ireten)
0217       if(ireten.eq.1)call utstop('idenct: idenco ret code = 1&')
0218       if(mod(ic(1),100).ne.0.or.mod(ic(2),100).ne.0)then
0219       id=9*10**8
0220       else
0221       id=8*10**8+ic(1)*100+ic(2)/100
0222       endif
0223            else
0224       call idtrbi(jc,ib1,ib2,ib3,ib4)
0225       id=id
0226      *+mod(jc(1,1)+jc(2,1)+jc(3,1)+jc(4,1),10**4)*10**4
0227      *+mod(jc(1,2)+jc(2,2)+jc(3,2)+jc(4,2),10**4)
0228            endif
0229       return
0230       end
0231
0232 c-----------------------------------------------------------------------
0233       subroutine idflav(id,ifl1,ifl2,ifl3,jspin,index)
0234 c     unpacks the ident code id=+/-ijkl
0235 c
0236 c          mesons--
0237 c          i=0, j<=k, +/- is sign for j
0238 c          id=110 for pi0, id=220 for eta, etc.
0239 c
0240 c          baryons--
0241 c          i<=j<=k in general
0242 c          j<i<k for second state antisymmetric in (i,j), eg. l = 2130
0243 c
0244 c          other--
0245 c          id=1,...,6 for quarks
0246 c          id=9 for gluon
0247 c          id=10 for photon
0248 c          id=11,...,16 for leptons
0249 c          i=17 for deuteron
0250 c          i=18 for triton
0251 c          i=19 for alpha
0252 c          id=20 for ks, id=-20 for kl
0253 c
0254 c          i=21...26 for scalar quarks
0255 c          i=29 for gluino
0256 c
0257 c          i=30 for h-dibaryon
0258 c
0259 c          i=31...36 for scalar leptons
0260 c          i=39 for wino
0261 c          i=40 for zino
0262 c
0263 c          id=80 for w+
0264 c          id=81,...,83 for higgs mesons (h0, H0, A0, H+)
0265 c          id=84,...,87 for excited bosons (Z'0, Z''0, W'+)
0266 c          id=90 for z0
0267 c
0268 c          diquarks--
0269 c          id=+/-ij00, i<j for diquark composed of i,j.
0270 c
0271 c
0272 c          index is a sequence number used internally
0273 c          (index=-1 if id doesn't exist)
0274 c
0275 c-----------------------------------------------------------------------
0276       parameter ( nqlep=41,nmes=2)
0277       ifl1=0
0278       ifl2=0
0279       ifl3=0
0280       jspin=0
0281       idabs=iabs(id)
0282       i=idabs/1000
0283       j=mod(idabs/100,10)
0284       k=mod(idabs/10,10)
0285       jspin=mod(idabs,10)
0286       if(id.ge.10000) goto 400
0287       if(id.ne.0.and.mod(id,100).eq.0) goto 300
0288       if(j.eq.0) goto 200
0289       if(i.eq.0) goto 100
0290 c          baryons
0291 c          only x,y baryons are qqx, qqy, q=u,d,s.
0292       ifl1=isign(i,id)
0293       ifl2=isign(j,id)
0294       ifl3=isign(k,id)
0295       if(k.le.6) then
0296         index=max0(i-1,j-1)**2+i+max0(i-j,0)+(k-1)*k*(2*k-1)/6
0297      1  +109*jspin+36*nmes+nqlep+11
0298       else
0299         index=max0(i-1,j-1)**2+i+max0(i-j,0)+9*(k-7)+91
0300      1  +109*jspin+36*nmes+nqlep+11
0301       endif
0302       return
0303 c          mesons
0304 100   continue
0305       ifl1=0
0306       ifl2=isign(j,id)
0307       ifl3=isign(k,-id)
0308       index=j+k*(k-1)/2+36*jspin+nqlep
0309       index=index+11
0310       return
0311 c          quarks, leptons, etc
0312 200   continue
0313       ifl1=0
0314       ifl2=0
0315       ifl3=0
0316       jspin=0
0317       index=idabs
0318       if(idabs.lt.20) return
0319 c          define index=20 for ks, index=21 for kl
0320       index=idabs+1
0321       if(id.eq.20) index=20
0322 c          index=nqlep+1,...,nqlep+11 for w+, higgs, z0
0323       if(idabs.lt.80) return
0324       index=nqlep+idabs-79
0325       return
0326 300   ifl1=isign(i,id)
0327       ifl2=isign(j,id)
0328       ifl3=0
0329       jspin=0
0330       index=0
0331       return
0332  400  index=-1
0333       return
0334       end
0335
0336 c-----------------------------------------------------------------------
0337       subroutine idqufl(n,id,nqu,nqd,nqs)
0338 c     unpacks the ident code of particle (n) and give the number of
0339 c     quarks of each flavour(only u,d,s)
0340 c-----------------------------------------------------------------------
0341       include 'epos.inc'
0342       include 'epos.incems'
0343       integer jc(nflav,2),ic(2)
0344
0345       nqu=0
0346       nqd=0
0347       nqs=0
0348       if(iabs(id).ge.7.and.iabs(id).lt.100.and.iabs(id).ne.20)return
0349       if(iabs(id)/10.eq.11.or.iabs(id)/10.eq.22)return
0350       if(iabs(id).eq.20)then
0351         if(iorptl(n).gt.0)then
0352           if(idptl(iorptl(n)).gt.0)then
0353             nqd=1
0354             nqs=-1
0355           else
0356             nqd=-1
0357             nqs=1
0358           endif
0359         else
0360           if(ish.ge.4)write(ifch,*)'Cannot count the number of quark'
0361         endif
0362         return
0363       endif
0364       if(id.ne.0.and.mod(id,100).eq.0.and.id.le.10**8) goto 300
0365       if(id/10**8.ne.7)then
0366         call idtr4(id,ic)
0367         call iddeco(ic,jc)
0368       else
0369         call idtrb(ibptl(1,n),ibptl(2,n),ibptl(3,n),ibptl(4,n),jc)
0370       endif
0371       nqu=jc(1,1)-jc(1,2)
0372       nqd=jc(2,1)-jc(2,2)
0373       nqs=jc(3,1)-jc(3,2)
0374       return
0375  300  i=iabs(id)/1000
0376       j=mod(iabs(id)/100,10)
0377       ifl1=isign(i,id)
0378       ifl2=isign(j,id)
0379       if(iabs(ifl1).eq.1)nqu=isign(1,ifl1)
0380       if(iabs(ifl1).eq.2)nqd=isign(1,ifl1)
0381       if(iabs(ifl1).eq.3)nqs=isign(1,ifl1)
0382       if(iabs(ifl2).eq.1)nqu=nqu+isign(1,ifl2)
0383       if(iabs(ifl2).eq.2)nqd=nqd+isign(1,ifl2)
0384       if(iabs(ifl2).eq.3)nqs=nqs+isign(1,ifl2)
0385 c      write(ifch,*)'id',id,ifl1,ifl2,nqu,nqd,nqs
0386       return
0387       end
0388
0389 c-----------------------------------------------------------------------
0390       function idlabl(id)
0391 c     returns the character*8 label for the particle id
0392 c-----------------------------------------------------------------------
0393       parameter ( nqlep=41,nmes=2)
0394 c
0395       character*8 idlabl
0396       character*8 llep,lmes0,lmes1,lbar0,labar0,lbar1,labar1
0397       character*8 lqq,laqq
0398       dimension llep(104)
0399       dimension lmes0(64),lmes1(64)
0400       dimension lbar0(109),labar0(109),lbar1(109),labar1(109)
0401       dimension lqq(21),laqq(21)
0402 c          diquark labels
0403       data lqq/
0404      1'uu0. ','ud0. ','dd0. ','us0. ','ds0. ','ss0. ','uc0. ','dc0. ',
0405      2'sc0. ','cc0. ','ub0. ','db0. ','sb0. ','cb0. ','bb0. ','ut0. ',
0406      3'dt0. ','st0. ','ct0. ','bt0. ','tt0. '/
0407       data laqq/
0408      1'auu0.','aud0.','add0.','aus0.','ads0.','ass0.','auc0.','adc0.',
0409      2'asc0.','acc0.','aub0.','adb0.','asb0.','acb0.','abb0.','aut0.',
0410      3'adt0.','ast0.','act0.','abt0.','att0.'/
0411 c          quark and lepton labels
0412       data llep/
0413      *'     ','up   ','ub   ','dn   ','db   ','st   ','sb   ','ch   ',
0414      *'cb   ','bt   ','bb   ','tp   ','tb   ','y    ','yb   ','x    ',
0415      *'xb   ','gl   ','err  ','gm   ','err  ','nue  ','anue ','e-   ',
0416      *'e+   ','num  ','anum ','mu-  ','mu+  ','nut  ','anut ','tau- ',
0417      *'tau+ ','deut ','adeut','trit ','atrit','alph ','aalph','ks   ',
0418      *'err  ','err  ','kl   ',
0419      *'upss ','ubss ','dnss ','dbss ','stss ','sbss ','chss ','cbss ',
0420      *'btss ','bbss ','tpss ','tbss ','err  ','err  ','err  ','err  ',
0421      *'glss ','err  ','hdiba','err  ','ness ','aness','e-ss ','e+ss ',
0422      *'nmss ','anmss','mu-ss','mu+ss','ntss ','antss','t-ss ','t+ss ',
0423      *'err  ','err  ','err  ','err  ','w+ss ','w-ss ','z0ss ','err  ',
0424      *'w+   ','w-   ','h0   ','ah0  ','H0   ','aH0  ','A0   ','aA0  ',
0425      *'H+   ','H-   ','Zp0  ','aZp0 ','Zpp0 ','aZpp0','Wp+  ','Wp-  ',
0426      *'err  ','err  ','err  ','err  ','z0   '/
0427 c          0- meson labels
0428       data lmes0/
0429      1'pi0  ','pi+  ','eta  ','pi-  ','k+   ','k0   ','etap ','ak0  ',
0430      2'k-   ','ad0  ','d-   ','f-   ','etac ','f+   ','d+   ','d0   ',
0431      2'ub.  ','db.  ','sb.  ','cb.  ','bb.  ','bc.  ','bs.  ','bd.  ',
0432      3'bu.  ','ut.  ','dt.  ','st.  ','ct.  ','bt.  ','tt.  ','tb.  ',
0433      4'tc.  ','ts.  ','td.  ','tu.  ','uy.  ','dy.  ','sy.  ','cy.  ',
0434      5'by.  ','ty.  ','yy.  ','yt.  ','yb.  ','yc.  ','ys.  ','yd.  ',
0435      6'yu.  ','ux.  ','dx.  ','sx.  ','cx.  ','bx.  ','tx.  ','yx.  ',
0436      7'xx.  ','xy.  ','xt.  ','xb.  ','xc.  ','xs.  ','xd.  ','xu.  '/
0437 c          1- meson labels
0438       data lmes1/
0439      1'rho0 ','rho+ ','omeg ','rho- ','k*+  ','k*0  ','phi  ','ak*0 ',
0440      2'k*-  ','ad*0 ','d*-  ','f*-  ','jpsi ','f*+  ','d*+  ','d*0  ',
0441      3'ub*  ','db*  ','sb*  ','cb*  ','upsl ','bc*  ','bs*  ','bd*  ',
0442      4'bu*  ','ut*  ','dt*  ','st*  ','ct*  ','bt*  ','tt*  ','tb*  ',
0443      5'tc*  ','ts*  ','td*  ','tu*  ','uy*  ','dy*  ','sy*  ','cy*  ',
0444      6'by*  ','ty*  ','yy*  ','yt*  ','yb*  ','yc*  ','ys*  ','yd*  ',
0445      7'yu*  ','ux*  ','dx*  ','sx*  ','cx*  ','bx*  ','tx*  ','yx*  ',
0446      8'xx*  ','xy*  ','xt*  ','xb*  ','xc*  ','xs*  ','xd*  ','xu*  '/
0447 c          1/2+ baryon labels
0448       data lbar0/
0449      1'err  ','p    ','n    ','err  ','err  ','s+   ','s0   ','s-   ',
0450      2'l    ','xi0  ','xi-  ','err  ','err  ','err  ','sc++ ','sc+  ',
0451      3'sc0  ','lc+  ','usc. ','dsc. ','ssc. ','sdc. ','suc. ','ucc. ',
0452      4'dcc. ','scc. ','err  ','err  ','err  ','err  ','uub. ','udb. ',
0453      5'ddb. ','dub. ','usb. ','dsb. ','ssb. ','sdb. ','sub. ','ucb. ',
0454      6'dcb. ','scb. ','ccb. ','csb. ','cdb. ','cub. ','ubb. ','dbb. ',
0455      7'sbb. ','cbb. ','err  ','err  ','err  ','err  ','err  ','utt. ',
0456      8'udt. ','ddt. ','dut. ','ust. ','dst. ','sst. ','sdt. ','sut. ',
0457      9'uct. ','dct. ','sct. ','cct. ','cst. ','cdt. ','cut. ','ubt. ',
0458      1'dbt. ','sbt. ','cbt. ','bbt. ','bct. ','bst. ','bdt. ','but. ',
0459      2'utt. ','dtt. ','stt. ','ctt. ','btt. ','err  ','err  ','err  ',
0460      3'err  ','err  ','err  ','uuy. ','udy. ','ddy. ','duy. ','usy. ',
0461      4'dsy. ','ssy. ','sdy. ','suy. ','uux. ','udx. ','ddx. ','dux. ',
0462      5'usx. ','dsx. ','ssx. ','sdx. ','sux. '/
0463       data labar0/
0464      1'err  ','ap   ','an   ','err  ','err  ','as-  ','as0  ','as+  ',
0465      2'al   ','axi0 ','axi+ ','err  ','err  ','err  ','asc--','asc- ',
0466      3'asc0 ','alc- ','ausc.','adsc.','assc.','asdc.','asuc.','aucc.',
0467      4'adcc.','ascc.','err  ','err  ','err  ','err  ','auub.','audb.',
0468      5'addb.','adub.','ausb.','adsb.','assb.','asdb.','asub.','aucb.',
0469      6'adcb.','ascb.','accb.','acsb.','acdb.','acub.','aubb.','adbb.',
0470      7'asbb.','acbb.','err  ','err  ','err  ','err  ','err  ','autt.',
0471      8'audt.','addt.','adut.','aust.','adst.','asst.','asdt.','asut.',
0472      9'auct.','adct.','asct.','acct.','acst.','acdt.','acut.','aubt.',
0473      1'adbt.','asbt.','acbt.','abbt.','abct.','abst.','abdt.','abut.',
0474      2'autt.','adtt.','astt.','actt.','abtt.','err  ','err  ','err  ',
0475      3'err  ','err  ','err  ','auuy.','audy.','addy.','aduy.','ausy.',
0476      4'adsy.','assy.','asdy.','asuy.','auux.','audx.','addx.','adux.',
0477      5'ausx.','adsx.','assx.','asdx.','asux.'/
0478 c          3/2+ baryon labels
0479       data lbar1/
0480      1'dl++ ','dl+  ','dl0  ','dl-  ','err  ','s*+  ','s*0  ','s*-  ',
0481      2'err  ','xi*0 ','xi*- ','om-  ','err  ','err  ','uuc* ','udc* ',
0482      3'ddc* ','err  ','usc* ','dsc* ','ssc* ','err  ','err  ','ucc* ',
0483      4'dcc* ','scc* ','ccc* ','err  ','err  ','err  ','uub* ','udb* ',
0484      5'ddb* ','err  ','usb* ','dsb* ','ssb* ','err  ','err  ','ucb* ',
0485      6'dcb* ','scb* ','ccb* ','err  ','err  ','err  ','ubb* ','dbb* ',
0486      7'sbb* ','cbb* ','bbb* ','err  ','err  ','err  ','err  ','utt* ',
0487      8'udt* ','ddt* ','err  ','ust* ','dst* ','sst* ','err  ','err  ',
0488      9'uct* ','dct* ','sct* ','cct* ','err  ','err  ','err  ','ubt* ',
0489      1'dbt* ','sbt* ','cbt* ','bbt* ','err  ','err  ','err  ','err  ',
0490      2'utt* ','dtt* ','stt* ','ctt* ','btt* ','ttt* ','err  ','err  ',
0491      3'err  ','err  ','err  ','uuy* ','udy* ','ddy* ','err  ','usy* ',
0492      4'dsy* ','ssy* ','err  ','err  ','uux* ','udx* ','ddx* ','err  ',
0493      5'usx* ','dsx* ','ssx* ','err  ','err  '/
0494       data labar1/
0495      1'adl--','adl- ','adl0 ','adl+ ','err  ','as*- ','as*0 ','as*+ ',
0496      2'err  ','axi*0','axi*+','aom+ ','err  ','err  ','auuc*','audc*',
0497      3'addc*','err  ','ausc*','adsc*','assc*','err  ','err  ','aucc*',
0498      4'adcc*','ascc*','accc*','err  ','err  ','err  ','auub*','audb*',
0499      5'addb*','err  ','ausb*','adsb*','assb*','err  ','err  ','aucb*',
0500      6'adcb*','ascb*','accb*','err  ','err  ','err  ','aubb*','adbb*',
0501      7'asbb*','acbb*','abbb*','err  ','err  ','err  ','err  ','autt*',
0502      8'audt*','addt*','err  ','aust*','adst*','asst*','err  ','err  ',
0503      9'auct*','adct*','asct*','acct*','err  ','err  ','err  ','aubt*',
0504      1'adbt*','asbt*','acbt*','abbt*','err  ','err  ','err  ','err  ',
0505      2'autt*','adtt*','astt*','actt*','abtt*','attt*','err  ','err  ',
0506      3'err  ','err  ','err  ','auuy*','audy*','addy*','err  ','ausy*',
0507      4'adsy*','assy*','err  ','err  ','auux*','audx*','addx*','err  ',
0508      5'ausx*','adsx*','assx*','err  ','err  '/
0509 c          entry
0510       call idflav(id,ifl1,ifl2,ifl3,jspin,ind)
0511       if(iabs(id).lt.100) goto200
0512       if(iabs(id).lt.1000) goto100
0513       if(id.ne.0.and.mod(id,100).eq.0) goto300
0514 c          baryons
0515       ind=ind-109*jspin-36*nmes-nqlep
0516       ind=ind-11
0517       if(jspin.eq.0.and.id.gt.0) idlabl=lbar0(ind)
0518       if(jspin.eq.0.and.id.lt.0) idlabl=labar0(ind)
0519       if(jspin.eq.1.and.id.gt.0) idlabl=lbar1(ind)
0520       if(jspin.eq.1.and.id.lt.0) idlabl=labar1(ind)
0521       return
0522 c          mesons
0523 100   continue
0524       i=max0(ifl2,ifl3)
0525       j=-min0(ifl2,ifl3)
0526       ind=max0(i-1,j-1)**2+i+max0(i-j,0)
0527       if(jspin.eq.0) idlabl=lmes0(ind)
0528       if(jspin.eq.1) idlabl=lmes1(ind)
0529       return
0530 c          quarks, leptons, etc.
0531 200   continue
0532       ind=2*ind
0533       if(id.le.0) ind=ind+1
0534       idlabl=llep(ind)
0535       return
0536 300   i=iabs(ifl1)
0537       j=iabs(ifl2)
0538       ind=i+j*(j-1)/2
0539       if(id.gt.0) idlabl=lqq(ind)
0540       if(id.lt.0) idlabl=laqq(ind)
0541       return
0542       end
0543
0544 c-----------------------------------------------------------------------
0545       subroutine idmass(idi,amass)
0546 c     returns the mass of the particle with ident code id.
0547 c     (deuteron, triton and alpha mass come from Gheisha ???)
0548 c-----------------------------------------------------------------------
0549       dimension ammes0(15),ammes1(15),ambar0(30),ambar1(30)
0550       dimension amlep(52)
0551       parameter ( nqlep=41,nmes=2)
0552 c-c   data amlep/.3,.3,.5,1.6,4.9,30.,-1.,-1.,0.,0.,
0553       data amlep/.005,.009,.180,1.6,4.9,170.,-1.,-1.,0.,0.,0.
0554      *     ,.5109989e-3,0.,.105658,0.,1.777,1.87656,2.8167,3.755,.49767
0555      *     ,.49767,100.3,100.3,100.5,101.6,104.9,130.,2*-1.,100.,0.,
0556      *     100.,100.005,100.,100.1,100.,101.8,2*-1.,100.,100.,
0557      *     11*0./
0558 c          0- meson mass table
0559       data ammes0/.1349766,.13957018,.547853       !pi0,pi+-,eta
0560      *           ,.493677,.497614,.95778           !K+-, K0,eta'
0561      *    ,1.86483,1.86960,1.96847,2.9803          !D0,D+-,Ds,etac
0562      1    ,5.27917,5.27950,5.3663,6.277,9.390/     !B+-,B0,Bs,Bc,etab
0563 c     1- meson mass table
0564       data ammes1/.77549,.77549,.78265             !rho0,rho+-,omega
0565      *           ,.889166,.89594,1.019455          !K*+-,K0*,phi
0566      1     ,2.00693,2.01022,2.1123,3.096916        !D0*,D*+-,D*s,j/psi
0567      *     ,5.3251,5.3251,5.4154,6.610,9.46030/    !B*+-,B0*,B*s,B*c,upsilon
0568 c     1/2+ baryon mass table
0569       data ambar0/-1.,.93828,.93957,2*-1.,1.1894,1.1925,1.1974
0570      1     ,1.1156,1.3149,1.3213,3*-1.
0571      $     ,2.453               !15          sigma_c++!
0572      $     ,2.454               !            sigma_c+
0573      $     ,2.452               !            sigma_c0
0574      $     ,2.286               !            lambda_c+
0575      2     ,2.576               !19  1340   !Xi'_c+
0576      $     ,2.578               !20  2340   !Xi'_c0
0577      $     ,2.695               !21  3340   !omegac0
0578      $     ,2.471               !22  3240   !Xi_c0
0579      $     ,2.468               !23  3140   !Xi_c+
0580      $     ,3.55                !24  1440
0581      $     ,3.55                !25  2440
0582      $     ,3.70                !26  3440
0583      $     ,4*-1./
0584 c     3/2+ baryon mass table
0585       data ambar1/1.232,1.232,1.232,1.232,-1.,1.3823,1.3820
0586      1     ,1.3875,-1.,1.5318,1.5350,1.6722,2*-1.
0587      2     ,2.519               !15          sigma_c++
0588      $     ,2.52                !            sigma_c+
0589      $     ,2.517               !            sigma_c0
0590      $     ,-1.
0591      $     ,2.645
0592      $     ,2.644
0593      $     ,2.80
0594      $     ,2*-1.
0595      $     ,3.75
0596      $     ,3.75
0597      3     ,3.90
0598      $     ,4.80
0599      $     ,3*-1./
0600 c     entry
0601       id=idi
0602       amass=0.
0603 ctp060829      if(iabs(id).eq.30)then
0604 ctp060829        amass=amhdibar
0605 ctp060829        return
0606 ctp060829      endif
0607       if(idi.eq.0)then
0608         id=1120                 !for air target
0609       elseif(abs(idi).ge.1000000000)then
0610         goto 500                !nucleus
0611       endif
0612       if(idi.gt.10000)return
0613       call idflav(id,ifl1,ifl2,ifl3,jspin,ind)
0614       if(id.ne.0.and.mod(id,100).eq.0) goto400
0615       if(iabs(ifl1).ge.5.or.iabs(ifl2).ge.5.or.iabs(ifl3).ge.5)
0616      1     goto300
0617       if(ifl2.eq.0) goto200
0618       if(ifl1.eq.0) goto100
0619 c          baryons
0620       ind=ind-109*jspin-36*nmes-nqlep
0621       ind=ind-11
0622       amass=(1-jspin)*ambar0(ind)+jspin*ambar1(ind)
0623       return
0624 c          mesons
0625 100   continue
0626       ind=ind-36*jspin-nqlep
0627       ind=ind-11
0628       amass=(1-jspin)*ammes0(ind)+jspin*ammes1(ind)
0629       return
0630 c          quarks and leptons (+deuteron, triton, alpha, Ks and Kl)
0631 200   continue
0632       amass=amlep(ind)
0633       return
0634 c          b and t particles
0635 300   continue
0636       amass=amlep(iabs(ifl2))+amlep(iabs(ifl3))+1.07+.045*jspin
0637       if(ifl1.ne.0) amass=amass+amlep(iabs(ifl1))
0638       return
0639 c          diquarks
0640 400   amass=amlep(iabs(ifl1))+amlep(iabs(ifl2))
0641       return
0642 c          nuclei
0643 500   nbrpro=mod(abs(id/10000),1000)
0644       nbrneu=mod(abs(id/10),1000)-nbrpro
0645       amass=nbrpro*ambar0(2)+nbrneu*ambar0(3)
0646       return
0647       end
0648
0649 cc-----------------------------------------------------------------------
0650 c      subroutine idmix(ic,jspin,icm,idm)
0651 cc     accounts for flavour mixing
0652 cc-----------------------------------------------------------------------
0653 c      parameter (nflav=6)
0654 c      real pmix1(3,2),pmix2(3,2)
0655 c      integer ic(2),icm(2)
0656 c      data pmix1/.25,.25,.5,0.,.5,1./,pmix2/.5,.5,1.,0.,0.,1./
0657 c      icm(1)=0
0658 c      icm(2)=0
0659 c      idm=0
0660 c      i=ic(1)
0661 c      if(i.ne.ic(2))return
0662 c      id=0
0663 c      if(i.eq.100000)id=1
0664 c      if(i.eq. 10000)id=2
0665 c      if(i.eq.  1000)id=3
0666 c      if(id.eq.0)return
0667 c      rnd=rangen()
0668 c      idm=int(pmix1(id,jspin+1)+rnd)+int(pmix2(id,jspin+1)+rnd)+1
0669 c      icm(1)=10**(nflav-idm)
0670 c      icm(2)=ic(1)
0671 c      idm=idm*100+idm*10+jspin
0672 c      return
0673 c      end
0674 c
0675 cc-----------------------------------------------------------------------
0676 c      subroutine idcleanjc(jc)
0677 cc-----------------------------------------------------------------------
0678 c      parameter (nflav=6)
0679 c      integer jc(nflav,2)
0680 c      ns=0
0681 c      do n=1,nflav
0682 c        jj=min(jc(n,1),jc(n,2))
0683 c        jc(n,1)=jc(n,1)-jj
0684 c        jc(n,2)=jc(n,2)-jj
0685 c        ns=ns+jc(n,1)+jc(n,2)
0686 c      enddo
0687 c      if(ns.eq.0)then
0688 c        jc(1,1)=1
0689 c        jc(1,2)=1
0690 c      endif
0691 c      end
0692 c
0693 c-----------------------------------------------------------------------
0694       subroutine idquacjc(jc,nqu,naq)
0695 c     returns quark content of jc
0696 c        jc(nflav,2) = jc-type particle identification code.
0697 c        nqu = # quarks
0698 c        naq = # antiquarks
0699 c-----------------------------------------------------------------------
0700       parameter (nflav=6)
0701       integer jc(nflav,2)
0702       nqu=0
0703       naq=0
0704       do 53 n=1,nflav
0705         nqu=nqu+jc(n,1)
0706 53      naq=naq+jc(n,2)
0707       return
0708       end
0709
0710 c-----------------------------------------------------------------------
0711       subroutine idquacic(ic,nqu,naq)
0712 c     returns quark content of ic
0713 c        ic(2) = ic-type particle identification code.
0714 c        nqu = # quarks
0715 c        naq = # antiquarks
0716 c-----------------------------------------------------------------------
0717       parameter (nflav=6)
0718       integer jc(nflav,2),ic(2)
0719       nqu=0
0720       naq=0
0721       call iddeco(ic,jc)
0722       do 53 n=1,nflav
0723         nqu=nqu+jc(n,1)
0724 53      naq=naq+jc(n,2)
0725       return
0726       end
0727
0728 c-----------------------------------------------------------------------
0729       subroutine idquac(i,nq,ns,na,jc)
0730 c     returns quark content of ptl i from /cptl/ .
0731 c        nq = # quarks - # antiquarks
0732 c        ns = # strange quarks - # strange antiquarks
0733 c        na = # quarks + # antiquarks
0734 c        jc(nflav,2) = jc-type particle identification code.
0735 c-----------------------------------------------------------------------
0736       include 'epos.inc'
0737       include 'epos.incems'
0738       integer jc(nflav,2),ic(2)
0739
0740
0741       if(iabs(idptl(i)).eq.20)then
0742       idptl(i)=230
0743       if(rangen().lt..5)idptl(i)=-230
0744       goto9999
0745       endif
0746
0747       if(iabs(idptl(i)).lt.100)then
0748       nq=0
0749       ns=0
0750       do 1 n=1,nflav
0751       jc(n,1)=0
0752 1     jc(n,2)=0
0753       return
0754       endif
0755 9999  if(idptl(i)/10**8.ne.7)then
0756       call idtr4(idptl(i),ic)
0757       call iddeco(ic,jc)
0758       else
0759       call idtrb(ibptl(1,i),ibptl(2,i),ibptl(3,i),ibptl(4,i),jc)
0760       endif
0761       na=0
0762       nq=0
0763       do 53 n=1,nflav
0764       na=na+jc(n,1)+jc(n,2)
0765 53    nq=nq+jc(n,1)-jc(n,2)
0766       ns=   jc(3,1)-jc(3,2)
0767       return
0768       end
0769
0770 cc-----------------------------------------------------------------------
0771 c      subroutine idquad(i,nq,na,jc)
0772 cc-----------------------------------------------------------------------
0773 cc     returns quark content of ptl i from /cptl/ .
0774 cc        nq = # quarks - # antiquarks
0775 cc        na = # quarks + # antiquarks
0776 cc        jc(nflav,2) = jc-type particle identification code.
0777 cc-----------------------------------------------------------------------
0778 c      include 'epos.inc'
0779 c      integer jc(nflav,2),ic(2)
0780 c
0781 c      id=idptl(i)
0782 c      if(iabs(id).eq.20)then
0783 c      id=230
0784 c      if(rangen().lt..5)id=-230
0785 c      goto9999
0786 c      endif
0787 c
0788 c      if(iabs(id).lt.100)then
0789 c      nq=0
0790 cc      ns=0
0791 c      do 1 n=1,nflav
0792 c      jc(n,1)=0
0793 c1     jc(n,2)=0
0794 c      return
0795 c      endif
0796 c
0797 c9999  if(id/10**8.ne.7)then
0798 c      call idtr4(id,ic)
0799 c      call iddeco(ic,jc)
0800 c      else
0801 c      call idtrb(ibptl(1,i),ibptl(2,i),ibptl(3,i),ibptl(4,i),jc)
0802 c      endif
0803 c      na=0
0804 c      nq=0
0805 c      do 53 n=1,nflav
0806 c      na=na+jc(n,1)+jc(n,2)
0807 c53    nq=nq+jc(n,1)-jc(n,2)
0808 cc      ns=   jc(3,1)-jc(3,2)
0809 c      return
0810 c      end
0811 c
0812 c-----------------------------------------------------------------------
0813       integer function idraflx(proba,xxx,qqs,icl,jc,jcval,j,iso,c)
0814 c-----------------------------------------------------------------------
0815 c     returns random flavor, according to jc and GRV structure function
0816 c for x(=xxx) and Q2(=qqs) for valence quark (jcval) and sea quarks
0817 c and update jc with quark-antiquark cancellation
0818 c             jc : quark content of remnant
0819 c     j=1 quark, j=2 antiquark,
0820 c     iso : isospin
0821 c     proba : probability of the selected quark (output)
0822 c     c     : "v" is to choose a valence quark, "s" a sea or valence quark
0823 c-----------------------------------------------------------------------
0824       include 'epos.inc'
0825       integer jc(nflav,2),jcval(nflav,2)
0826       double precision s,puv,pdv,psv,pcv,pus,pds,pss,pcs,piso,proba
0827       character c*1
0828
0829        if(ish.ge.8)then
0830          write(ifch,10)c,j,xxx,qqs,jc,jcval
0831        endif
0832  10    format('entry idraflx, j,x,q: ',a1,i2,2g13.5,/
0833      &  ,15x,'jc :',2(1x,6i2),/,14x,'jcv :',2(1x,6i2))
0834
0835        puv=dble(jcval(1,j))
0836        pdv=dble(jcval(2,j))
0837        psv=dble(jcval(3,j))
0838        pcv=dble(jcval(4,j))
0839        if(c.eq."v")then
0840          pus=0d0
0841          pds=0d0
0842          pss=0d0
0843          pcs=0d0
0844        else
0845          pus=dble(jc(1,j))-puv
0846          pds=dble(jc(2,j))-pdv
0847          pss=dble(jc(3,j))-psv
0848          pcs=dble(jc(4,j))-pcv
0849        endif
0850
0851       if(ish.ge.8)then
0852         write(ifch,'(a,4f6.3)')'idraflx valence:',puv,pdv,psv,pcv
0853         write(ifch,'(a,4f6.3)')'idraflx sea:',pus,pds,pss,pcs
0854       endif
0855
0856        qq=0.
0857        if(iso.gt.0)then
0858          if(icl.eq.2)puv=puv*0.5d0   !because GRV already take into account the fact that there is 2 u quark in a proton
0859          puv=puv*dble(psdfh4(xxx,qqs,qq,icl,1))
0860          pus=pus*dble(psdfh4(xxx,qqs,qq,icl,-1))
0861          pdv=pdv*dble(psdfh4(xxx,qqs,qq,icl,2))
0862          pds=pds*dble(psdfh4(xxx,qqs,qq,icl,-2))
0863        elseif(iso.lt.0)then
0864          puv=puv*dble(psdfh4(xxx,qqs,qq,icl,2))
0865          pus=pus*dble(psdfh4(xxx,qqs,qq,icl,-2))
0866          if(icl.eq.2)pdv=pdv*0.5d0
0867          pdv=pdv*dble(psdfh4(xxx,qqs,qq,icl,1))
0868          pds=pds*dble(psdfh4(xxx,qqs,qq,icl,-1))
0869        else
0870          piso=(dble(psdfh4(xxx,qqs,qq,icl,1))
0871      &        +dble(psdfh4(xxx,qqs,qq,icl,2)))
0872          if(icl.eq.2)then   !3 quarks
0873            piso=piso/3d0
0874          else               !2 quarks
0875            piso=piso/2d0
0876          endif
0877          puv=puv*piso
0878          pdv=pdv*piso
0879          piso=0.5d0*(dble(psdfh4(xxx,qqs,qq,icl,-1))
0880      &              +dble(psdfh4(xxx,qqs,qq,icl,-2)))
0881          pus=pus*piso
0882          pds=pds*piso
0883        endif
0884        psv=psv*dble(psdfh4(xxx,qqs,qq,icl,3))
0885        pss=pss*dble(psdfh4(xxx,qqs,qq,icl,-3))
0886        if(nrflav.ge.4)then
0887          pcv=pcv*dble(psdfh4(xxx,qqs,qq,icl,4))
0888          pcs=pcs*dble(psdfh4(xxx,qqs,qq,icl,-4))
0889        else
0890          pcv=0d0
0891          pcs=0d0
0892        endif
0893
0894       if(ish.ge.8)then
0895         write(ifch,'(a,4f6.3)')'idraflx P(valence):',puv,pdv,psv,pcv
0896         write(ifch,'(a,4f6.3)')'idraflx P(sea):',pus,pds,pss,pcs
0897       endif
0898
0899       s=puv+pdv+psv+pcv+pus+pds+pss+pcs
0900 c... initialize
0901       i=0.0
0902       if(s.gt.0.)then
0903        r=rangen()*s
0904        if(r.gt.(pdv+pus+pds+pss+psv+pcv+pcs).and.puv.gt.0.)then
0905         i=1
0906         jcval(i,j)=jcval(i,j)-1
0907         proba=puv
0908        elseif(r.gt.(pus+pds+pss+psv+pcv+pcs).and.pdv.gt.0.)then
0909         i=2
0910         jcval(i,j)=jcval(i,j)-1
0911         proba=pdv
0912        elseif(r.gt.(pds+pss+psv+pcv+pcs).and.pus.gt.0.)then
0913         i=1
0914         proba=pus
0915        elseif(r.gt.(pss+psv+pcv+pcs).and.pds.gt.0.)then
0916         i=2
0917         proba=pds
0918        elseif(r.gt.(psv+pcv+pcs).and.pss.gt.0.)then
0919         i=3
0920         proba=pss
0921        elseif(r.gt.(pcv+pcs).and.psv.gt.0.)then
0922         i=3
0923         jcval(i,j)=jcval(i,j)-1
0924         proba=psv
0925        elseif(r.gt.pcs.and.pcv.gt.0.)then
0926         i=4
0927         jcval(i,j)=jcval(i,j)-1
0928         proba=pcv
0929        elseif(pcs.gt.0.)then
0930         i=4
0931         proba=pcs
0932        else
0933         call utstop("Problem in idraflx, should not be !&")
0934        endif
0935       else
0936         i=idrafl(icl,jc,j,"v",0,iretso)      !no update of jc here
0937         if(jc(i,j)-jcval(i,j).lt.1)jcval(i,j)=jcval(i,j)-1
0938         proba=0d0
0939       endif
0940       idraflx=i
0941
0942       if(ish.ge.8)then
0943         write(ifch,'(a,2(1x,6i2))')'jc before updating:',jc
0944         write(ifch,20)i,j,jcval,proba
0945       endif
0946  20   format('i,j|jcval|P:',2i3,' |',2(1x,6i2),' |',g15.3)
0947
0948       call idsufl3(i,j,jc)
0949
0950       if(ish.ge.8)
0951      & write(ifch,'(a,2(1x,6i2))')'jc after updating:',jc
0952
0953       return
0954       end
0955
0956 c-----------------------------------------------------------------------
0957       integer function idrafl(icl,jc,j,c,imod,iretso)
0958 c-----------------------------------------------------------------------
0959 c     returns random flavor,
0960 c     if : c='v' : according to jc
0961 c          c='s' : from sea
0962 c          c='r' : from sea (always without c quarks)
0963 c          c='d' : from sea for second quark in diquark
0964 c          c='c' : take out c quark first
0965 c             jc : quark content of remnant
0966 c     j=1 quark, j=2 antiquark,
0967 c     imod=0     : returns random flavor of a quark
0968 c     imod=1     : returns random flavor of a quark and update jc
0969 c                 (with quark-antiquark cancellation)
0970 c     imod=2     : returns random flavor of a quark and update jc
0971 c                 (without quark-antiquak cancellation -> accumulate quark)
0972 c     imod=3     : returns random flavor of a quark and update jc with
0973 c                  the corresponding quark-antiquark pair
0974 c                 (without quark-antiquak cancellation -> accumulate quark)
0975 c
0976 c     iretso=0   : ok
0977 c           =1   : more than 9 quarks of same flavor attempted
0978 c-----------------------------------------------------------------------
0979       include 'epos.inc'
0980       integer jc(nflav,2),ic(2)
0981       character c*1
0982
0983 c       write(ifch,*)'entry idrafl, j,imod,c: ',j,imod,' ',c
0984
0985       pui=1.
0986       if(c.eq.'s')then
0987         pu=pui
0988         pd=pui*exp(-pi*difud/fkappa)
0989         ps=pui*exp(-pi*difus/fkappa)
0990         pc=pui*exp(-pi*difuc/fkappa)
0991         pu=rstrau(icl)*pu
0992         pd=rstrad(icl)*pd
0993         ps=rstras(icl)*ps
0994         pc=rstrac(icl)*pc
0995       elseif(c.eq.'d')then
0996         pu=pui*exp(-pi*difuuu/fkappa)
0997         pd=pui*exp(-pi*difudd/fkappa)
0998         ps=pui*exp(-pi*difuss/fkappa)
0999         pc=pui*exp(-pi*difucc/fkappa)
1000         pu=pu*rstrau(icl)
1001         pd=pd*rstrad(icl)
1002         ps=ps*rstras(icl)
1003         pc=pc*rstrac(icl)
1004       elseif(c.eq.'v')then
1005         pu=float(jc(1,j))
1006         pd=float(jc(2,j))
1007         ps=float(jc(3,j))
1008         pc=float(jc(4,j))
1009       elseif(c.eq.'r')then
1010         pu=1.
1011         pd=1.
1012         ps=1.
1013         pc=0.
1014         pu=rstrau(icl)*pu
1015         pd=rstrad(icl)*pd
1016         ps=rstras(icl)*ps
1017       elseif(c.eq.'c')then
1018         pu=0.
1019         pd=0.
1020         ps=0.
1021         pc=1.
1022       else
1023         stop'idrafl: dunnowhatodo'
1024       endif
1025
1026 c      write(ifch,*)'idrafl',pu,pd,ps
1027
1028       s=pu+pd+ps+pc
1029       if(s.gt.0.)then
1030        r=rangen()*s
1031        if(r.gt.(pu+pd+ps).and.pc.gt.0d0)then
1032         i=4
1033        elseif(r.gt.(pu+pd).and.ps.gt.0d0)then
1034         i=3
1035        elseif(r.gt.pu.and.pd.gt.0d0)then
1036         i=2
1037        else
1038         i=1
1039        endif
1040       elseif(iremn.le.1.or.c.ne.'v')then
1041         i=1+min(2,int((2.+rstras(icl))*rangen()))
1042       else
1043         idrafl=0
1044         return
1045       endif
1046       idrafl=i
1047
1048 c      write(ifch,*)'jc before updating',jc
1049 c      write(ifch,*)'i,j,jc',i,j,jc
1050
1051       if(imod.eq.1)then
1052         if(iLHC.eq.0.and.iremn.eq.2)then
1053           call idsufl3(i,j,jc)
1054 c   be sure that jc is not empty
1055           if(jc(i,j).eq.0)then
1056             call idenco(jc,ic,iret)
1057             if(iret.eq.0.and.ic(1).eq.0.and.ic(2).eq.0)then
1058               jc(i,j)=jc(i,j)+1
1059               jc(i,3-j)=jc(i,3-j)+1
1060               iretso=1
1061             endif
1062           endif
1063         elseif(iremn.ge.2)then
1064           call idsufl3(i,j,jc)
1065         else
1066           call idsufl(i,j,jc,iretso)
1067           if(iretso.ne.0.and.ish.ge.2)then
1068             call utmsg('idrafl')
1069             write(ifmt,*)'iret none 0 in idrafl',iretso
1070             write(ifch,*)'iret none 0 in idrafl',iretso
1071             call utmsgf
1072           endif
1073         endif
1074       elseif(imod.eq.2)then
1075         call idsufl2(i,j,jc)    !do not cancel antiquarks with quarks
1076       elseif(imod.eq.3)then
1077         call idsufl2(i,1,jc)    !do not cancel antiquarks with quarks
1078         call idsufl2(i,2,jc)    !do not cancel antiquarks with quarks
1079       endif
1080
1081
1082 c      write(ifch,*)'jc after updating',jc
1083
1084       return
1085       end
1086
1087
1088 c-----------------------------------------------------------------------
1089       integer function idraflz(jc,j)
1090 c-----------------------------------------------------------------------
1091       include 'epos.inc'
1092       integer jc(nflav,2)
1093
1094       pu=float(jc(1,j))
1095       pd=float(jc(2,j))
1096       ps=float(jc(3,j))
1097       pc=float(jc(4,j))
1098
1099       s=pu+pd+ps+pc
1100       if(s.gt.0.)then
1101        r=rangen()*s
1102        if(r.gt.(pu+pd+ps).and.pc.gt.0d0)then
1103         i=4
1104        elseif(r.gt.(pu+pd).and.ps.gt.0d0)then
1105         i=3
1106        elseif(r.gt.pu.and.pd.gt.0d0)then
1107         i=2
1108        else
1109         i=1
1110        endif
1111       else
1112        stop'in idraflz (1)                      '
1113       endif
1114       idraflz=i
1115
1116       if(jc(i,j).lt.1)stop'in idraflz (2)              '
1117       jc(i,j)=jc(i,j)-1
1118
1119       end
1120
1121 c-----------------------------------------------------------------------
1122       subroutine idsufl(i,j,jc,iretso)
1123 c-----------------------------------------------------------------------
1124 c subtract flavor i, j=1 quark, j=2 antiquark
1125 c add antiflavor if jc(i,j)=0
1126 c iretso=0  ok
1127 c       =1 : more than 9 quarks of same flavor attempted
1128 c-----------------------------------------------------------------------
1129       integer jc(6,2),ic(2)
1130
1131       if(jc(i,j).gt.0)then
1132        jc(i,j)=jc(i,j)-1
1133        call idenco(jc,ic,iret)
1134        if(ic(1).eq.0.and.ic(2).eq.0)then
1135          jc(i,j)=jc(i,j)+1
1136          if(jc(i,3-j).lt.9.and.iret.eq.0)then
1137            jc(i,3-j)=jc(i,3-j)+1
1138          else
1139            iretso=1
1140          endif
1141        endif
1142       else
1143         if(j.eq.1)then
1144           if(jc(i,2).lt.9)then
1145             jc(i,2)=jc(i,2)+1
1146           else
1147             iretso=1
1148           endif
1149         else
1150           if(jc(i,1).lt.9)then
1151             jc(i,1)=jc(i,1)+1
1152           else
1153             iretso=1
1154           endif
1155         endif
1156       endif
1157
1158       return
1159       end
1160
1161 c-----------------------------------------------------------------------
1162       subroutine idsufl2(i,j,jc)
1163 c-----------------------------------------------------------------------
1164 c substract flavor i, by adding antiquark i, j=1 quark, j=2 antiquark
1165 c Can replace idsufl if we don't want to cancel quarks and antiquarks
1166 c
1167 c Warning : No protection against jc(i,j)>9 ! should not encode jc without test
1168 c
1169 c-----------------------------------------------------------------------
1170       parameter(nflav=6)
1171       integer jc(nflav,2)
1172
1173       jc(i,3-j)=jc(i,3-j)+1
1174
1175       return
1176       end
1177
1178 c-----------------------------------------------------------------------
1179       subroutine idsufl3(i,j,jc)
1180 c-----------------------------------------------------------------------
1181 c subtract flavor i, j=1 quark, j=2 antiquark
1182 c add antiflavor if jc(i,j)=0
1183 c
1184 c Warning : No protection against jc(i,j)>9 ! should not encode jc without test
1185 c
1186 c-----------------------------------------------------------------------
1187       parameter(nflav=6)
1188       integer jc(nflav,2)
1189
1190       if(jc(i,j).gt.0)then
1191         jc(i,j)=jc(i,j)-1
1192       else
1193         jc(i,3-j)=jc(i,3-j)+1
1194       endif
1195
1196       return
1197       end
1198
1199 cc-----------------------------------------------------------------------
1200 c      subroutine idchfl(jc1,jc2,iret)
1201 cc-----------------------------------------------------------------------
1202 cc checks whether jc1 and jc2 have the same number of quarks and antiquarks
1203 cc if yes: iret=0, if no: iret=1
1204 cc-----------------------------------------------------------------------
1205 c      integer jc1(6,2),jc2(6,2)
1206 c
1207 c      iret=0
1208 c
1209 c      do j=1,2
1210 c       n1=0
1211 c       n2=0
1212 c       do i=1,6
1213 c        n1=n1+jc1(i,j)
1214 c        n2=n2+jc2(i,j)
1215 c       enddo
1216 c       if(n1.ne.n2)then
1217 c        iret=1
1218 c        return
1219 c       endif
1220 c      enddo
1221 c
1222 c      end
1223 c
1224 c
1225 c-----------------------------------------------------------------------
1226       subroutine idres(idi,am,idr,iadj)
1227 c     returns resonance id idr corresponding to mass am.
1228 c     performs mass adjustment, if necessary (if so iadj=1, 0 else)
1229 c     (only for mesons and baryons, error (stop) otherwise)
1230 c-----------------------------------------------------------------------
1231       include 'epos.inc'
1232       parameter (mxindx=1000,mxre=100,mxma=11,mxmx=6)
1233       common/crema/indx(mxindx),rema(mxre,mxma),rewi(mxre,mxma)
1234      *,idmx(mxma,mxmx),icre1(mxre,mxma),icre2(mxre,mxma)
1235       character cad*10
1236
1237       write(cad,'(i10)')idi
1238       iadj=0
1239       idr=0
1240       if(abs(idi).lt.20)then
1241         idr=idi
1242         return
1243       endif
1244       if(abs(am).lt.1.e-5)am=1e-5
1245       id=idi
1246       ami=am
1247       if(am.lt.0.)then
1248         call idmass(id,am)
1249         iadj=1
1250         if(am.le.0.)then
1251         write(ifch,*)'*****  warning in idres (0): '
1252      *,'neg mass returned from idmass'
1253         write(ifch,*)'id,am(input):',idi,ami
1254         am=1e-5
1255         endif
1256       endif
1257
1258       if(id.eq.0)goto 9999
1259       if(abs(id).eq.20)id=sign(230,idi)
1260       m1=1
1261       if(iabs(id).ge.1000)m1=3
1262       m2=2
1263       if(iabs(id).ge.1000)m2=mxmx
1264       do 3 k=m1,m2
1265       do 3 m=2,mxma
1266         if(iabs(id).eq.idmx(m,k)) then
1267           id=idmx(1,k)*10*id/iabs(id)
1268           goto 43
1269         endif
1270  3    continue
1271  43   continue
1272       ix=iabs(id)/10
1273       if(ix.lt.1.or.ix.gt.mxindx)then
1274         call utstop('idres: ix out of range. id='//cad//'&')
1275       endif
1276       i=indx(ix)
1277       if(i.lt.1.or.i.gt.mxre)then
1278         write(ifch,*)'idres problem',id,am
1279         call utstop('idres: particle not in table&')
1280       endif
1281       do 1 j=1,mxma-1
1282       if(am.ge.rema(i,j).and.am.le.rema(i,j+1))then
1283       if(j-1.gt.9)call utstop('idres: spin > 9&')
1284       idr=id/10*10+(j-1)*id/iabs(id)
1285       goto 2
1286       endif
1287 1     continue
1288       goto 9999
1289 2     continue
1290
1291       do 4 k=1,mxmx
1292       if(ix.eq.idmx(1,k))then
1293       if(j.lt.1.or.j.gt.mxma-1)
1294      *call utstop('idres: index j out of range&')
1295       if(idmx(j+1,k).ne.0)idr=idmx(j+1,k)*id/iabs(id)
1296       endif
1297 4     continue
1298
1299       iy=mod(iabs(idr),10)
1300       if(iy.gt.maxres)then
1301       iadj=0
1302       idr=0
1303       goto 9999
1304       endif
1305
1306       if(iy.ne.0.and.iy.ne.1)goto 9999
1307
1308       call idmass(idr,am)
1309       if(am.lt.0.)then
1310       write(ifch,*)'*****  error in idres: '
1311      *,'neg mass returned from idmass'
1312       write(ifch,*)'id,am(input):',idi,ami
1313       write(ifch,*)'idr,am:',idr,am
1314       call utstop('idres: neg mass returned from idmass&')
1315       endif
1316       del=max(1.e-3,2.*rewi(i,j))
1317       if(abs(ami-am).gt.del)iadj=1
1318 c      write(ifch,*)'res:',id,idr,ami,am,rewi(i,j),iadj
1319
1320 9999  if(.not.(ish.ge.8))return
1321       write(ifch,*)'return from idres. id,ami,am,idr,iadj:'
1322       write(ifch,*)idi,ami,am,idr,iadj
1323       return
1324       end
1325
1326 c-----------------------------------------------------------------------
1327       subroutine idresi
1328 c-----------------------------------------------------------------------
1329 c  initializes /crema/
1330 c  masses are limit between stable state (so the average between 2 mass states)
1331 c  width=hbar(6.582e-25 GeV.s)/tau for 151, 251, 351, 451 arbitrary
1332 c  (no data found) !!!!!!!!!!!
1333 c-----------------------------------------------------------------------
1334
1335       parameter (mxindx=1000,mxre=100,mxma=11,mxmx=6)
1336       common/crema/indx(mxindx),rema(mxre,mxma),rewi(mxre,mxma)
1337      *,idmx(mxma,mxmx),icre1(mxre,mxma),icre2(mxre,mxma)
1338       parameter (n=35)
1339       dimension remai(n,mxma),rewii(n,mxma),idmxi(mxma,mxmx)
1340      *,icrei(n,2*mxma)
1341
1342       data (idmxi(j,1),j=1,mxma)/ 11, 110, 111,   0,   0,   0,   0, 4*0/
1343       data (idmxi(j,2),j=1,mxma)/ 22, 220, 330, 331,   0,   0,   0, 4*0/
1344       data (idmxi(j,3),j=1,mxma)/123,2130,1230,1231,1233,1234,1235, 4*0/
1345       data (idmxi(j,4),j=1,mxma)/124,2140,1240,1241,   0,   0,   0, 4*0/
1346       data (idmxi(j,5),j=1,mxma)/134,3140,1340,1341,   0,   0,   0, 4*0/
1347       data (idmxi(j,6),j=1,mxma)/234,3240,2340,2341,   0,   0,   0, 4*0/
1348
1349       data ((icrei(k,m),m=1,2*mxma),k=1,10)/
1350      *111,000000, 9*300000,    11*0,
1351      *222,000000, 9*030000,    11*0,
1352      *112,       10*210000,    11*0,
1353      *122,       10*120000,    11*0,
1354      *113,       10*201000,    11*0,
1355      *223,       10*021000,    11*0,
1356      *123,       10*111000,    11*0,
1357      *133,       10*102000,    11*0,
1358      *233,       10*012000,    11*0,
1359      *333,000000, 9*003000,    11*0/
1360       data ((icrei(k,m),m=1,2*mxma),k=11,20)/
1361      *114,       10*200100,    11*0,
1362      *124,       10*110100,    11*0,
1363      *224,       10*020100,    11*0,
1364      *134,       10*101100,    11*0,
1365      *234,       10*011100,    11*0,
1366      *334,       10*002100,    11*0,
1367      *144,       10*100200,    11*0,
1368      *244,       10*010200,    11*0,
1369      *344,       10*001200,    11*0,
1370      *444,000000, 9*000300,    11*0/
1371       data ((icrei(k,m),m=1,2*mxma),k=21,29)/
1372      * 11,  10*100000,    0,   10*100000,
1373      * 22,  10*001000,    0,   10*001000,
1374      * 12,  10*100000,    0,   10*010000,
1375      * 13,  10*100000,    0,   10*001000,
1376      * 23,  10*010000,    0,   10*001000,
1377      * 14,  10*100000,    0,   10*000100,
1378      * 24,  10*010000,    0,   10*000100,
1379      * 34,  10*001000,    0,   10*000100,
1380      * 44,  10*000100,    0,   10*000100/
1381       data ((icrei(k,m),m=1,2*mxma),k=30,35)/
1382      * 15,  10*100000,    0,   10*000010,
1383      * 25,  10*010000,    0,   10*000010,
1384      * 35,  10*001000,    0,   10*000010,
1385      * 45,  10*000100,    0,   10*000010,
1386      * 55,  10*000010,    0,   10*000010,
1387      *  3,  10*222000,    0,   10*000010/
1388
1389       data ((remai(k,m),m=1,mxma),k=1,10)/
1390      *111.,0.000,1.425,1.660,1.825,2.000,0.000,0.000,0.000,0.000,0.000,
1391      *222.,0.000,1.425,1.660,1.825,2.000,0.000,0.000,0.000,0.000,0.000,
1392      *112.,1.080,1.315,1.485,1.575,1.645,1.685,1.705,1.825,2.000,0.000,
1393      *122.,1.080,1.315,1.485,1.575,1.645,1.685,1.705,1.825,2.000,0.000,
1394      *113.,1.300,1.500,1.700,1.850,2.000,0.000,0.000,0.000,0.000,0.000,
1395      *223.,1.300,1.500,1.700,1.850,2.000,0.000,0.000,0.000,0.000,0.000,
1396      *123.,1.117,1.300,1.395,1.465,1.540,1.655,1.710,1.800,1.885,2.000,
1397 c     *123.,1.154,1.288,1.395,1.463,1.560,1.630,1.710,1.800,1.885,2.000,
1398      *133.,1.423,2.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1399      *233.,1.428,2.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1400 c     *133.,1.423,1.638,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1401 c     *233.,1.427,1.634,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1402      *333.,0.000,2.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000/
1403       data ((remai(k,m),m=1,mxma),k=11,20)/
1404      *114.,2.530,2.730,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1405      *124.,2.345,2.530,2.730,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1406      *224.,2.530,2.730,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1407      *134.,2.450,2.600,2.800,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1408      *234.,2.450,2.600,2.800,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1409      *334.,2.700,2.900,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1410      *144.,3.650,3.850,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1411      *244.,3.650,3.850,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1412      *344.,3.800,4.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1413      *444.,0.000,5.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000/
1414       data ((remai(k,m),m=1,mxma),k=21,29)/
1415      * 11.,0.450,0.950,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1416      * 22.,0.750,0.965,1.500,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1417      * 12.,0.450,0.950,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1418      * 13.,0.500,1.075,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1419      * 23.,0.500,1.075,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1420      * 14.,1.900,2.250,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1421      * 24.,1.900,2.250,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1422      * 34.,2.050,2.500,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1423      * 44.,3.037,3.158,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000/
1424       data ((remai(k,m),m=1,mxma),k=30,35)/
1425      * 15.,5.300,5.400,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1426      * 25.,5.300,5.400,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1427      * 35.,5.396,5.500,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1428      * 45.,6.450,7.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1429      * 55.,9.660,9.999,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,
1430      *  3.,2.230,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000/
1431
1432       data ((rewii(k,m),m=1,mxma),k=1,5)/
1433      *111.,0.000e+00,0.115e+00,0.140e+00,0.250e+00,0.250e+00,
1434      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1435      *222.,0.000e+00,0.115e+00,0.140e+00,0.250e+00,0.250e+00,
1436      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1437      *112.,0.000e+00,0.115e+00,0.200e+00,0.140e+00,0.140e+00,
1438      *     0.145e+00,0.250e+00,0.140e+00,0.250e+00,0.000e+00,
1439      *122.,7.451e-28,0.115e+00,0.200e+00,0.140e+00,0.140e+00,
1440      *     0.145e+00,0.250e+00,0.140e+00,0.250e+00,0.000e+00,
1441      *113.,0.824e-14,0.036e+00,0.080e+00,0.100e+00,0.170e+00,
1442      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00/
1443       data ((rewii(k,m),m=1,mxma),k=6,10)/
1444      *223.,0.445e-14,0.039e+00,0.080e+00,0.100e+00,0.170e+00,
1445      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1446      *123.,0.250e-14,0.890e-05,0.036e+00,0.040e+00,0.016e+00,
1447      *     0.090e+00,0.080e+00,0.100e+00,0.145e+00,0.170e+00,
1448      *133.,0.227e-14,0.009e+00,0.000e+00,0.000e+00,0.000e+00,
1449      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1450      *233.,0.400e-14,0.010e+00,0.000e+00,0.000e+00,0.000e+00,
1451      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1452      *333.,0.000e+00,0.800e-14,0.000e+00,0.000e+00,0.000e+00,
1453      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00/
1454       data ((rewii(k,m),m=1,mxma),k=11,15)/
1455      *114.,0.400e-11,0.010e+00,0.000e+00,0.000e+00,0.000e+00,
1456      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1457      *124.,0.400e-11,0.400e-11,0.010e+00,0.000e+00,0.000e+00,
1458      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1459      *224.,0.400e-11,0.010e+00,0.010e+00,0.000e+00,0.000e+00,
1460      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1461      *134.,0.150e-11,0.400e-11,0.010e+00,0.000e+00,0.000e+00,
1462      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1463      *234.,0.150e-11,0.400e-11,0.010e+00,0.000e+00,0.000e+00,
1464      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00/
1465       data ((rewii(k,m),m=1,mxma),k=16,20)/
1466      *334.,0.400e-11,0.010e+00,0.010e+00,0.000e+00,0.000e+00,
1467      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1468      *144.,0.400e-11,0.010e+00,0.010e+00,0.000e+00,0.000e+00,
1469      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1470      *244.,0.400e-11,0.010e+00,0.010e+00,0.000e+00,0.000e+00,
1471      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1472      *344.,0.400e-11,0.010e+00,0.010e+00,0.000e+00,0.000e+00,
1473      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1474      *444.,0.400e-11,0.010e+00,0.010e+00,0.000e+00,0.000e+00,
1475      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00/
1476       data ((rewii(k,m),m=1,mxma),k=21,25)/
1477      * 11.,7.849e-09,0.153e+00,0.057e+00,0.000e+00,0.000e+00,
1478      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1479      * 22.,0.130e-05,0.210e-03,0.034e+00,0.004e+00,0.000e+00,
1480      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1481      * 12.,2.524e-17,0.153e+00,0.057e+00,0.000e+00,0.000e+00,
1482      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1483      * 13.,5.307e-17,0.051e+00,0.000e+00,0.000e+00,0.000e+00,
1484      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1485      * 23.,0.197e-02,0.051e+00,0.000e+00,0.000e+00,0.000e+00,
1486      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00/
1487       data ((rewii(k,m),m=1,mxma),k=26,29)/
1488      * 14.,0.154e-11,0.002e+00,0.000e+00,0.000e+00,0.000e+00,
1489      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1490      * 24.,0.615e-12,0.002e+00,0.000e+00,0.000e+00,0.000e+00,
1491      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1492      * 34.,0.133e-11,0.002e+00,0.000e+00,0.000e+00,0.000e+00,
1493      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1494      * 44.,0.010e+00,0.068e-03,0.000e+00,0.000e+00,0.000e+00,
1495      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00/
1496       data ((rewii(k,m),m=1,mxma),k=30,35)/
1497      * 15.,0.402e-12,0.010e+00,0.000e+00,0.000e+00,0.000e+00,
1498      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1499      * 25.,0.430e-12,0.010e+00,0.000e+00,0.000e+00,0.000e+00,
1500      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1501      * 35.,0.448e-12,0.010e+00,0.000e+00,0.000e+00,0.000e+00,
1502      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1503      * 45.,0.143e-13,0.010e+00,0.000e+00,0.000e+00,0.000e+00,
1504      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1505      * 55.,0.525e-04,0.010e+00,0.000e+00,0.000e+00,0.000e+00,
1506      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1507      *  3.,0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00,
1508      *     0.000e+00,0.000e+00,0.000e+00,0.000e+00,0.000e+00/
1509
1510       do 3 i=1,mxindx
1511 3     indx(i)=0
1512       do 4 k=1,mxre
1513       do 4 m=1,mxma
1514 4     rema(k,m)=0
1515
1516       do 2 j=1,mxma
1517       do 2 i=1,mxmx
1518 2     idmx(j,i)=idmxi(j,i)
1519
1520       ntec=n
1521       if(ntec.gt.mxre)call utstop('idresi: dimension mxre too small&')
1522       do 1 k=1,n
1523       ix=nint(remai(k,1))
1524       ix2=nint(rewii(k,1))
1525       ix3=icrei(k,1)
1526       if(ix.ne.ix2)call utstop('idresi: ix /= ix2&')
1527       if(ix.ne.ix3)call utstop('idresi: ix /= ix3&')
1528       if(ix.lt.1.or.ix.gt.mxindx)
1529      *call utstop('idresi: ix out of range.&')
1530       indx(ix)=k
1531       rema(k,1)=0.
1532       rewi(k,1)=0.
1533       icre1(k,1)=0
1534       icre2(k,1)=0
1535       do 1 m=2,mxma
1536       rema(k,m)=remai(k,m)
1537       rewi(k,m)=rewii(k,m)
1538       icre1(k,m)=icrei(k,m)
1539 1     icre2(k,m)=icrei(k,mxma+m)
1540
1541       indx(33) =indx(22)
1542       indx(213)=indx(123)
1543       indx(214)=indx(124)
1544       indx(314)=indx(134)
1545       indx(324)=indx(234)
1546
1547       return
1548       end
1549
1550 cc-----------------------------------------------------------------------
1551 c      integer function idsgl(ic,gen,cmp)
1552 cc     returns 1 for singlets (qqq or qqbar) 0 else.
1553 cc-----------------------------------------------------------------------
1554 c      parameter (nflav=6)
1555 c      integer ic(2),jcx(nflav,2),icx(2)
1556 c      character gen*6,cmp*6
1557 c      idsgl=0
1558 c      if(cmp.eq.'cmp-ys')then
1559 c      call idcomi(ic,icx)
1560 c      else
1561 c      icx(1)=ic(1)
1562 c      icx(2)=ic(2)
1563 c      endif
1564 c      call iddeco(icx,jcx)
1565 c      nq=0
1566 c      na=0
1567 c      do 1 i=1,nflav
1568 c      nq=nq+jcx(i,1)
1569 c1     na=na+jcx(i,2)
1570 c      if(nq.eq.0.and.na.eq.0)return
1571 c      if(gen.eq.'gen-no')then
1572 c      if(nq.eq.3.and.na.eq.0.or.nq.eq.1.and.na.eq.1
1573 c     *.or.nq.eq.0.and.na.eq.3)idsgl=1
1574 c      elseif(gen.eq.'gen-ys')then
1575 c      if(mod(nq-na,3).eq.0)idsgl=1
1576 c      endif
1577 c      return
1578 c      end
1579 c
1580 c-----------------------------------------------------------------------
1581       subroutine idtau(id,p4,p5,taugm)
1582 c     returns lifetime(c*tau(fm))*gamma for id with energy p4, mass p5
1583 c-----------------------------------------------------------------------
1584       include 'epos.inc'
1585       parameter (mxindx=1000,mxre=100,mxma=11,mxmx=6)
1586       common/crema/indx(mxindx),rema(mxre,mxma),rewi(mxre,mxma)
1587      *,idmx(mxma,mxmx),icre1(mxre,mxma),icre2(mxre,mxma)
1588            if(iabs(id).eq.14)then
1589       wi=.197/658.654e15
1590            elseif(iabs(id).eq.16)then
1591       wi=.197/87.11e9
1592            elseif(id.eq.-20)then
1593       wi=.197/15.34e15
1594            elseif(id.eq.20)then
1595       wi=.197/2.6842e13
1596            elseif((iabs(id).lt.100.and.id.ne.20)
1597      *         .or.id.gt.1e9)then
1598       wi=0
1599            elseif(iabs(id).lt.1e8)then
1600       ix=iabs(id)/10
1601       if(ix.lt.1.or.ix.gt.mxindx)then
1602         write(ifch,*)'id:',id
1603         call utstop('idtau: ix out of range.&')
1604       endif
1605       ii=indx(ix)
1606       jj=mod(iabs(id),10)+2
1607
1608       m1=1
1609       if(iabs(id).ge.1000)m1=3
1610       m2=2
1611       if(iabs(id).ge.1000)m2=mxmx
1612       do 75 imx=m1,m2
1613       do 75 ima=2,mxma
1614       if(iabs(id).eq.idmx(ima,imx))then
1615         jj=ima
1616         goto 75
1617       endif
1618 75    continue
1619       if(ii.lt.1.or.ii.gt.mxre.or.jj.lt.1.or.jj.gt.mxma)then
1620       write(ifch,*)'id,ii,jj:',id,'   ',ii,jj
1621       call utstop('idtau: ii or jj out of range&')
1622       endif
1623       wi=rewi(ii,jj)
1624            else
1625       tauz=taunll
1626 c-c   tauz=amin1(9./p5**2,tauz)
1627 c-c   tauz=amax1(.2,tauz)
1628       wi=.197/tauz
1629            endif
1630       if(wi.eq.0.)then
1631       tau=ainfin
1632       else
1633       tau=.197/wi
1634       endif
1635       if(p5.ne.0.)then
1636       gm=p4/p5
1637       else
1638       gm=ainfin
1639       endif
1640       if(tau.ge.ainfin.or.gm.ge.ainfin)then
1641       taugm=ainfin
1642       else
1643       taugm=tau*gm
1644       endif
1645       return
1646       end
1647
1648 c-----------------------------------------------------------------------
1649       subroutine idtr4(id,ic)
1650 c     transforms generalized paige_id -> werner_id  (for < 4 flv)
1651 c-----------------------------------------------------------------------
1652       include 'epos.inc'
1653       parameter (mxindx=1000,mxre=100,mxma=11,mxmx=6)
1654       common/crema/indx(mxindx),rema(mxre,mxma),rewi(mxre,mxma)
1655      *     ,idmx(mxma,mxmx),icre1(mxre,mxma),icre2(mxre,mxma)
1656       integer ic(2)
1657
1658       ic(1)=000000
1659       ic(2)=000000
1660       if(mod(abs(id),100).eq.99)return !not a particle
1661       if(iabs(id).lt.20)then
1662         if(id.eq.1)then
1663           ic(1)=100000
1664           ic(2)=000000
1665         elseif(id.eq.-1)then
1666           ic(1)=000000
1667           ic(2)=100000
1668         elseif(id.eq.2)then
1669           ic(1)=010000
1670           ic(2)=000000
1671         elseif(id.eq.-2)then
1672           ic(1)=000000
1673           ic(2)=010000
1674         elseif(id.eq.3)then
1675           ic(1)=001000
1676           ic(2)=000000
1677         elseif(id.eq.-3)then
1678           ic(1)=000000
1679           ic(2)=001000
1680         elseif(id.eq.4)then
1681           ic(1)=000100
1682           ic(2)=000000
1683         elseif(id.eq.-4)then
1684           ic(1)=000000
1685           ic(2)=000100
1686         elseif(id.eq.5)then
1687           ic(1)=000010
1688           ic(2)=000000
1689         elseif(id.eq.-5)then
1690           ic(1)=000000
1691           ic(2)=000010
1692         elseif(id.eq.17)then
1693           ic(1)=330000
1694           ic(2)=000000
1695         elseif(id.eq.-17)then
1696           ic(1)=000000
1697           ic(2)=330000
1698         elseif(id.eq.18)then
1699           ic(1)=450000
1700           ic(2)=000000
1701         elseif(id.eq.-18)then
1702           ic(1)=000000
1703           ic(2)=450000
1704         elseif(id.eq.19)then
1705           ic(1)=660000
1706           ic(2)=000000
1707         elseif(id.eq.-19)then
1708           ic(1)=000000
1709           ic(2)=660000
1710         endif
1711         return
1712       endif
1713       if(id.eq.30)then
1714          ic(1)=222000
1715          ic(2)=000000
1716          return
1717       endif
1718       if(iabs(id).lt.1e8)then
1719         ix=iabs(id)/10
1720         if(ix.lt.1.or.ix.gt.mxindx)goto9999
1721         ii=indx(ix)
1722         if(ii.eq.0)goto9998
1723         jj=mod(iabs(id),10)+2
1724         do 27 imx=1,mxmx
1725           do 27 ima=2,mxma
1726             if(iabs(id).eq.idmx(ima,imx))jj=ima
1727  27     continue
1728         if(id.gt.0)then
1729           ic(1)=icre1(ii,jj)
1730           ic(2)=icre2(ii,jj)
1731         else
1732           ic(2)=icre1(ii,jj)
1733           ic(1)=icre2(ii,jj)
1734         endif
1735         if(ic(1).eq.100000.and.ic(2).eq.100000.and.rangen().lt.0.5)
1736      $       then
1737           ic(1)=010000
1738           ic(2)=010000
1739         endif
1740       elseif(mod(id/10**8,10).eq.8)then
1741         ic(1)=mod(id,10**8)/10000*100
1742         ic(2)=mod(id,10**4)*100
1743       elseif(id/10**9.eq.1.and.mod(id,10).eq.0)then   !nuclei
1744         nstr=mod(id,10**8)/10000000
1745         npro=mod(id,10**7)/10000
1746         nneu=mod(id,10**4)/10
1747         ic(1)=(2*npro+nneu)*10**5+(2*nneu+npro)*10**4+nstr*10**3
1748         ic(2)=0
1749       else
1750         write(ifch,*)'***** id: ',id
1751         call utstop('idtr4: unrecognized id&')
1752       endif
1753       return
1754
1755  9998 continue
1756       write(ifch,*)'id: ',id
1757       call utstop('idtr4: indx=0.&')
1758
1759  9999 continue
1760       write(ifch,*)'id: ',id
1761       call utstop('idtr4: ix out of range.&')
1762       end
1763
1764 c-----------------------------------------------------------------------
1765       integer function idtra(ic,ier,ires,imix)
1766 c-----------------------------------------------------------------------
1767 c     tranforms from werner-id to paige-id
1768 c         ier .... error message (1) or not (0) in case of problem
1769 c         ires ... dummy variable, not used  any more !!!!
1770 c         imix ... 1 not supported any more
1771 c                  2 010000 010000 -> 110, 001000 000100 -> 110
1772 c                  3 010000 010000 -> 110, 001000 000100 -> 220
1773 c-----------------------------------------------------------------------
1774       include 'epos.inc'
1775       parameter (nidt=54)
1776       integer idt(3,nidt),ic(2)!,icm(2)
1777       data idt/
1778      * 100000,100000, 110   ,100000,010000, 120   ,010000,010000, 220
1779      *,100000,001000, 130   ,010000,001000, 230   ,001000,001000, 330
1780      *,100000,000100, 140   ,010000,000100, 240   ,001000,000100, 340
1781      *,000100,000100, 440
1782      *,100000,000010, 150   ,010000,000010, 250   ,001000,000010, 350
1783      *,000100,000010, 450   ,000010,000010, 550   ,100000,000000,   1
1784      *,010000,000000,   2   ,001000,000000,   3   ,000100,000000,   4
1785      *,000010,000000,   5
1786      *,200000,000000,1100   ,110000,000000,1200   ,020000,000000,2200
1787      *,101000,000000,1300   ,011000,000000,2300   ,002000,000000,3300
1788      *,100100,000000,1400   ,010100,000000,2400   ,001100,000000,3400
1789      *,000200,000000,4400
1790      *,330000,000000,  17   ,450000,000000,  18   ,660000,000000,  19
1791      *,300000,000000,1111   ,210000,000000,1120   ,120000,000000,1220
1792      *,030000,000000,2221   ,201000,000000,1130   ,111000,000000,1230
1793      *,000001,000000,   6
1794      *,021000,000000,2230   ,102000,000000,1330   ,012000,000000,2330
1795      *,003000,000000,3331   ,200100,000000,1140   ,110100,000000,1240
1796      *,020100,000000,2240   ,101100,000000,1340   ,011100,000000,2340
1797      *,002100,000000,3340
1798      *,100200,000000,1440   ,010200,000000,2440   ,001200,000000,3440
1799      *,000300,000000,4441/
1800
1801       idtra=0
1802       if(ic(1).eq.0.and.ic(2).eq.0)return
1803       i=1
1804       do while(i.le.nidt.and.idtra.eq.0)
1805         if(ic(2).eq.idt(1,i).and.ic(1).eq.idt(2,i))idtra=-idt(3,i)
1806         if(ic(1).eq.idt(1,i).and.ic(2).eq.idt(2,i))idtra=idt(3,i)
1807         i=i+1
1808       enddo
1809       isi=1
1810       if(idtra.ne.0)isi=idtra/iabs(idtra)
1811
1812       jspin=0
1813
1814       if(imix.eq.1)stop'imix=1 no longer supported'
1815       if(imix.eq.2)then
1816       if(idtra.eq.220)idtra=110
1817       if(idtra.eq.330)idtra=110
1818       elseif(imix.eq.3)then
1819       if(idtra.eq.220)idtra=110
1820       if(idtra.eq.330)idtra=220
1821       endif
1822
1823       if(idtra.ne.0)idtra=idtra+jspin*isi
1824
1825       if(idtra.ne.0)return
1826       if(ier.ne.1)return
1827       write(ifch,*)'idtra: ic = ',ic,ires
1828       call utstop('idtra: unknown code&')
1829
1830       entry idtrai(num,id,ier)
1831       idtrai=0
1832       if(iabs(id).eq.20)then
1833         j=5
1834       elseif(iabs(id).eq.110.or.iabs(id).eq.220)then
1835         j=1+2*int(2.*rangen())
1836       else
1837         j=0
1838         do i=1,nidt
1839           if(iabs(id).eq.idt(3,i))then
1840             j=i
1841             goto 2
1842           endif
1843         enddo
1844  2      continue
1845       endif
1846       if(j.ne.0)then
1847         if(id.lt.0)then
1848           idtrai=idt(3-num,j)
1849         else
1850           idtrai=idt(num,j)
1851         endif
1852         return
1853       endif
1854       if(ier.ne.1)return
1855       write(ifch,*)'idtrai: id = ',id
1856       call utstop('idtrai: unknown code&')
1857       end
1858
1859 c-----------------------------------------------------------------------
1860       subroutine idtrb(ib1,ib2,ib3,ib4,jc)
1861 c     id transformation ib -> jc
1862 c-----------------------------------------------------------------------
1863       parameter (nflav=6)
1864       integer jc(nflav,2)
1865       jc(1,1)=ib1/10**4
1866       jc(2,1)=ib2/10**4
1867       jc(3,1)=ib3/10**4
1868       jc(4,1)=ib4/10**4
1869       jc(5,1)=0
1870       jc(6,1)=0
1871       jc(1,2)=mod(ib1,10**4)
1872       jc(2,2)=mod(ib2,10**4)
1873       jc(3,2)=mod(ib3,10**4)
1874       jc(4,2)=mod(ib4,10**4)
1875       jc(5,2)=0
1876       jc(6,2)=0
1877       return
1878       end
1879
1880 c-----------------------------------------------------------------------
1881       subroutine idtrbi(jc,ib1,ib2,ib3,ib4)
1882 c     id transformation jc -> ib
1883 c-----------------------------------------------------------------------
1884       include 'epos.inc'
1885       integer jc(nflav,2)
1886       ib1=jc(1,1)*10**4+jc(1,2)
1887       ib2=jc(2,1)*10**4+jc(2,2)
1888       ib3=jc(3,1)*10**4+jc(3,2)
1889       ib4=jc(4,1)*10**4+jc(4,2)
1890       ib5=jc(5,1)*10**4+jc(5,2)
1891       ib6=jc(6,1)*10**4+jc(6,2)
1892       if(ib5.ne.0.or.ib6.ne.0)then
1893       write(ifch,*)'***** error in idtrbi: bottom or top quarks'
1894       write(ifch,*)'jc:'
1895       write(ifch,*)jc
1896       call utstop('idtrbi: bottom or top quarks&')
1897       endif
1898       return
1899       end
1900
1901 c------------------------------------------------------------------------------
1902       integer function idtrafo(code1,code2,idi)
1903 c------------------------------------------------------------------------------
1904 c.....tranforms id of code1 (=idi) into id of code2 (=idtrafo)
1905 c.....supported codes:
1906 c.....'nxs' = epos
1907 c.....'pdg' = PDG 1996
1908 c.....'qgs' = QGSJet
1909 c.....'ghe' = Gheisha
1910 c.....'sib' = Sibyll
1911 c.....'cor' = Corsika (GEANT)
1912 c.....'flk' = Fluka
1913
1914 C --- ighenex(I)=EPOS CODE CORRESPONDING TO GHEISHA CODE I ---
1915
1916       common /ighnx/ ighenex(35)
1917       data ighenex/
1918      $               10,   11,   -12,    12,   -14,    14,   120,   110,
1919      $             -120,  130,    20,   -20,  -130,  1120, -1120,  1220,
1920      $            -1220, 2130, -2130,  1130,  1230,  2230, -1130, -1230,
1921      $            -2230, 1330,  2330, -1330, -2330,    17,    18,    19,
1922      $            3331, -3331,  30/
1923
1924 C --- DATA STMTS. FOR GEANT/GHEISHA PARTICLE CODE CONVERSIONS ---
1925 C --- KIPART(I)=GHEISHA CODE CORRESPONDING TO GEANT   CODE I ---
1926 C --- IKPART(I)=GEANT   CODE CORRESPONDING TO GHEISHA CODE I ---
1927       DIMENSION        KIPART(48)!,IKPART(35)
1928       DATA KIPART/
1929      $               1,   3,   4,   2,   5,   6,   8,   7,
1930      $               9,  12,  10,  13,  16,  14,  15,  11,
1931      $              35,  18,  20,  21,  22,  26,  27,  33,
1932      $              17,  19,  23,  24,  25,  28,  29,  34,
1933      $              35,  35,  35,  35,  35,  35,  35,  35,
1934      $              35,  35,  35,  35,  30,  31,  32,  35/
1935
1936 c      DATA IKPART/
1937 c     $               1,   4,   2,   3,   5,   6,   8,   7,
1938 c     $               9,  11,  16,  10,  12,  14,  15,  13,
1939 c     $              25,  18,  26,  19,  20,  21,  27,  28,
1940 c     $              29,  22,  23,  30,  31,  45,  46,  47,
1941 c     $              24,  32,  48/
1942       INTEGER          ICFTABL(200),IFCTABL(-6:100)
1943 C  ICTABL CONVERTS CORSIKA PARTICLES INTO FLUKA PARTICLES
1944 C  FIRST TABLE ONLY IF CHARMED PARTICLES CAN BE TREATED
1945       DATA ICFTABL/
1946      *   7,   4,   3,   0,  10,  11,  23,  13,  14,  12,  ! 10
1947      *  15,  16,   8,   1,   2,  19,   0,  17,  21,  22,  ! 20
1948      *  20,  34,  36,  38,   9,  18,  31,  32,  33,  34,  ! 30
1949      *  37,  39,  24,  25, 6*0,
1950      *  0,    0,   0,   0,  -3,  -4,  -6,  -5,   0,   0,  ! 50
1951      *  10*0,
1952      *   0,   0,   0,   0,   0,   5,   6,  27,  28,   0,  ! 70
1953      *  10*0,
1954      *  10*0,
1955      *  10*0,                                             !100
1956      *  10*0,
1957      *   0,   0,   0,   0,   0,  47,  45,  46,  48,  49,  !120
1958      *  50,   0,   0,   0,   0,   0,   0,   0,   0,   0,  !130
1959      *  41,  42,  43,  44,   0,   0,  51,  52,  53,   0,  !140
1960      *   0,   0,  54,  55,  56,   0,   0,   0,  57,  58,  !150
1961      *  59,   0,   0,   0,  60,  61,  62,   0,   0,   0,  !160
1962      *  40*0/
1963 C  IFCTABL CONVERTS FLUKA PARTICLES INTO CORSIKA PARTICLES
1964       DATA IFCTABL/
1965      *                402, 302, 301, 201,   0,   0,   0,
1966      *  14,  15,   3,   2,  66,  67,   1,  13,  25,   5,
1967      *   6,  10,   8,   9,  11,  12,  18,  26,  16,  21,
1968      *  19,  20,   7,  33,  34,   0,  68,  69,   0,   0,
1969      *  27,  28,  29,  22,  30,  23,  31,  24,  32,   0,
1970      * 131, 132, 133, 134, 117, 118, 116, 119, 120, 121,
1971      * 137, 138, 139, 143, 144, 145, 149, 150, 151, 155,
1972      * 156, 157,   0,   0,   36*0/
1973 c-------------------------------------------------------------------------------
1974
1975       character*3 code1,code2
1976       parameter (ncode=5,nidt=338)
1977       integer idt(ncode,nidt)
1978       double precision drangen,dummy
1979
1980 c            nxs|pdg|qgs|cor|sib
1981       data ((idt(i,j),i=1,ncode),j= 1,68)/
1982      *          1,2,99,99,99             !u quark
1983      *     ,    2,1,99,99,99             !d
1984      *     ,    3,3,99,99,99             !s
1985      *     ,    4,4,99,99,99             !c
1986      *     ,    5,5,99,99,99             !b
1987      *     ,    6,6,99,99,99             !t
1988      *     ,   10,22,99,1,1              !gamma
1989      *     ,   9 ,21,99,99,99            !gluon
1990      *     ,   12,11,11,3,3              !e-
1991      *     ,  -12,-11,-11,2,2            !e+
1992      *     ,   11,12,99,66,15            !nu_e-
1993      *     ,  -11,-12,99,67,16           !nu_e+
1994      *     ,   14,13,99,6,5              !mu-
1995      *     ,  -14,-13,99,5,4             !mu+
1996      *     ,   13,14,99,68,17            !nu_mu-
1997      *     ,  -13,-14,99,69,18           !nu_mu+
1998      *     ,   16,15,99,132,19           !tau-
1999      *     ,  -16,-15,99,131,-19         !tau+
2000      *     ,   15,16,99,133,20           !nu_tau-
2001      *     ,  -15,-16,99,134,-20         !nu_tau+
2002      *     ,  110,111,0,7,6              !pi0
2003      *     ,  120,211,1,8,7              !pi+
2004      *     , -120,-211,-1,9,8            !pi-
2005      *     ,  220,221,10,17,23           !eta
2006      *     ,  130,321,4,11,9             !k+
2007      *     , -130,-321,-4,12,10          !k-
2008      *     ,  230,311,5,33,21            !k0
2009      *     , -230,-311,-5,34,22          !k0b
2010      *     ,   20,310,5,16,12            !kshort
2011      *     ,  -20,130,-5,10,11           !klong
2012      *     ,  330,331,99,99,24           !etaprime
2013      *     ,  111,113,19,51,27           !rho0
2014      *     ,  121,213,99,52,25           !rho+
2015      *     , -121,-213,99,53,26          !rho-
2016      *     ,  221,223,99,50,32           !omega
2017      *     ,  131,323,99,63,28           !k*+
2018      *     , -131,-323,99,64,29          !k*-
2019      *     ,  231,313,99,62,30           !k*0
2020      *     , -231,-313,99,65,31          !k*0b
2021      *     ,  331,333,99,99,33           !phi
2022      *     , -140,421,8,116,99           !D0(1.864)
2023      *     ,  140,-421,8,119,99          !D0b(1.864)
2024      *     , -240,411,7,117,99           !D(1.869)+
2025      *     ,  240,-411,7,118,99          !Db(1.869)-
2026      *     , 1120,2212,2,14,13           !proton
2027      *     , 1220,2112,3,13,14           !neutron
2028      *     , 2130,3122,6,18,39           !lambda
2029      *     , 1130,3222,99,19,34          !sigma+
2030      *     , 1230,3212,99,20,35          !sigma0
2031      *     , 2230,3112,99,21,36          !sigma-
2032      *     , 1330,3322,99,22,37          !xi0
2033      *     , 2330,3312,99,23,38          !xi-
2034      *     , 1111,2224,99,54,40          !delta++
2035      *     , 1121,2214,99,55,41          !delta+
2036      *     , 1221,2114,99,56,42          !delta0
2037      *     , 2221,1114,99,57,43          !delta-
2038      *     , 1131,3224,99,99,44          !sigma*+
2039      *     , 1231,3214,99,99,45          !sigma*0
2040      *     , 2231,3114,99,99,46          !sigma*-
2041      *     , 1331, 3324,99,99,47         !xi*0
2042      *     , 2331, 3314,99,99,48         !xi*-
2043      *     , 3331, 3334,99,24,49         !omega-
2044      *     , 2140, 4122,9,137,99         !LambdaC(2.285)+
2045      *     ,17,1000010020,99,201,1002            !  Deuteron
2046      *     ,18,1000010030,99,301,1003            !  Triton
2047      *     ,19,1000020040,99,402,1004            !  Alpha
2048      *     ,0,0,99,0,0                  !  Air
2049      *     ,99,99,99,99,99 /             !  unknown
2050       data ((idt(i,j),i=1,ncode),j= 69,89)/
2051      $      -340,431,99,120,99           !  Ds+
2052      $     ,340,-431,99,121,99           !  Ds-
2053      $     ,-241,413,99,124,99           !  D*+
2054      $     ,241,-413,99,125,99           !  D*-
2055      $     ,-141,423,99,123,99           !  D*0
2056      $     ,141,-423,99,126,99           !  D*0b
2057      $     ,-341,433,99,127,99           !  Ds*+
2058      $     ,341,-433,99,128,99           !  Ds*-
2059      $     ,440,441,99,122,99            !  etac
2060      $     ,441,443,99,130,99            !  J/psi
2061      $     ,2240,4112,99,142,99          !  sigmac0
2062      $     ,1240,4212,99,141,99          !  sigmac+
2063      $     ,1140,4222,99,140,99          !  sigmac++
2064      $     ,2241,4114,99,163,99          !  sigma*c0
2065      $     ,1241,4214,99,162,99          !  sigma*c+
2066      $     ,1141,4224,99,161,99          !  sigma*c++
2067      $     ,3240,4132,99,139,99          !  Xic0
2068      $     ,2340,4312,99,144,99          !  Xi'c0
2069      $     ,3140,4232,99,138,99          !  Xic+
2070      $     ,1340,4322,99,143,99          !  Xi'c+
2071      $     ,3340,4332,99,145,99 /        !  omegac0
2072       data ((idt(i,j),i=1,ncode),j= 90,nidt)/
2073      $       1112,32224,99,99,99         !  Delta(1600)++
2074      $     , 1112, 2222,99,99,99         !  Delta(1620)++
2075      $     , 1113,12224,99,99,99         !  Delta(1700)++
2076      $     , 1114,12222,99,99,99         !  Delta(1900)++
2077      $     , 1114, 2226,99,99,99         !  Delta(1905)++
2078      $     , 1114,22222,99,99,99         !  Delta(1910)++
2079      $     , 1114,22224,99,99,99         !  Delta(1920)++
2080      $     , 1114,12226,99,99,99         !  Delta(1930)++
2081      $     , 1114, 2228,99,99,99         !  Delta(1950)++
2082      $     , 2222,31114,99,99,99         !  Delta(1600)-
2083      $     , 2222, 1112,99,99,99         !  Delta(1620)-
2084      $     , 2223,11114,99,99,99         !  Delta(1700)-
2085      $     , 2224,11112,99,99,99         !  Delta(1900)-
2086      $     , 2224, 1116,99,99,99         !  Delta(1905)-
2087      $     , 2224,21112,99,99,99         !  Delta(1910)-
2088      $     ,2224,21114,99,99,99          !  Delta(1920)-
2089      $     ,2224,11116,99,99,99          !  Delta(1930)-
2090      $     ,2224, 1118,99,99,99          !  Delta(1950)-
2091      $     ,1122,12212,99,99,99          !  N(1440)+
2092      $     ,1123, 2124,99,99,99          !  N(1520)+
2093      $     ,1123,22212,99,99,99          !  N(1535)+
2094      $     ,1124,32214,99,99,99          !  Delta(1600)+
2095      $     ,1124, 2122,99,99,99          !  Delta(1620)+
2096      $     ,1125,32212,99,99,99          !  N(1650)+
2097      $     ,1125, 2216,99,99,99          !  N(1675)+
2098      $     ,1125,12216,99,99,99          !  N(1680)+
2099      $     ,1126,12214,99,99,99          !  Delta(1700)+
2100      $     ,1127,22124,99,99,99          !  N(1700)+
2101      $     ,1127,42212,99,99,99          !  N(1710)+
2102      $     ,1127,32124,99,99,99          !  N(1720)+
2103      $     ,1128,12122,99,99,99          !  Delta(1900)+
2104      $     ,1128, 2126,99,99,99          !  Delta(1905)+
2105      $     ,1128,22122,99,99,99          !  Delta(1910)+
2106      $     ,1128,22214,99,99,99          !  Delta(1920)+
2107      $     ,1128,12126,99,99,99          !  Delta(1930)+
2108      $     ,1128, 2218,99,99,99          !  Delta(1950)+
2109      $     ,1222,12112,99,99,99          !  N(1440)0
2110      $     ,1223, 1214,99,99,99          !  N(1520)0
2111      $     ,1223,22112,99,99,99          !  N(1535)0
2112      $     ,1224,32114,99,99,99          !  Delta(1600)0
2113      $     ,1224, 1212,99,99,99          !  Delta(1620)0
2114      $     ,1225,32112,99,99,99          !  N(1650)0
2115      $     ,1225, 2116,99,99,99          !  N(1675)0
2116      $     ,1225,12116,99,99,99          !  N(1680)0
2117      $     ,1226,12114,99,99,99          !  Delta(1700)0
2118      $     ,1227,21214,99,99,99          !  N(1700)0
2119      $     ,1227,42112,99,99,99          !  N(1710)0
2120      $     ,1227,31214,99,99,99          !  N(1720)0
2121      $     ,1228,11212,99,99,99          !  Delta(1900)0
2122      $     ,1228, 1216,99,99,99          !  Delta(1905)0
2123      $     ,1228,21212,99,99,99          !  Delta(1910)0
2124      $     ,1228,22114,99,99,99          !  Delta(1920)0
2125      $     ,1228,11216,99,99,99          !  Delta(1930)0
2126      $     ,1228, 2118,99,99,99          !  Delta(1950)0
2127      $     ,1233,13122,99,99,99          !  Lambda(1405)0
2128      $     ,1234, 3124,99,99,99          !  Lambda(1520)0
2129      $     ,1235,23122,99,99,99          !  Lambda(1600)0
2130      $     ,1235,33122,99,99,99          !  Lambda(1670)0
2131      $     ,1235,13124,99,99,99          !  Lambda(1690)0
2132      $     ,1236,13212,99,99,99          !  Sigma(1660)0
2133      $     ,1236,13214,99,99,99          !  Sigma(1670)0
2134      $     ,1237,23212,99,99,99          !  Sigma(1750)0
2135      $     ,1237, 3216,99,99,99          !  Sigma(1775)0
2136      $     ,1238,43122,99,99,99          !  Lambda(1800)0
2137      $     ,1238,53122,99,99,99          !  Lambda(1810)0
2138      $     ,1238, 3126,99,99,99          !  Lambda(1820)0
2139      $     ,1238,13126,99,99,99          !  Lambda(1830)0
2140      $     ,1238,23124,99,99,99          !  Lambda(1890)0
2141      $     ,1239,13216,99,99,99          !  Sigma(1915)0
2142      $     ,1239,23214,99,99,99          !  Sigma(1940)0
2143      $     ,1132,13222,99,99,99          !  Sigma(1660)+
2144      $     ,1132,13224,99,99,99          !  Sigma(1670)+
2145      $     ,1133,23222,99,99,99          !  Sigma(1750)+
2146      $     ,1133,3226,99,99,99           !  Sigma(1775)+
2147      $     ,1134,13226,99,99,99          !  Sigma(1915)+
2148      $     ,1134,23224,99,99,99          !  Sigma(1940)+
2149      $     ,2232,13112,99,99,99          !  Sigma(1660)-
2150      $     ,2232,13114,99,99,99          !  Sigma(1670)-
2151      $     ,2233,23112,99,99,99          !  Sigma(1750)-
2152      $     ,2233,3116,99,99,99           !  Sigma(1775)-
2153      $     ,2234,13116,99,99,99          !  Sigma(1915)-
2154      $     ,2234,23114,99,99,99          !  Sigma(1940)-
2155      $     ,5,7,99,99,99                 !  quark b'
2156      $     ,6,8,99,99,99                 !  quark t'
2157      $     ,16,17,99,99,99               !  lepton tau'
2158      $     ,15,18,99,99,99               !  lepton nu' tau
2159      $     ,90,23,99,99,99               !  Z0
2160      $     ,80,24,99,99,99               !  W+
2161      $     ,81,25,99,99,99               !  h0
2162      $     ,85,32,99,99,99               !  Z'0
2163      $     ,86,33,99,99,99               !  Z''0
2164      $     ,87,34,99,99,99               !  W'+
2165      $     ,82,35,99,99,99               !  H0
2166      $     ,83,36,99,99,99               !  A0
2167      $     ,84,37,99,99,99               !  H+
2168      $     ,1200,2101,99,99,99           !  diquark ud_0
2169      $     ,2300,3101,99,99,99           !  diquark sd_0
2170      $     ,1300,3201,99,99,99           !  diquark su_0
2171      $     ,2400,4101,99,99,99           !  diquark cd_0
2172      $     ,1400,4201,99,99,99           !  diquark cu_0
2173      $     ,3400,4301,99,99,99           !  diquark cs_0
2174      $     ,2500,5101,99,99,99           !  diquark bd_0
2175      $     ,1500,5201,99,99,99           !  diquark bu_0
2176      $     ,3500,5301,99,99,99           !  diquark bs_0
2177      $     ,4500,5401,99,99,99           !  diquark bc_0
2178      $     ,2200,1103,99,99,99           !  diquark dd_1
2179      $     ,1200,2103,99,99,99           !  diquark ud_1
2180      $     ,1100,2203,99,99,99           !  diquark uu_1
2181      $     ,2300,3103,99,99,99           !  diquark sd_1
2182      $     ,1300,3203,99,99,99           !  diquark su_1
2183      $     ,3300,3303,99,99,99           !  diquark ss_1
2184      $     ,2400,4103,99,99,99           !  diquark cd_1
2185      $     ,1400,4203,99,99,99           !  diquark cu_1
2186      $     ,3400,4303,99,99,99           !  diquark cs_1
2187      $     ,4400,4403,99,99,99           !  diquark cc_1
2188      $     ,2500,5103,99,99,99           !  diquark bd_1
2189      $     ,1500,5203,99,99,99           !  diquark bu_1
2190      $     ,3500,5303,99,99,99           !  diquark bs_1
2191      $     ,4500,5403,99,99,99           !  diquark bc_1
2192      $     ,5500,5503,99,99,99           !  diquark bb_1
2193      $     ,800000091,91,99,99,99        !  parton system in cluster fragmentation  (pythia)
2194      $     ,800000092,92,99,99,99        !  parton system in string fragmentation  (pythia)
2195      $     ,800000093,93,99,99,99        !  parton system in independent system  (pythia)
2196      $     ,800000094,94,99,99,99        !  CMshower (pythia)
2197      $     ,250,511,99,99,99             !  B0
2198      $     ,150,521,99,99,99             !  B+
2199      $     ,350,531,99,99,99             !  B0s+
2200      $     ,450,541,99,99,99             !  Bc+
2201      $     ,251,513,99,99,99             !  B*0
2202      $     ,151,523,99,99,99             !  B*+
2203      $     ,351,533,99,99,99             !  B*0s+
2204      $     ,451,543,99,99,99             !  B*c+
2205      $     ,550,551,99,99,99             !  etab
2206      $     ,551,553,99,99,99             !  Upsilon
2207      $     ,2341,4314,99,99,99           !  Xi*c0
2208      $     ,1341,4324,99,99,99           !  Xi*c+
2209      $     ,3341,4334,99,99,99           !  omega*c0
2210      $     ,2440,4412,99,99,99           !  dcc
2211      $     ,2441,4414,99,99,99           !  dcc*
2212      $     ,1440,4422,99,99,99           !  ucc
2213      $     ,1441,4424,99,99,99           !  ucc*
2214      $     ,3440,4432,99,99,99           !  scc
2215      $     ,3441,4434,99,99,99           !  scc*
2216      $     ,4441,4444,99,99,99           !  ccc*
2217      $     ,2250,5112,99,99,99           !  sigmab-
2218      $     ,2150,5122,99,99,99           !  lambdab0
2219      $     ,3250,5132,99,99,99           !  sdb
2220      $     ,4250,5142,99,99,99           !  cdb
2221      $     ,1250,5212,99,99,99           !  sigmab0
2222      $     ,1150,5222,99,99,99           !  sigmab+
2223      $     ,3150,5232,99,99,99           !  sub
2224      $     ,4150,5242,99,99,99           !  cub
2225      $     ,2350,5312,99,99,99           !  dsb
2226      $     ,1350,5322,99,99,99           !  usb
2227      $     ,3350,5332,99,99,99           !  ssb
2228      $     ,4350,5342,99,99,99           !  csb
2229      $     ,2450,5412,99,99,99           !  dcb
2230      $     ,1450,5422,99,99,99           !  ucb
2231      $     ,3450,5432,99,99,99           !  scb
2232      $     ,4450,5442,99,99,99           !  ccb
2233      $     ,2550,5512,99,99,99           !  dbb
2234      $     ,1550,5522,99,99,99           !  ubb
2235      $     ,3550,5532,99,99,99           !  sbb
2236      $     ,3550,5542,99,99,99           !  scb
2237      $     ,2251,5114,99,99,99           !  sigma*b-
2238      $     ,1251,5214,99,99,99           !  sigma*b0
2239      $     ,1151,5224,99,99,99           !  sigma*b+
2240      $     ,2351,5314,99,99,99           !  dsb*
2241      $     ,1351,5324,99,99,99           !  usb*
2242      $     ,3351,5334,99,99,99           !  ssb*
2243      $     ,2451,5414,99,99,99           !  dcb*
2244      $     ,1451,5424,99,99,99           !  ucb*
2245      $     ,3451,5434,99,99,99           !  scb*
2246      $     ,4451,5444,99,99,99           !  ccb*
2247      $     ,2551,5514,99,99,99           !  dbb*
2248      $     ,1551,5524,99,99,99           !  ubb*
2249      $     ,3551,5534,99,99,99           !  sbb*
2250      $     ,4551,5544,99,99,99           !  cbb*
2251      $     ,5551,5554,99,99,99           !  bbb*
2252      $     ,123,10213,99,99,99           !  b1
2253      $     ,122,10211,99,99,99           !  a0+
2254      $     ,233,10313,99,99,99           !  K0_1
2255      $     ,232,10311,99,99,99           !  K*0_1
2256      $     ,133,10323,99,99,99           !  K+_1
2257      $     ,132,10321,99,99,99           !  K*+_1
2258      $     ,143,10423,99,99,99           !  D0_1
2259      $     ,132,10421,99,99,99           !  D*0_1
2260      $     ,243,10413,99,99,99           !  D+_1
2261      $     ,242,10411,99,99,99           !  D*+_1
2262      $     ,343,10433,99,99,99           !  D+s_1
2263      $     ,342,10431,99,99,99           !  D*0s+_1
2264      $     ,223,10113,99,99,99           !  b_10
2265      $     ,222,10111,99,99,99           !  a_00
2266      $     ,113,10223,99,99,99           !  h_10
2267      $     ,112,10221,99,99,99           !  f_00
2268      $     ,333,10333,99,99,99           !  h'_10
2269      $     ,332,10331,99,99,99           !  f'_00
2270      $     ,443,10443,99,99,99           !  h_1c0
2271      $     ,442,10441,99,99,99           !  Xi_0c0
2272      $     ,444,10443,99,99,99           !  psi'
2273      $     ,253,10513,99,99,99           !  db_10
2274      $     ,252,10511,99,99,99           !  db*_00
2275      $     ,153,10523,99,99,99           !  ub_10
2276      $     ,152,10521,99,99,99           !  ub*_00
2277      $     ,353,10533,99,99,99           !  sb_10
2278      $     ,352,10531,99,99,99           !  sb*_00
2279      $     ,453,10543,99,99,99           !  cb_10
2280      $     ,452,10541,99,99,99           !  cb*_00
2281      $     ,553,10553,99,99,99           !  Upsilon'
2282      $     ,552,10551,99,99,99           !  Upsilon'*
2283      $     ,124,20213,99,99,99           !  a_1+
2284      $     ,125,215,99,99,99             !  a_2+
2285      $     ,234,20313,99,99,99           !  K*0_1
2286      $     ,235,315,99,99,99             !  K*0_2
2287      $     ,134,20323,99,99,99           !  K*+_1
2288      $     ,135,325,99,99,99             !  K*+_2
2289      $     ,144,20423,99,99,99           !  D*_10
2290      $     ,135,425,99,99,99             !  D*_20
2291      $     ,244,20413,99,99,99           !  D*_1+
2292      $     ,245,415,99,99,99             !  D*_2+
2293      $     ,344,20433,99,99,99           !  D*_1s+
2294      $     ,345,435,99,99,99             !  D*_2s+
2295      $     ,224,20113,99,99,99           !  a_10
2296      $     ,225,115,99,99,99             !  a_20
2297      $     ,114,20223,99,99,99           !  f_10
2298      $     ,115,225,99,99,99             !  f_20
2299      $     ,334,20333,99,99,99           !  f'_10
2300      $     ,335,335,99,99,99             !  f'_20
2301      $     ,444,20443,99,99,99           !  Xi_1c0
2302      $     ,445,445,99,99,99             !  Xi_2c0
2303      $     ,254,20513,99,99,99           !  db*_10
2304      $     ,255,515,99,99,99             !  db*_20
2305      $     ,154,20523,99,99,99           !  ub*_10
2306      $     ,155,525,99,99,99             !  ub*_20
2307      $     ,354,20533,99,99,99           !  sb*_10
2308      $     ,355,535,99,99,99             !  sb*_20
2309      $     ,454,20543,99,99,99           !  cb*_10
2310      $     ,455,545,99,99,99             !  cb*_20
2311      $     ,554,20553,99,99,99           !  bb*_10
2312      $     ,555,555,99,99,99             !  bb*_20
2313      $     ,11099,9900110,99,99,99       !  diff pi0 state
2314      $     ,12099,9900210,99,99,99       !  diff pi+ state
2315      $     ,22099,9900220,99,99,99       !  diff omega state
2316      $     ,33099,9900330,99,99,99       !  diff phi state
2317      $     ,44099,9900440,99,99,99       !  diff J/psi state
2318      $     ,112099,9902210,99,99,99      !  diff proton state
2319      $     ,122099,9902110,99,99,99      !  diff neutron state
2320      $     ,800000110,110,99,99,99       !  Reggeon
2321      $     ,800000990,990,99,99,99 /     !  Pomeron
2322
2323
2324 c      print *,'idtrafo',' ',code1,' ',code2,idi
2325
2326       nidtmx=68
2327       id1=idi
2328       if(code1.eq.'nxs')then
2329         i=1
2330       elseif(code1.eq.'pdg')then
2331         i=2
2332       elseif(code1.eq.'qgs')then
2333         i=3
2334         if(id1.eq.-10)id1=19
2335       elseif(code1.eq.'cor')then
2336         i=4
2337       elseif(code1.eq.'sib')then
2338         i=5
2339       elseif(code1.eq.'ghe')then
2340         id1=ighenex(id1)
2341         i=1
2342       elseif(code1.eq.'flk')then
2343         id1=IFCTABL(id1)          !convert to corsika code
2344         i=4
2345       else
2346         stop "unknown code in idtrafo"
2347       endif
2348       if(code2.eq.'nxs')then
2349         j=1
2350         ji=j
2351         if(i.eq.2.and.id1.gt.1000000000)then   !nucleus from PDG
2352           idtrafo=id1
2353           return
2354         elseif(i.eq.4.and.id1.gt.402)then               !nucleus from Corsika
2355           idtrafo=1000000000+mod(id1,100)*10000+(id1/100)*10
2356           return
2357         elseif(i.eq.5.and.id1.gt.1004)then               !nucleus from Sibyll
2358           id1=(id1-1000)
2359           idtrafo=1000000000+id1/2*10000+id1*10
2360           return
2361         elseif(id1.eq.130.and.i.eq.2)then
2362           idtrafo=-20
2363           return
2364         endif
2365         if(i.eq.2) nidtmx=nidt
2366         if(i.eq.4) nidtmx=89
2367       elseif(code2.eq.'pdg')then
2368         j=2
2369         ji=j
2370         if(i.eq.1.and.id1.gt.1000000000)then !nucleus from NEXUS
2371           idtrafo=id1
2372           return
2373         elseif(i.eq.4.and.id1.gt.402)then               !nucleus from Corsika
2374           idtrafo=1000000000+mod(id1,100)*10000+(id1/100)*10
2375           return
2376         elseif(i.eq.5.and.id1.gt.1004)then               !nucleus from Sibyll
2377           id1=(id1-1000)
2378           idtrafo=1000000000+id1/2*10000+id1*10
2379           return
2380         elseif(id1.eq.-20.and.i.eq.1)then
2381           idtrafo=130
2382           return
2383         endif
2384         if(i.eq.1) nidtmx=nidt
2385         if(i.eq.4) nidtmx=89
2386        elseif(code2.eq.'qgs')then
2387         j=3
2388         ji=j
2389       elseif(code2.eq.'cor')then
2390         j=4
2391         ji=j
2392       elseif(code2.eq.'sib')then
2393         j=5
2394         ji=j
2395       elseif(code2.eq.'ghe')then
2396         j=4
2397         ji=6
2398       elseif(code2.eq.'flk')then
2399         j=4
2400         ji=7
2401         if(i.le.2) nidtmx=89
2402        else
2403         stop "unknown code in idtrafo"
2404       endif
2405       if(i.eq.4)then !corsika  id always >0 so convert antiparticles
2406         iadtr=id1
2407         if(iadtr.eq.25)then
2408           id1=-13
2409         elseif(iadtr.eq.15)then
2410           id1=-14
2411         elseif(iadtr.ge.26.and.iadtr.le.32)then
2412           id1=-iadtr+8
2413         elseif(iadtr.ge.58.and.iadtr.le.61)then
2414           id1=-iadtr+4
2415         elseif(iadtr.ge.149.and.iadtr.le.157)then
2416           id1=-iadtr+12
2417         elseif(iadtr.ge.171.and.iadtr.le.173)then
2418           id1=-iadtr+10
2419         endif
2420       endif
2421       iad1=abs(id1)
2422       isi=sign(1,id1)
2423
2424       if(i.ne.j)then
2425       do n=1,nidtmx
2426         if(iad1.eq.abs(idt(i,n)))then
2427           m=1
2428           if(n+m.lt.nidt)then
2429             do while(abs(idt(i,n+m)).eq.iad1)
2430               m=m+1
2431             enddo
2432           endif
2433           mm=0
2434           if(m.gt.1)then
2435             if(m.eq.2.and.idt(i,n)*idt(i,n+1).lt.0)then
2436               if(id1.eq.idt(i,n+1))mm=1
2437               isi=1
2438             else
2439               mm=int(drangen(dummy)*dble(m))
2440             endif
2441           endif
2442           idtrafo=idt(j,n+mm)*isi
2443           if(abs(idtrafo).eq.99)call utstop('New particle not allowed ')
2444           if(idtrafo.lt.0.and.j.eq.4)then           !corsika  id always >0
2445             iadtr=abs(idtrafo)
2446             if(iadtr.eq.13)then
2447               idtrafo=25
2448             elseif(iadtr.eq.14)then
2449               idtrafo=15
2450             elseif(iadtr.ge.18.and.iadtr.le.24)then
2451               idtrafo=iadtr+8
2452             elseif(iadtr.ge.54.and.iadtr.le.57)then
2453               idtrafo=iadtr+4
2454             elseif(iadtr.ge.137.and.iadtr.le.145)then
2455               idtrafo=iadtr+12
2456             elseif(iadtr.ge.161.and.iadtr.le.163)then
2457               idtrafo=iadtr+10
2458             else
2459               idtrafo=iadtr
2460             endif
2461           elseif(idtrafo.eq.19.and.j.eq.3)then
2462             idtrafo=-10
2463           endif
2464           if(j.ne.ji)goto 100
2465           return
2466         endif
2467       enddo
2468       else
2469         idtrafo=id1
2470         if(j.ne.ji)goto 100
2471         return
2472       endif
2473
2474       print *, 'idtrafo: ',code1,' -> ', code2,id1,' not found.   '
2475       stop
2476 c      idtrafocx=0
2477 c      return
2478
2479  100  if(j.eq.4)then            !corsika
2480         if(idtrafo.eq.201)then
2481           idtrafo=45
2482         elseif(idtrafo.eq.301)then
2483           idtrafo=46
2484         elseif(idtrafo.eq.402)then
2485           idtrafo=47
2486         elseif(idtrafo.eq.302)then
2487           idtrafo=48
2488         endif
2489         if(idtrafo.ne.0)then      !air
2490           if(ji.eq.6)then
2491             idtrafo=kipart(idtrafo)
2492           elseif(ji.eq.7)then
2493             idtrafo=ICFTABL(idtrafo)
2494           endif
2495         endif
2496         return
2497       else
2498         call utstop('Should not happen in idtrafo !&')
2499       endif
2500
2501       end
2502