IAP GITLAB

Commit 8a1eb2c3 authored by Tanguy Pierog's avatar Tanguy Pierog

add corrections to avoid warnings during compilation

add table production control (not finalized yet)
revision similar to r3250 for epos/branches/1.99.LHC



git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@3252 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 8fadb57b
......@@ -34,7 +34,7 @@ CRMC<OutputPolicy>::CRMC(const CRMCoptions& cfg)
template<class OutputPolicy>
bool
CRMC<OutputPolicy>::init(const int& ipout)
CRMC<OutputPolicy>::init(const int& iptab, const int& ipout)
{
setbuf(stdout, 0); // set output to unbuffered
......@@ -49,7 +49,7 @@ CRMC<OutputPolicy>::init(const int& ipout)
fCfg.fProjectileId,
fCfg.fTargetId,
fCfg.fHEModel,
ipout,
iptab, ipout,
fCfg.GetOutputFileName().c_str(),
fCfg.fParamFileName.c_str());
......
......@@ -33,7 +33,7 @@ class CRMC : public OutputPolicy {
public:
CRMC(const CRMCoptions& cfg);
bool init(const int&);
bool init(const int&, const int&);
bool run(const int&);
bool finish(const int&);
......
......@@ -80,7 +80,7 @@ bool CRMCinterface::init(int HEmodel)
crmc_init_f_ = (void(*)( const int&, const int&, const double&, const double&,
const int&, const int&, const int&, const int&,
const char*, const char*)) dlsym(fLibrary, "crmc_init_f_");
const int&, const char*, const char*)) dlsym(fLibrary, "crmc_init_f_");
if(crmc_init_f_ == NULL)
{
ostringstream errMsg;
......
......@@ -2,7 +2,7 @@
c 15.01.2009 Simplified Main program and random number generator for epos
subroutine crmc_init_f(iEvent,iSeed,pproj,ptarg,
$ ipart,itarg,model,iout,output,param)
$ ipart,itarg,model,itab,iout,output,param)
***************************************************************
*
......@@ -14,6 +14,7 @@ c 15.01.2009 Simplified Main program and random number generator for epos
* ipart - primary particle type
* itarg - target particle type
* model - HE model switch
* itab - force tables production or stop if missing
* iout - output type
* output - output file name
* param - param file name
......@@ -32,9 +33,16 @@ c Input values
double precision ycm2det
common/boostvars/ycm2det
common/producetab/ producetables !used to link CRMC
logical producetables !with EPOS and QII
c Set parameters to default value
call aaset(0)
c Stop program if missing tables (after aaset)
producetables=.false.
if(itab.eq.1)producetables=.true.
c Calculations of energy of the center-of-mass in the detector frame
call idmass(1120,m2) !target mass = proton
m1=m2 !projectile mass
......@@ -109,7 +117,7 @@ c Output quantities
double precision ycm2det
common/boostvars/ycm2det !is this common block needed?
integer i,k
integer i!,k
c Calculate an inelastic event
call aepos(-1)
......@@ -185,8 +193,10 @@ c-----------------------------------------------------------------------
if(iModel.eq.0)then
call LHCparameters !LHC tune for EPOS
lhct=".lhc"
iadd=4
else
lhct=" "
lhct=""
iadd=0
endif
nfnii=15 ! epos tab file name lenght
......@@ -389,6 +399,8 @@ c PDG
if ( model.eq.1 ) then !model variable has no eposLHC
istmax = 1 !final and mother particles (istmax=1 includes
!mother particles)
else
istmax=0
endif
end
......
......@@ -25,7 +25,7 @@ main(int argc, char **argv)
case CRMCoptions::eROOT:
{
CRMC2ROOT crmc(cfg);
crmc.init(0);
crmc.init(itab,0);
crmc.run(0);
crmc.finish(0);
}
......@@ -35,7 +35,7 @@ main(int argc, char **argv)
case CRMCoptions::eHepMCGZ:
{
CRMC2HepMC crmc(cfg);
crmc.init(0);
crmc.init(itab,0);
crmc.run(0);
crmc.finish(0);
}
......@@ -45,7 +45,7 @@ main(int argc, char **argv)
case CRMCoptions::eLHEGZ:
{
CRMC2LHE crmc(cfg);
crmc.init(1);
crmc.init(itab,1);
crmc.run(1);
crmc.finish(1);
}
......
......@@ -4,7 +4,7 @@
switch fusion on !nuclear effects due to high density
!more realistic but slow (can be switched off)
!set istmax 1 !uncomment to get the full particle list (decayed particles)
!set istmax 0 !uncomment to get only final particles with EPOS
!!Set up particle Decays
!switch decay off !no decay at all
......
......@@ -17,6 +17,8 @@ c-----------------------------------------------------------------------
common/ghecsquel/anquasiel,iquasiel
common/cbincond/nozero,ibmin,ibmax
common/photrans/phoele(4),ebeam
common/producetab/ producetables !used to link CRMC
logical producetables !with EPOS and QII
!gauss weights
data (tgss(2,i),i=1,7)/ .3399810436,.8611363116 ,5*0. /
data (wgss(2,i),i=1,7)/ .6521451549,.3478548451 ,5*0. /
......@@ -85,6 +87,8 @@ c file names and units
hydt='---'
producetables=.true.
c initial seed
c following number should be less than kseq in RMMARD (=2 in EPOS)
......@@ -754,7 +758,7 @@ c other
!or multiquark remnant with valence quark conservation and inelastic remnant as low mass droplet only (2)
!or suppression of multiquark remnant with different string end flavors (3)
ntrymx=10 !try-again parameter
istmax=0 !analyse only istptl <= istmax
istmax=1 !analyse only istptl <= istmax
irdmpr=0 !random sign for projectile if 1
ilprtg=1 !consider leading particle in projectile (1)
!or target (-1) side
......@@ -1361,6 +1365,7 @@ c copy first list in final list to define daughters of beam particles
nhep=0
lrcor0=.true. !link spectator remnants to core only once
lrcore=.false.
do j=1,maproj+matarg
......@@ -2096,6 +2101,7 @@ c PARAMETER (MAXNUP=500)
if(istore.eq.-1.or.iappl.eq.-1)then
ifinp=ifdt
if(ichkfile.eq.0)then
if(iappl.eq.-1)then
inquire(file=fnin(1:nfnin),exist=info)
......@@ -2105,8 +2111,6 @@ c PARAMETER (MAXNUP=500)
else
call utstop('Cannot open file for conversion !&')
endif
else
ifinp=ifdt
endif
endif
......@@ -3837,6 +3841,7 @@ ccc if(nopen.eq.0.and.iprmpt.eq.-1)iprmpt=1
if(linex(ix:jx).eq.'engy')then
call idmass(idproj,amproj)
call idmass(idtarg,amtarg)
xxengy=0.
if(engy.gt.0.)then
xxengy=engy
elseif(ecms.gt.0.)then
......@@ -4439,8 +4444,14 @@ ctp290806 nrow=1+(nel-1)/4
c if(ne.eq.0.and.iprmpt.gt.0)
c * write(ifmt,'(a)')'kinks: icbac1 icbac2 icfor1 icfor2?'
call utword(line,i,j,0)
if(line(i:j).eq.'event')ir=1
if(line(i:j).eq.'particle')ir=2
ir=0
if(line(i:j).eq.'event')then
ir=1
elseif(line(i:j).eq.'particle')then
ir=2
else
call utstop("Wrong definition for record!&")
endif
maxrec(ir)=0
20 call utworn(line,j,ne)
c if(ne.eq.0.and.iprmpt.gt.0)
......@@ -4664,8 +4675,14 @@ c if(line(i:j).eq.'off')iorsce=0
elseif(line(i:j).eq.'write'.or.line(i:j).eq.'writex')then
if(line(i:j).eq.'write')ii=1
if(line(i:j).eq.'writex')ii=2
ii=0
if(line(i:j).eq.'write')then
ii=1
elseif(line(i:j).eq.'writex')then
ii=2
else
call utstop("Wrong definition for write!&")
endif
call utword(line,i,j,0)
idol=0
if(line(i:j).eq.'$')then
......@@ -4795,6 +4812,7 @@ c if(line(i:j).eq.'off')iorsce=0
else
ioint=0
iocontr=0
ncontr=0
if(line(i:j).eq.'int')then
ioint=1
call utword(line,i,j,0)
......@@ -5136,6 +5154,7 @@ c for Air target, set the target nucleus
endif
ntry=0
nptly=0
if(nin.le.1)bimevt=-1
c save statistic at last inelastic event
ntevt0=ntevt
......
......@@ -498,6 +498,8 @@ c-----------------------------------------------------------------------
x=0
av=0
su=0
p1=0.
p2=0.
do k=1,kmx
xb=x
x=nbb(k)/float(nnn)
......@@ -1038,6 +1040,8 @@ c-----------------------------------------------------------------------
real x(n),y(n),z(n)
character ch*1
massnr=0
iii=0
if(ch.eq.'p')then
massnr=maproj
iii=1
......@@ -1157,6 +1161,7 @@ c-----------------------------------------------------------------------
if(idtargin.eq.0)imax=max(imax,40)
do massnr=1,mamxx
dif=0.54
rad=0.
if(massnr.gt.imax.or.massnr.eq.1)then
dif=0
rad=0
......@@ -1185,6 +1190,7 @@ c iii = 1 (proj) or 2 (targ)
c-----------------------------------------------------------------------
include 'epos.inc'
if(model.ne.1)return
massnr=1
if(iii.eq.1)then
massnr=maproj
elseif(iii.eq.2)then
......@@ -1288,6 +1294,7 @@ c-----------------------------------------------------------------------
include 'epos.inc'
conrho=1.
if(model.ne.1)return
massnr=1
if(iii.eq.1)then
massnr=maproj
elseif(iii.eq.2)then
......
......@@ -324,6 +324,7 @@ c select one of the decay channel
naddptl=0
sum=0.
c nstart=nptl+1 !?????????????????unused
new=0
do 110 i=1,5 !store id and mass of products
if(mode(i,ipoint).eq.0) goto 110
if(nptl+naddptl+1.gt.mxptl) goto 9999
......@@ -412,6 +413,7 @@ c --------------------------------------
goto1000
endif
rnd(1)=1.
jsave=1
do 310 i=2,naddptl1
rnew=rangen()
i1=i-1
......
......@@ -631,8 +631,11 @@ c-------------------------------------------------------------------------------
if(ntau.eq.1.or.ntau.eq.ntauho)ftau=1
do meta=-netaho+netastep,netaho,netastep
do ii=1,2
if(ii.eq.1)mt=meta-netastep
if(ii.eq.2)mt=meta
if(ii.eq.1)then
mt=meta-netastep
else !if(ii.eq.2)
mt=meta
endif
wwii(ii)=g1*waa(n1,mt,ntau,nphi) +g2*waa(n2,mt,ntau,nphi)
wwii(ii)=max(0.,wwii(ii))
wxii(ii)=g1*vaa(1,n1,mt,ntau,nphi) +g2*vaa(1,n2,mt,ntau,nphi)
......@@ -716,6 +719,7 @@ c-------------------------------------------------------------------------------
pz2p=pz/sqrt(pz**2+ptr**2)
nrap=1+(rap-rapmin)/delrax
npha=1+pha/(2*pi/numiv)
sap=0.501
if(pz2p.le.-1.)then
nsap=0
elseif(pz2p.ge.1.)then
......@@ -1396,6 +1400,7 @@ c----------------------------------------------------------------------
300 continue
endif
else
phinull=0.
do n=1,np
vradi(n)=0.
vtang(n)=0.
......
......@@ -3165,6 +3165,9 @@ c-----------------------------------------------------------------------
jrem=1
elseif(ir.eq.-1)then
jrem=2
else
jrem=0
call utstop("Wrong ir in WriteZZ !&")
endif
do li=1,lremn(irem,jrem)
......@@ -3182,6 +3185,7 @@ c-----------------------------------------------------------------------
c write(ifch,*)'remn',irem,' (',jrem,' ) pom',npom
c & ,' ',zzremn(irem,jrem)
ie=0
is1=0
if(ifrptl(1,npom).gt.0)then
do is=ifrptl(1,npom),ifrptl(2,npom)
if(ie.eq.0)is1=is
......@@ -3278,6 +3282,8 @@ ccc ntrymx=1 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
c initial definitions
ntry=0
iremo1=0
jremo=0
if(ir.eq.1)then
cremn='targ'
jrem=1
......@@ -5147,6 +5153,7 @@ c-----------------------------------------------------------------------
enddo
enddo
flow=1.
if(ir.eq.1)then
c if(kolp(m).le.0)goto1000
if(iep(m).le.-1)goto1000
......@@ -5156,7 +5163,6 @@ c if(kolp(m).le.0)goto1000
zz=zzremn(m,1)
iclpt=iclpro
isopt=isoproj
flow=1.
if(iremn.ge.2)then !number of valence quarks still in proj
if((iept.eq.3.or.iept.eq.5).and.yrmaxi.gt.1.e-5)
& flow=1./fradflii**2
......@@ -5185,7 +5191,6 @@ c if(kolt(m).le.0)goto1000
zz=zzremn(m,2)
iclpt=icltar
isopt=isotarg
flow=1.
if(iremn.ge.2)then !number of valence quarks still in proj
if((iept.eq.3.or.iept.eq.5).and.yrmaxi.gt.1.e-5)
& flow=1./fradflii**2
......@@ -5212,7 +5217,7 @@ c if(kolt(m).le.0)goto1000
&write(ifch,*)'remnant particle index:',mm,m,iclpt,isopt
if(ish.ge.8)call alist('ProRef&',1,nptl)
antotre=antotre+1
antotre=antotre+1.
mmini=mm
nptlini=nptl
......@@ -8380,6 +8385,9 @@ c ------------
if(iokoll.eq.1)then ! precisely matarg collisions
c nothing to do
ntg=0
npj=0
ncoli=0
else
......@@ -8536,6 +8544,8 @@ c if(kolz.eq.0)call utstop(' kolz=0 (proj)&')
xorptl(4,nptl)=t
tivptl(1,nptl)=t
tivptl(2,nptl)=t
naq=0
nqu=0
if(iremn.ge.2)then !update icproj
idp(i)=abs(idp(i))
......@@ -8682,6 +8692,8 @@ c if(kolz.eq.0)call utstop(' kolz=0 (targ)&')
xorptl(4,nptl)=t
tivptl(1,nptl)=t
tivptl(2,nptl)=t
naq=0
nqu=0
if(iremn.ge.2)then !update ictarg
idt(j)=abs(idt(j))
......@@ -8978,6 +8990,7 @@ c Projectile fragment(s)
irest = maproj*100+abs(laproj)
inew=0
idrest=0
mapro=maproj
xmean=0d0
ymean=0d0
......@@ -9913,8 +9926,13 @@ c-----------------------------------------------------------------------
elseif(iii.eq.2)then
do j=1,nid
if(j.eq.1)iclrem=iclpro
if(j.eq.2)iclrem=icltar
if(j.eq.1)then
iclrem=iclpro
elseif(j.eq.2)then
iclrem=icltar
else
iclrem=0
endif
write(ifhi,'(a)') '!----------------------------------'
write(ifhi,'(a)') '! remnant xp distribution '
write(ifhi,'(a)') '!----------------------------------'
......@@ -11326,6 +11344,7 @@ c-----------------------------------------------------------------------
if(xmc.lt.xu)return
if(ptmc.eq.0.)return
ymc=0.
if(iqq.eq.1)ymc=0.5*alog(xmc*s/ptmc)*ih
if(iqq.eq.2)ymc=0.5*alog(xmc/ptmc)
i=1+int(alog(xmc/xu)/alog(xo/xu)*nbix)
......
......@@ -136,6 +136,7 @@ c & p1(5)=dble(pptl(5,iorrem))
c.....pbreak ................................
pbi=pbreak
zz=0.
if(iLHC.eq.1)then
x=p1(5) !sub string mass (energy in cms) for e+e- and pp
if(pbi.gt.0.d0.or.(pbi.gt.-0.99999d0.and.kky.eq.0))then
......@@ -159,7 +160,6 @@ c the ratio Energy_tot / Energy_sub is used to quantify the modification of fkap
x0=x !pomeron energy
x=p1(4) !sub string energy
endif
zz=0.
c if(pbi.gt.0.d0.or.(pbi.gt.-0.99999d0.and.krm.eq.1))then
if(pbi.gt.0.d0.or.(pbi.gt.-0.99999d0.and.kky.eq.0))then
c around 0.4 for soft because of pi0 spectra for CR and pp multiplicity
......@@ -2414,6 +2414,9 @@ C...aplanarity and the related event axes. stolen from jetset ;-)
C...Calculate matrix to be diagonalized.
NP=0
JA=0
JB=0
JC=0
c MSTU41=2
DO 110 J1=1,3
DO 100 J2=J1,3
......@@ -2552,6 +2555,7 @@ C...and the related event axes. stolen from jetset ;-)
dimension thru(4,3)
C...Take copy of particles that are to be considered in thrust analysis.
IAGR=0
NP=0
PS=0.
c MSTU41=2
......@@ -2809,6 +2813,8 @@ c
endif
enddo
iflag=0
i1=0
j1=0
do while (iflag.eq.0.and.nn.ge.2)
a2min=ycut
iflag=1
......
......@@ -3033,6 +3033,7 @@ c----------------------------------------------------------------------
return
endif
x=0.
if(tem.ne.0.)x=chemgc(i)/tem
if(x.le.70.)then
......@@ -3690,6 +3691,8 @@ c --------------------------
c determine pair, construct new pair, update ident
c ------------------------------------------------
xab=1
xba=1
if(iopair.eq.1)then
c (single pair method)
call hnbpad(1,n1,n2,n3,n4,mm,jc)
......@@ -4668,10 +4671,8 @@ c----------------------------------------------------------------------
parameter(maxp=500)
common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
if(k.eq.2)then
k1=n1
k2=n2
endif
c determine n1,n2 and mm
c ----------------------
......@@ -6364,6 +6365,8 @@ c----------------------------------------------------------------------
common/chise/hise(mspecs,nhise)
de=2./nhise/2.
j=0
if(iii.gt.0)then
i=iii
......@@ -6973,6 +6976,7 @@ c-----------------------------------------------------------------------
if(ist.eq.0.and.iostat.eq.1)ist=1
pn=0.0
id=0
jx=100
ymin=1./nevent/10.
......@@ -6998,6 +7002,7 @@ c-----------------------------------------------------------------------
x1=max(x1,0.0)
dx=(x2-x1)/2./jx
x0=x1+dx
pn=0.0
do j=1,jx
datx(j)=x0+(j-1)*dx*2.
......@@ -7115,6 +7120,7 @@ c-----------------------------------------------------------------------
ltyp=nint(xpar6)
ist=nint(xpar7)
if(ist.eq.0.and.iostat.eq.1)ist=1
pn=0.
id=0
......@@ -8041,6 +8047,7 @@ c-----------------------------------------------------------------------
if(iii.lt.0)jjj=nint(xpar1)
id=0
ish0=ish
c ish=98
......@@ -8344,6 +8351,7 @@ c------------------------------------------------------------------------
if(iii.gt.0)stop'STOP: iii>0 not supported'
jjj=nint(xpar1)
id=0
c ----------------
if(iii.eq.0)then
......
......@@ -2157,6 +2157,7 @@ c if(istptl(i).eq.0.or.istptl(i).eq.3)esu=esu+pptl(4,i)
c enddo
c write(ifmt,'(a,41x,f15.1)')' +++++Eajint+++++',esu
delen=0.
if(ish.ge.6)write(ifch,*)'check high pt segments'
if(fploss.gt.0.0)then
ein=0
......@@ -2207,9 +2208,11 @@ c write(ifmt,'(a,41x,f15.1)')' +++++Eajint+++++',esu
va=ya
wa=xa
endif
if(is.eq.-1)imax=1
if(is.eq. 1)imax=m1grid
delen=0
if(is.eq.-1)then
imax=1
else !if(is.eq. 1)
imax=m1grid
endif
vr=sqrt(vx**2+vy**2)
ratx=vz/vr
qu=fploss*delxgr*sqrt(1+rat**2)*sqrt(1+ratx**2)
......@@ -2362,6 +2365,7 @@ c enddo
c...identify clusters
ifirst=0
if(ish.ge.6)write(ifch,*)'identify clusters'
do k=1,m3grid !~~~~~~k-loop
jjj=0
......@@ -3405,6 +3409,7 @@ c & ,pptl(4,nh)
c...ranphi
ranphi=0
rini=0.
if(mjjsegsum.gt.0)then
xx=xx/float(mjjsegsum)
yy=yy/float(mjjsegsum)
......@@ -3621,6 +3626,7 @@ c if(nclu.lt.50)lcont=.false.
ntot=nclu
c prepare debug output for flow
nall=0
if(ish.ge.5)then
nall=nmax-nmin+1
do ii=1,nall
......@@ -3641,6 +3647,8 @@ c initialization for rescaling on ipart particles
ipart=0
tecm0=0.
tecm=0.
ini0=0
ifi0=0
c do 900 while (ncl.le.nptlb-1)
do 900 while (ntot.gt.0)
......@@ -4497,8 +4505,14 @@ ctp to avoid precision problem, replace abs(p3)/p4 by sqrt(1-(pt2+m2)/E2)
xo3=dble(xorptl(3,i))
xo4=dble(xorptl(4,i))
zza=xo3-xo4*vv
if(iopt.eq.0)ti1=dble(tivptl(1,i))
if(iopt.eq.1)ti1=dble(xo4)
if(iopt.eq.0)then
ti1=dble(tivptl(1,i))
elseif(iopt.eq.1)then
ti1=dble(xo4)
else
ti1=0
call utstop("Wrong iopt in jtain !&")
endif
ti2=dble(tivptl(2,i))
if(ti1.gt.ti2)then
......
......@@ -112,15 +112,15 @@ c----------------------------------------------------------------------
! Z_parton_projectile (zparpro) and Z_parton_target (zpartar)-----------
ztav=0.
zpav=0.
if(iscreen.eq.1.or.isplit.eq.1)then
b2x=b2xscr
alpfom=alpfomi*fegypp
bcut=0.
ztav=0.
zpav=0.
rho0p=conrho(1,0.)
rho0t=conrho(2,0.)