From de4b745b472127f5c4f50ff2be9e3ea3c82752c2 Mon Sep 17 00:00:00 2001 From: morsch Date: Mon, 13 Oct 2014 11:11:24 +0200 Subject: [PATCH] TGenerator interface to DIME generator --- DIME/AliDimeRndm.cxx | 81 ++ DIME/AliDimeRndm.h | 44 + DIME/CMakeLists.txt | 10 + DIME/CMakelibTDime.pkg | 37 + DIME/CMakelibdime.pkg | 43 + DIME/DCommon.h | 123 +++ DIME/TDime.cxx | 213 ++++ DIME/TDime.h | 43 + DIME/TDimeLinkDef.h | 8 + DIME/dimeLinkDef.h | 8 + DIME/dimemcv1.05.f | 2320 ++++++++++++++++++++++++++++++++++++++++ 11 files changed, 2930 insertions(+) create mode 100644 DIME/AliDimeRndm.cxx create mode 100644 DIME/AliDimeRndm.h create mode 100644 DIME/CMakeLists.txt create mode 100644 DIME/CMakelibTDime.pkg create mode 100644 DIME/CMakelibdime.pkg create mode 100644 DIME/DCommon.h create mode 100644 DIME/TDime.cxx create mode 100644 DIME/TDime.h create mode 100644 DIME/TDimeLinkDef.h create mode 100644 DIME/dimeLinkDef.h create mode 100644 DIME/dimemcv1.05.f diff --git a/DIME/AliDimeRndm.cxx b/DIME/AliDimeRndm.cxx new file mode 100644 index 00000000000..e34e957ae09 --- /dev/null +++ b/DIME/AliDimeRndm.cxx @@ -0,0 +1,81 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +//----------------------------------------------------------------------------- +// Class: AliDimeRndm +// Responsibilities: Interface to Root random number generator +// from Fortran (re-implements FINCTION RLU_HIJING +// from HIJING) +// Note: Since AliGenDime belongs to another module (TDime) one cannot +// pass the ponter to the generator via static variable +// Collaborators: AliGenDime class +// Example: +// +// root> AliGenDime *gener = new AliGenDime(-1); +// root> AliDimeRndm::SetDimeRandom(new TRandom3()); +// root> AliDimeRndm::GetDimeRandom()->SetSeed(0); +// root> cout<<"Seed "<< AliDimeRndm::GetDimeRandom()->GetSeed() < + +#include "AliDimeRndm.h" + +TRandom * AliDimeRndm::fgDimeRandom=0; + +ClassImp(AliDimeRndm) + +//_______________________________________________________________________ +void AliDimeRndm::SetDimeRandom(TRandom *ran) { + // + // Sets the pointer to an existing random numbers generator + // + if(ran) fgDimeRandom=ran; + else fgDimeRandom=gRandom; +} + +//_______________________________________________________________________ +TRandom * AliDimeRndm::GetDimeRandom() { + // + // Retrieves the pointer to the random numbers generator + // + return fgDimeRandom; +} + +//_______________________________________________________________________ +//# define rluget_hijing rluget_hijing_ +//# define rluset_hijing rluset_hijing_ +//# define rlu_hijing rlu_hijing_ + +extern "C" { + + // void rluget_hijing(Int_t & /*lfn*/, Int_t & /*move*/) + //{printf("Dummy version of rluget_hijing reached\n");} + + //void rluset_hijing(Int_t & /*lfn*/, Int_t & /*move*/) + //{printf("Dummy version of rluset_hijing reached\n");} + + //Float_t rlu_hijing(Int_t & /*idum*/) + // { + // Wrapper to FINCTION RLU_HIJING from HIJING + // Uses static method to retrieve the pointer to the (C++) generator + // Double_t r; + // do r=AliDimeRndm::GetDimeRandom()->Rndm(); + // while(0 >= r || r >= 1); + // return (Float_t)r; + // } +} diff --git a/DIME/AliDimeRndm.h b/DIME/AliDimeRndm.h new file mode 100644 index 00000000000..10a9e5ae50e --- /dev/null +++ b/DIME/AliDimeRndm.h @@ -0,0 +1,44 @@ +#ifndef ALIDIMERNDM_H +#define ALIDIMERNDM_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + +#include +#include + +class TRandom; + +class AliDimeRndm { + public: + AliDimeRndm() { + // Default constructor. The static data member is initialized + // in the implementation file + } + AliDimeRndm(const AliDimeRndm &/*rn*/) { + // Copy constructor: no copy allowed for the object + ::Fatal("Copy constructor","Not allowed\n"); + } + virtual ~AliDimeRndm() { + // Destructor + fgDimeRandom=0; + } + AliDimeRndm & operator=(const AliDimeRndm& /*rn*/) { + // Assignment operator: no assignment allowed + ::Fatal("Assignment operator","Not allowed\n"); + return (*this); + } + + static void SetDimeRandom(TRandom *ran=0); + static TRandom * GetDimeRandom(); + +private: + + static TRandom * fgDimeRandom; //! pointer to the random number generator + + ClassDef(AliDimeRndm,0) //Random Number generator wrapper (non persistent) +}; + +#endif + diff --git a/DIME/CMakeLists.txt b/DIME/CMakeLists.txt new file mode 100644 index 00000000000..ca8d8125701 --- /dev/null +++ b/DIME/CMakeLists.txt @@ -0,0 +1,10 @@ +# AliRoot Build System CMakeLists for HIJING +# +# Author: Johny Jose m(johny.jose@cern.ch) +# Port of previous Makefile build to cmake + +cmake_minimum_required(VERSION 2.8.4 FATAL_ERROR) + +file(GLOB PACKAGES CMake*.pkg) + +ALICE_BuildModule() diff --git a/DIME/CMakelibTDime.pkg b/DIME/CMakelibTDime.pkg new file mode 100644 index 00000000000..c6e219a1d17 --- /dev/null +++ b/DIME/CMakelibTDime.pkg @@ -0,0 +1,37 @@ +# -*- mode: CMake -*- +#--------------------------------------------------------------------------------# +# Package File for THijing # +# Author : Johny Jose (johny.jose@cern.ch) # +# Variables Defined : # +# # +# SRCS - C++ source files # +# HDRS - C++ header files # +# DHDR - ROOT Dictionary Linkdef header file # +# CSRCS - C source files # +# CHDRS - C header files # +# EINCLUDE - Include directories # +# EDEFINE - Compiler definitions # +# ELIBS - Extra libraries to link # +# ELIBSDIR - Extra library directories # +# PACKFFLAGS - Fortran compiler flags for package # +# PACKCXXFLAGS - C++ compiler flags for package # +# PACKCFLAGS - C compiler flags for package # +# PACKSOFLAGS - Shared library linking flags # +# PACKLDFLAGS - Module linker flags # +# PACKBLIBS - Libraries to link (Executables only) # +# EXPORT - Header files to be exported # +# CINTHDRS - Dictionary header files # +# CINTAUTOLINK - Set automatic dictionary generation # +# ARLIBS - Archive Libraries and objects for linking (Executables only) # +# SHLIBS - Shared Libraries and objects for linking (Executables only) # +#--------------------------------------------------------------------------------# + +set ( SRCS TDime.cxx) + +string(REPLACE ".cxx" ".h" HDRS "${SRCS}") + +set ( DHDR TDimeLinkDef.h) + +set ( EXPORT TDime.h) + +set ( EINCLUDE HIJING EVGEN STEER/STEER STEER/STEERBase ) diff --git a/DIME/CMakelibdime.pkg b/DIME/CMakelibdime.pkg new file mode 100644 index 00000000000..171711d471d --- /dev/null +++ b/DIME/CMakelibdime.pkg @@ -0,0 +1,43 @@ +#--------------------------------------------------------------------------------# +# Package File for hijing # +# Author : Johny Jose (johny.jose@cern.ch) # +# Variables Defined : # +# # +# SRCS - C++ source files # +# HDRS - C++ header files # +# DHDR - ROOT Dictionary Linkdef header file # +# CSRCS - C source files # +# CHDRS - C header files # +# EINCLUDE - Include directories # +# EDEFINE - Compiler definitions # +# ELIBS - Extra libraries to link # +# ELIBSDIR - Extra library directories # +# PACKFFLAGS - Fortran compiler flags for package # +# PACKCXXFLAGS - C++ compiler flags for package # +# PACKCFLAGS - C compiler flags for package # +# PACKSOFLAGS - Shared library linking flags # +# PACKLDFLAGS - Module linker flags # +# PACKBLIBS - Libraries to link (Executables only) # +# EXPORT - Header files to be exported # +# CINTHDRS - Dictionary header files # +# CINTAUTOLINK - Set automatic dictionary generation # +# ARLIBS - Archive Libraries and objects for linking (Executables only) # +# SHLIBS - Shared Libraries and objects for linking (Executables only) # +#--------------------------------------------------------------------------------# + +set ( SRCS AliDimeRndm.cxx) + +string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" ) + +set ( DHDR dimeLinkDef.h) + +set ( EXPORT ) + +set ( FSRCS dimemcv1.05.f) + + +if( ALICE_TARGET STREQUAL "macosxicc") + + string(REGEX REPLACE "-O[^ ]*" "" PACKFFLAGS "${FFLAGS}") + +endif( ALICE_TARGET STREQUAL "macosxicc") diff --git a/DIME/DCommon.h b/DIME/DCommon.h new file mode 100644 index 00000000000..0aea28e5b8f --- /dev/null +++ b/DIME/DCommon.h @@ -0,0 +1,123 @@ +#ifndef ROOT_DCommon +#define ROOT_DCommon + +#ifndef __CFORTRAN_LOADED +//*KEEP,cfortran. +#include "cfortran.h" +//*KEND. +#endif + +/* +ccccc hepevt output + integer nmxhep,kk + parameter (nmxhep=4000) + integer nevhep,nhep,isthep,idhep,jmohep,jdahep + double precision phep,vhep + common /hepevt/ nevhep,nhep,isthep(nmxhep),idhep(nmxhep), + &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep) +*/ + const Int_t nmxhep = 4000; + typedef struct { + Int_t nevhep; + Int_t nhep; + Int_t isthep[nmxhep]; + Int_t idhep[nmxhep]; + Int_t jmohep[nmxhep][2]; + Int_t jdahep[nmxhep][2]; + Double_t phep[nmxhep][5]; + Double_t vhep[nmxhep][4]; + } hepevtCommon; + +#define HEPEVT COMMON_BLOCK(HEPEVT, hepevt) +COMMON_BLOCK_DEF(hepevtCommon, HEPEVT); + +/* +ccccc Les Houches Event Common Block + INTEGER MAXNUP + PARAMETER (MAXNUP=500) + INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP + DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP + COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP, + & IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP), + & ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP), + & SPINUP(MAXNUP) + */ + extern "C" { + const Int_t MAXNUP = 500; + typedef struct { + Int_t NUP; + Int_t IDPRUP; + Double_t XWGTUP; + Double_t SCALUP; + Double_t AQEDUP; + Double_t AQCDUP; + Int_t IDUP[MAXNUP]; + Int_t ISTUP[MAXNUP]; + Int_t MOTHUP[MAXNUP][2]; + Int_t ICOLUP[MAXNUP][2]; + Double_t PUP[MAXNUP][5]; + Double_t VTIMUP[MAXNUP]; + Double_t SPINUP[MAXNUP]; + } hepeupCommon; + +#define HEPEUP COMMON_BLOCK(HEPEUP, hepeup) +COMMON_BLOCK_DEF(hepeupCommon, HEPEUP); + + /* + common/vars/s,rts,mmes,yx + */ + typedef struct { + Double_t s; + Double_t rts; + Double_t mmes; + Double_t yx; + } varsCommon; + +#define VARS COMMON_BLOCK(VARS, vars) + COMMON_BLOCK_DEF(varsCommon, VARS); + + + //ccc + // common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut + typedef struct { + Double_t etaelmax; + Double_t etaelmin; + Double_t ptelmin; + Double_t ptphmin; + Double_t ecut; + Double_t rmax; + Double_t rmin; + Double_t mcut; + } cutsCommon; +#define CUTS COMMON_BLOCK(CUTS, cuts) + COMMON_BLOCK_DEF(cutsCommon, CUTS); + + /* + character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10 + &,ppbar*10,output*10,mregge*10,cuts*10,unw*10 + common/flags/ iin, pflag, fsi, ppbar, output, cuts, unw + common/ff/formf + */ + typedef struct { + Int_t iin; + char pflag[10]; + char fsi[10]; + char ppbar[10]; + char output[10]; + char cuts[10]; + char unw[10]; + } flagsCommon; + +#define FLAGS COMMON_BLOCK(FLAGS, flags) + COMMON_BLOCK_DEF(flagsCommon, FLAGS); + + typedef struct { + char formf[10]; + } ffCommon; +#define FF COMMON_BLOCK(FF, ff) + COMMON_BLOCK_DEF(ffCommon, FF); + + +#endif +} + diff --git a/DIME/TDime.cxx b/DIME/TDime.cxx new file mode 100644 index 00000000000..e6c547ed73b --- /dev/null +++ b/DIME/TDime.cxx @@ -0,0 +1,213 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id:*/ + +#include +#include +#include +#include +#include "Dcommon.h" +#include "TDime.h" + +#ifndef WIN32 +# define dimeinit dimeinit_ +# define dimegenerate dimegenerate_ +# define type_of_call +#else +# define dimeinit DIMEINIT +# define dimegenerate DIMEGENERATE +# define type_of_call _stdcall +#endif + +#ifndef WIN32 +extern "C" void type_of_call dimeinit(); +extern "C" void type_of_call dimegenerate(); +#else +#endif + + + +ClassImp(TDime) + + +TDime::TDime(): + TGenerator("Dime","Dime"), + fEfrm(5500.) +{ +// Default constructor +} + +//______________________________________________________________________________ +TDime::TDime(Double_t efrm) : + TGenerator("Dime","Dime"), + fEfrm(efrm) +{ +// TDime constructor: +// Note that there may be only one functional TDime object +// at a time, so it's not use to create more than one instance of it. +} + +//______________________________________________________________________________ +TDime::~TDime() +{ +// Destructor +} + + void TDime::Initialize() +{ + VARS.rts = fEfrm; + CUTS.rmax = 1.8; + CUTS.rmin = -1.8; + CUTS.ecut = 2.; + strcpy(FLAGS.pflag, "rho "); + strcpy(FLAGS.fsi, "true "); + strcpy(FLAGS.ppbar, "false "); + strcpy(FLAGS.cuts, "true "); + strcpy(FLAGS.unw, "true "); + strcpy(FF.formf, "orexp "); + dimeinit(); +} + +void TDime::GenerateEvent() +{ + dimegenerate(); + //for (Int_t i = 0; i < HEPEUP.NUP; i++) { + //printf("%5d %5d %5d %5d %5d %13.3f %13.3f\n", i, + // HEPEUP.IDUP[i], HEPEUP.ISTUP[i], HEPEUP.MOTHUP[i][0], + // HEPEUP.ICOLUP[i][0], HEPEUP.PUP[i][2], HEPEUP.VTIMUP[i]); + //} + +} + +TObjArray* TDime::ImportParticles(Option_t *option) +{ +// +// Default primary creation method. It reads the /HEPEVT/ common block which +// has been filled by the GenerateEvent method. If the event generator does +// not use the HEPEVT common block, This routine has to be overloaded by +// the subclasses. +// The function loops on the generated particles and store them in +// the TClonesArray pointed by the argument particles. +// The default action is to store only the stable particles (ISTHEP = 1) +// This can be demanded explicitly by setting the option = "Final" +// If the option = "All", all the particles are stored. +// + fParticles->Clear(); + Int_t nump = 0; + + Int_t numpart = HEPEUP.NUP; + printf("\n TDime: DIME stack contains %d particles.", numpart); + + if (!strcmp(option,"") || !strcmp(option,"Final")) { + for (Int_t i = 0; i < numpart; i++) { + + if (HEPEUP.ISTUP[i] == 1) { +// +// Use the common block values for the TParticle constructor +// + nump++; + TParticle* p = new TParticle( + HEPEUP.IDUP[i], HEPEUP.MOTHUP[i][0], HEPEUP.MOTHUP[i][1] , + -1, -1, -1, + HEPEUP.PUP[i][0], HEPEUP.PUP[i][1], HEPEUP.PUP[i][2], HEPEUP.PUP[i][3] , + 0., 0., 0., 0.0); + fParticles->Add(p); + } + } + } + else if (!strcmp(option,"All")) { + nump = numpart; + for (Int_t i = 0; i < numpart; i++) { + Int_t iParent = HEPEUP.MOTHUP[i][0]-1; + if (iParent >= 0) { + TParticle *mother = (TParticle*) (fParticles->UncheckedAt(iParent)); + mother->SetLastDaughter(i); + if (mother->GetFirstDaughter()==-1) + mother->SetFirstDaughter(i); + } + + TParticle* p = new TParticle( + HEPEUP.IDUP[i], HEPEUP.MOTHUP[i][0]-1, HEPEUP.MOTHUP[i][1]-1 , + -1, -1, -1, + HEPEUP.PUP[i][0], HEPEUP.PUP[i][1], HEPEUP.PUP[i][2], HEPEUP.PUP[i][3] , + 0., 0., 0., 0.); + fParticles->Add(p); + } + } + return fParticles; +} + +Int_t TDime::ImportParticles(TClonesArray *particles, Option_t *option) +{ +// +// Default primary creation method. It reads the /HEPEVT/ common block which +// has been filled by the GenerateEvent method. If the event generator does +// not use the HEPEVT common block, This routine has to be overloaded by +// the subclasses. +// The function loops on the generated particles and store them in +// the TClonesArray pointed by the argument particles. +// The default action is to store only the stable particles (ISTHEP = 1) +// This can be demanded explicitly by setting the option = "Final" +// If the option = "All", all the particles are stored. +// + if (particles == 0) return 0; + TClonesArray &particlesR = *particles; + particlesR.Clear(); + Int_t nump = 0; + + Int_t numpart = HEPEUP.NUP; + printf("\n TDime: DIME stack contains %d particles.", numpart); + + if (!strcmp(option,"") || !strcmp(option,"Final")) { + for (Int_t i = 0; i < numpart; i++) { + if (HEPEUP.ISTUP[i] == 1) { +// +// Use the common block values for the TParticle constructor +// + nump++; + new(particlesR[i]) + TParticle( + HEPEUP.IDUP[i], HEPEUP.MOTHUP[i][0], HEPEUP.MOTHUP[i][1] , + -1, -1, -1, + HEPEUP.PUP[i][0], HEPEUP.PUP[i][1], HEPEUP.PUP[i][2], HEPEUP.PUP[i][3] , + 0., 0., 0., 0.0); + } + } + } + else if (!strcmp(option,"All")) { + nump = numpart; + for (Int_t i = 0; i < numpart; i++) { + + Int_t iParent = HEPEUP.MOTHUP[i][0]-1; + + if (iParent >= 0) { + TParticle *mother = (TParticle*) (particlesR.UncheckedAt(iParent)); + mother->SetLastDaughter(i); + if (mother->GetFirstDaughter()==-1) + mother->SetFirstDaughter(i); + } + + new(particlesR[i]) TParticle( + HEPEUP.IDUP[i], HEPEUP.MOTHUP[i][0], HEPEUP.MOTHUP[i][1] , + -1, -1, -1, + HEPEUP.PUP[i][0], HEPEUP.PUP[i][1], HEPEUP.PUP[i][2], HEPEUP.PUP[i][3] , + 0., 0., 0., 0.0 + ); + } + } + return nump; +} + diff --git a/DIME/TDime.h b/DIME/TDime.h new file mode 100644 index 00000000000..2afd26b0091 --- /dev/null +++ b/DIME/TDime.h @@ -0,0 +1,43 @@ +#ifndef TDIME_H +#define TDIME_H + + +////////////////////////////////////////////////////////////////////////// +// // +// TDime // +// // +// This class implements an interface to the DIME event generator. // +// // +////////////////////////////////////////////////////////////////////////// + +#ifndef ROOT_TGenerator +#include "TGenerator.h" +#endif +class TObjArray; + +class TDime : public TGenerator { +public: + + TDime(); + TDime(Double_t efrm); + virtual ~TDime(); + virtual void Initialize(); + virtual void GenerateEvent(); + virtual Int_t ImportParticles(TClonesArray *particles, Option_t *option=""); + virtual TObjArray* ImportParticles(Option_t *option=""); + //Parameters for the generation: + virtual void SetEnergyCMS(Float_t efrm) {fEfrm = efrm;} + virtual Float_t GetEnergyCMS() const {return fEfrm;} + protected: + Float_t fEfrm; // Energy in the centre of mass (CMS) or lab-frame (LAB) + ClassDef(TDime,1) //Interface to Dime Event Generator +}; + +#endif + + + + + + + diff --git a/DIME/TDimeLinkDef.h b/DIME/TDimeLinkDef.h new file mode 100644 index 00000000000..a3cfb80385a --- /dev/null +++ b/DIME/TDimeLinkDef.h @@ -0,0 +1,8 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class TDime+; +#endif diff --git a/DIME/dimeLinkDef.h b/DIME/dimeLinkDef.h new file mode 100644 index 00000000000..6ac223ba8da --- /dev/null +++ b/DIME/dimeLinkDef.h @@ -0,0 +1,8 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class AliDimeRndm+; +#endif diff --git a/DIME/dimemcv1.05.f b/DIME/dimemcv1.05.f new file mode 100644 index 00000000000..c2806a959a7 --- /dev/null +++ b/DIME/dimemcv1.05.f @@ -0,0 +1,2320 @@ + subroutine dimeinit + + implicit double precision(a-y) + implicit complex*16(z) + double precision pt1(2),pt2(2),ptx(2),q(4,20) + double precision pboo(4),pcm(4),plb(4) + double precision tvec(4),uvec(4),svec(4),vv1(4),vv2(4) + integer d,h,i,j,k,l,m,n,o,p,r,iset,af,nhist,ntotal,eflag,ichi, + &icut,ncut,id(20),id0,id1,id2,pdf,num,ptgen,mode,iord,ii,ip,ll,nev + &,hh,lmax + character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10 + &,ppbar*10,output*10,mregge*10,cuts*10,unw*10 + integer iin,nch + double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5) + common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm + common/ipars/nch + common/mom/q + common/int/ip + common/mompt/pt1,pt2,ptx + common/vars/s,rts,mmes,yx +ccc new + common/vars0/ mf127, mf1525, m0, mmes0, mmes1, mmes2, mp, + & mwidth, pi, rt2, ebeam, sum, sum1, weightm, + & bjac, bp + common /ivars/ ncut, ll, icut, nev, num +ccc new + common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut + common/flags/iin, pflag, fsi, ppbar, output, cuts, unw + common/gluons/g1,g2,kt1,kt2 + common/levicivita/epsilon + common/ang/cost,costa,costb + common/alph/alfas + common/hvars/sh,th,uh + common/wvars/sig0,bb,t1,t2 + common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m + common/sec/cpom,cf,crho,aff,ar + integer iii + common/ivar/iii + common/ff/formf + common/ffpars/bexp,ao,bo,aoo,boo + common/regge/mregge +ccccc hepevt output + integer nmxhep,kk + parameter (nmxhep=4000) + integer nevhep,nhep,isthep,idhep,jmohep,jdahep + double precision phep,vhep + common /hepevt/ nevhep,nhep,isthep(nmxhep),idhep(nmxhep), + &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep) +ccccc Les Houches Event Common Block + INTEGER MAXNUP + PARAMETER (MAXNUP=500) + INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP + DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP + COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP, + & IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP), + & ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP), + & SPINUP(MAXNUP) +ccccc rambo variables + integer npart,nrun + double precision pin(4),am(100),pout(4,100),wt,ein + common/ramboc/ pin, am, pout, wt, ein, + &npart, nrun +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c c +c Dime MC for the central exclusive production of c +c meson pairs by double Pomeron exchange: c +c c +c p(1) p(2) --> p(3) + M_1(5) M_2(6) + p(4) c +c c +c Momenta for each event in array q(i,j), where j is c +c the particle label and i is the 4-momentum c +c component, with: c +c c +c 1,2 = transverse components c +c 3 = beam component c +c 4 = energy c +c c +c PDG number for ith particle in arrary ID(i) c +c c +c Also gives event record in HEPEVT or LHE format c +c (others are available upon request) c +c c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c c +c Particles generated: c +c c +c pflag='pipm' : M_1=pi^+ M_2=pi^- c +c pflag='pi0' : M_1=M_2=pi^0 c +c pflag='kpkm' : M_1=K^+ M_2=K^- c +c pflag='ks' : M_1=M_2=K_0 c +c pflag='rho' : M_1=M_2=rho_0 c +c c +c with decay: rho(5) --> pi^+(7)pi^-(8) c +c rho(6) --> pi^+(9)pi^-(10) c +c according to phase space, with finite rho c +c width included c +c c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c c +c User defined cuts can readily be implemented c +c in subroutine 'icut' c +c c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c c +c This is version 1.05 : 17 Sep 2014 c +c c +c Comments etc to: lucian.harland-lang@durham.ac.uk c +c c +c If you use the code please reference (and see for c +c details of model used): c +c c +c L.A. Harland-Lang, V.A. Khoze, M.G. Ryskin c +c and W.J. Stirling arXiv:1105.1626 c +c c +c L.A. Harland-Lang, V.A. Khoze, M.G. Ryskin c +c and W.J. Stirling arXiv:1204.4803 c +c c +c L.A. Harland-Lang, V.A. Khoze and M.G. Ryskin c +c arXiv:1312.4553 c +c c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +ccAM open(35,file='exerec1.dat') ! event record + +ccAM rts=1.96d3 ! centre of mass energy + write(6,*) "energy", rts +cccc Some basic cuts, imposed in subtroutine 'icut'. Other user defined cuts can readily be implemented in subroutine +cccc note: in rhorho case these cuts are imposed on rho's, and not their decay productions. Cuts on the decay products +cccc can be imposed by editing 'icut' + +ccAM rmax=1.8d1 ! max meson rapidity +ccAM rmin=-1.8d1 ! min meson rapidity +ccAM ecut=2d0 ! min meson p_t + +c -- physics parameters + +ccAM pflag='rho' ! Process generated - see preamble for options +ccAM fsi='true' ! phenomenological model for exclusive suppression in Pom Pom --> meson pair. To turn on/off --> 'true/false' +ccAM formf='orexp' ! meson - pomeron form factor. +ccAM ppbar='false' ! set true if ppbar collisions + output='lhe' ! Event record style, HEPEVT or LHE +ccAM cuts='true' ! Impose cuts or not +ccAM unw='true' ! Set = 'true' for unweighted events +ccAM iin=1 ! Model for soft survival factor, as described in arXiv:1306.2149. Default = 1 + +cccccccc + +ccAM ntotal=1000000 ! no. of runs for weighted events +ccAM nev=100 ! no. of unweighted events generated to event record + +ccccccccc + + if(formf.eq.'exp')then ! Set parameters for Pomeron-Meson form factor in function 'fpi(x)' + bexp=1d0/2.2d0 + elseif(formf.eq.'orexp')then + bo=1d0/1.1d0 + ao=dsqrt(0.5d0) + elseif(formf.eq.'power')then + aoo=1.7d0 + endif + +ccccccccccccccccccccccccccccccccccccccc + +ccccccccc Start of main body of code + +cccccccccccccccccccccccccccccccccccccc + + call initpars(iin) ! Initialise soft survival parameters + call calcop ! proton opacity + call calcscreen ! screening amplitude + + if(pflag.eq.'pipm')then + mmes=0.13957018d0 ! pi+/- mass, PDG 2011 value + sig0=13.63d0 + alphapm=0.7d0 + alpha0m=-0.7d0*mmes**2 + cf=31.79d0/0.389d0 + crho=4.23d0/0.389d0 + + nup=6 + istup(5)=1 + idup(5)=211 + mothup(1,5)=1 + mothup(2,5)=2 + icolup(1,5)=0 + icolup(2,5)=0 + vtimup(5)=0 + spinup(5)=9. + + istup(6)=1 + idup(6)=-211 + mothup(1,6)=1 + mothup(2,6)=2 + icolup(1,6)=0 + icolup(2,6)=0 + vtimup(6)=0 + spinup(6)=9. + + elseif(pflag.eq.'pi0')then + mmes=0.1349766d0 ! pi0 mass, PDG 2011 value + sig0=13.63d0 + alphapm=0.7d0 + alpha0m=-0.7d0*mmes**2 + cf=31.79d0/0.389d0 + crho=0d0 + + nup=6 + istup(5)=1 + idup(5)=111 + mothup(1,5)=1 + mothup(2,5)=2 + icolup(1,5)=0 + icolup(2,5)=0 + vtimup(5)=0 + spinup(5)=9. + + istup(6)=1 + idup(6)=-111 + mothup(1,6)=1 + mothup(2,6)=2 + icolup(1,6)=0 + icolup(2,6)=0 + vtimup(6)=0 + spinup(6)=9. + + elseif(pflag.eq.'kpkm')then + mmes=0.493677d0 ! K+/- mass, PDG 2011 value + sig0=11.82d0 + cf=17.255d0/0.389d0 + crho=9.105d0/0.389d0 + + nup=6 + istup(5)=1 + idup(5)=321 + mothup(1,5)=1 + mothup(2,5)=2 + icolup(1,5)=0 + icolup(2,5)=0 + vtimup(5)=0 + spinup(5)=9. + + istup(6)=1 + idup(6)=-321 + mothup(1,6)=1 + mothup(2,6)=2 + icolup(1,6)=0 + icolup(2,6)=0 + vtimup(6)=0 + spinup(6)=9. + + elseif(pflag.eq.'ks')then + mmes=0.497614d0 ! K_0 mass, PDG 2011 value + sig0=11.82d0 + cf=17.255d0/0.389d0 + crho=0d0 + + nup=6 + istup(5)=1 + idup(5)=310 + mothup(1,5)=1 + mothup(2,5)=2 + icolup(1,5)=0 + icolup(2,5)=0 + vtimup(5)=0 + spinup(5)=9. + + istup(6)=1 + idup(6)=310 + mothup(1,6)=1 + mothup(2,6)=2 + icolup(1,6)=0 + icolup(2,6)=0 + vtimup(6)=0 + spinup(6)=9. + + elseif(pflag.eq.'rho')then + mmes0=0.77549d0 ! rho mass, PDG 2013 value + mwidth=0.1491d0 ! rho width PDG 2013 value + sig0=10d0 + cf=0d0 + crho=0d0 + + nup=10 + istup(5)=2 + idup(5)=113 + mothup(1,5)=1 + mothup(2,5)=2 + icolup(1,5)=0 + icolup(2,5)=0 + vtimup(5)=0 + spinup(5)=9. + + istup(6)=2 + idup(6)=113 + mothup(1,6)=1 + mothup(2,6)=2 + icolup(1,6)=0 + icolup(2,6)=0 + vtimup(6)=0 + spinup(6)=9 + + do k=7,10 + istup(k)=1 + mothup(2,k)=0 + icolup(1,k)=0 + icolup(2,k)=0 + vtimup(k)=0 + spinup(k)=9. + enddo + + idup(7)=211 + idup(8)=-211 + idup(9)=211 + idup(10)=-211 + mothup(1,7)=5 + mothup(1,8)=5 + mothup(1,9)=6 + mothup(1,10)=6 + + endif + +c -- other parameters + + ebeam=rts/2d0 + s=rts*rts + zi=(0d0,1d0) + rt2=dsqrt(2d0) + pi=dacos(-1d0) + bp=rts/dsqrt(2.0d0) + mp=0.93827d0 + beta=dsqrt(1d0-4d0*mp**2/s) + s0=1d0 + +c Pomeron + t-slope + + bb=4d0 + bjac=6d0 + bjac1=2d0 + + alphap=0.25d0 ! D-L '92 fit + alpha0=1.0808d0 + alphapr=0.93d0 + alpha0r=0.55d0 + + mf127=1.275d0 + mf1525=1.525d0 + + cpom=sig0/0.389d0 + aff=-0.860895d0 + ar=-1.16158d0 + + mmes1=mmes + mmes2=mmes + +ccccc Initialise rambo (rho0 decay) + + if(pflag.eq.'rho')then + + mmes=mmes0 + npart=2 + do j=1,4 + am(j)=0.13957018d0 + enddo + endif + +ccccccc + +c initialise counters etc + +c do hh=1,20 + + nhist=1 + sum=0d0 + sum1=0d0 + ncut=0 + + weightm=0d0 + +c initialise histograms + + do i=1,nhist + call histo3(i) + enddo + + num=0 + +cccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c Incoming protons to event record array + +cccccccccccccccccccccccccccccccccccccccccccccccccccccc + + ID(1)=2212 + q(1,1)=0d0 + q(2,1)=0d0 + q(3,1)=ebeam*beta + q(4,1)=ebeam + + istup(1)=-1 + idup(1)=2212 + mothup(1,1)=0 + mothup(2,1)=0 + icolup(1,1)=0 + icolup(2,1)=0 + do j=1,4 + pup(j,1)=q(j,1) + enddo + pup(5,1)=dsqrt(q(4,1)**2-q(3,1)**2-q(2,1)**2-q(1,1)**2) + vtimup(1)=0 + spinup(1)=9. + + q(1,2)=0d0 + q(2,2)=0d0 + q(3,2)=-ebeam*beta + q(4,2)=ebeam + + istup(2)=-1 + if(ppbar.eq.'true')then + idup(2)=-2212 + else + idup(2)=2212 + endif + mothup(1,2)=0 + mothup(2,2)=0 + icolup(1,2)=0 + icolup(2,2)=0 + do j=1,4 + pup(j,2)=q(j,2) + enddo + pup(5,2)=dsqrt(q(4,2)**2-q(3,2)**2-q(2,2)**2-q(1,2)**2) + vtimup(2)=0 + spinup(2)=9. + +cccc outgoing initial info + + istup(3)=1 + if(ppbar.eq.'true')then + idup(2)=-2212 + else + idup(2)=2212 + endif + idup(3)=2212 + mothup(1,3)=1 + mothup(2,3)=2 + icolup(1,3)=0 + icolup(2,3)=0 + vtimup(3)=0 + spinup(3)=9 + + istup(4)=1 + idup(4)=2212 + mothup(1,4)=1 + mothup(2,4)=2 + icolup(1,4)=0 + icolup(2,4)=0 + vtimup(4)=0 + spinup(4)=9 + +ccc HEPEVT + if(output.eq.'hepevt')then + + nhep=nup + write(6,*) "hepevt", nhep, nup + do k=1,5 + phep(k,1)=pup(k,1) + phep(k,2)=pup(k,2) + enddo + + do k=1,nup + isthep(k)=istup(k) + idhep(k)=idup(k) + jmohep(1,k)=mothup(1,k) + jmohep(2,k)=jmohep(1,k) + jdahep(1,k)=0 + jdahep(2,k)=0 + vhep(1,k)=0d0 + vhep(2,k)=0d0 + vhep(3,k)=0d0 + vhep(4,k)=0d0 + enddo + + if(pflag.eq.'rho')then + jdahep(1,5)=7 + jdahep(2,5)=8 + jdahep(1,6)=9 + jdahep(2,6)=10 + endif + + jmohep(2,5)=2 + jmohep(2,6)=2 + + endif + + if(unw.eq.'true')then + lmax=2 + else + lmax=1 + endif + +c do ll=1,lmax + + if(ll.eq.2)then +c ntotal=nev*10 + ntotal=1000000000 + endif + + ip=ntotal+1 + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c c +c Start of event loop c +c c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + if(ll.eq.1)then + write(6,*)'Generating weighted events...' + else + write(6,*)'Generating unweighted events...' + endif + return + end + + subroutine dimegenerate + implicit none + +ccccc hepevt output + integer nmxhep,kk + parameter (nmxhep=4000) + integer nevhep,nhep,isthep,idhep,jmohep,jdahep + double precision phep,vhep + common /hepevt/ nevhep,nhep,isthep(nmxhep),idhep(nmxhep), + &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep) +ccccc Les Houches Event Common Block + INTEGER MAXNUP + PARAMETER (MAXNUP=500) + INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP + DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP + COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP, + & IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP), + & ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP), + & SPINUP(MAXNUP) +ccccc rambo variables + integer npart,nrun + double precision pin(4), am(100), pout(4,100), wt, ein + common/ramboc/ pin, am, pout, wt, ein, + &npart, nrun + +ccccc + double precision etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax, + &rmin,mcut + common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut +ccccc + double precision q(4,20) + common /mom/ q +ccccc + double precision pt1(2), pt2(2), ptx(2) + common/mompt/pt1,pt2,ptx +ccccc + double precision s, rts, mmes, yx + double precision bjac, bp + double precision mf127, mf1525 + + double precision m0, mmes0, mmes1, mmes2, mp, mwidth, pi, rt2, + &ebeam, sh, th, uh, sum, sum1, weightm + + integer id(20), ncut, ll, icut, nev, num, iin + + common/vars/s,rts,mmes,yx + common/hvars/sh,th,uh + common/vars0/ mf127, mf1525, m0, mmes0, mmes1, mmes2, mp, + & mwidth, pi, rt2, ebeam, sum, sum1, weightm, + & bjac, bp + common /ivars/ ncut, ll, icut, nev, num +ccccc + character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10 + &,ppbar*10,output*10,mregge*10,cuts*10,unw*10 + common/flags/ iin, pflag, fsi, ppbar, output, cuts, unw + common/ff/formf + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +cccccc local variables ... + double precision aa1, aa2, al1, al2, almax, almin, c, cc1, cc2, + & msmax, msmin, cont, exn, mmesp, mwidth1, mwidth2, norm, + & p1m, p1p, p2m, p2p, phi1, phi2, phix1, pt1sq, pt1x, pt1y, + & pt2sq, pt2x, pt2y, ptxsq1, ptxsq2, ptxsqmax, ptxsqmin, + & ran0, ran1, ran2, ran3, ran4, ran5, ran6, ranhist, ranx1, + & ranx2, ranx3, ranx4, rm1, rm2, rmx1, rmx2, + & root1sq, root2sq, snp, t1, t2, weight, weight0, weight1, + & weight2, wthist, x1, x2, ymax1, ymax2, ymin1, ymin2, yx1, + & yx2 + double precision pboo(4), pcm(4), plb(4) + double precision svec(4) + complex*16 zmat + integer h, i, j, k, ntotal + weight=0d0 + call r2455(ran0) + call r2455(ran1) + call r2455(ran2) + call r2455(ran3) + call r2455(ran4) + call r2455(ran5) + call r2455(ran6) + +ccccccccccccccccccccccccccccccccccccccccccccccccccc + +c incoming protons + + ID(1)=2212 + q(1,1)=0d0 + q(2,1)=0d0 + q(3,1)=ebeam + q(4,1)=ebeam + + ID(2)=2212 + q(1,2)=0d0 + q(2,2)=0d0 + q(3,2)=-ebeam + q(4,2)=ebeam + + phi1=2d0*pi*ran0 + phi2=2d0*pi*ran1 + + pt1sq=-dlog(ran2)/(bjac) + pt2sq=-dlog(ran3)/(bjac) + + weight=(dexp(bjac*pt1sq)*dexp(bjac*pt2sq))/bjac**2 + + pt1(1)=dsqrt(pt1sq)*dsin(phi1) + pt1(2)=dsqrt(pt1sq)*dcos(phi1) + pt2(1)=dsqrt(pt2sq)*dsin(phi2) + pt2(2)=dsqrt(pt2sq)*dcos(phi2) + + pt1x=pt1(1) + pt1y=pt1(2) + pt2x=pt2(1) + pt2y=pt2(2) + +cccccccccccccccccccccccccccccccc + +ccc non-zero rho width + +cccccccccccccccccccccccccccccccc + + if(pflag.eq.'rho')then ! new part here + + 677 call r2455(rm1) + call r2455(rm2) + + msmax=mmes0+4d0*mwidth + msmin=2d0*0.13957018d0 + +cccccccccccccccc + + almin=datan(-(mmes0**2-msmin**2)/mwidth/mmes0) + almax=datan(-(mmes0**2-msmax**2)/mwidth/mmes0) + + al1=almin+(almax-almin)*rm1 + al2=almin+(almax-almin)*rm2 + + mmes1=dsqrt(dtan(al1)*mmes0*mwidth+mmes0**2) + mmes2=dsqrt(dtan(al2)*mmes0*mwidth+mmes0**2) + + weight=weight*(almax-almin)**2 + weight=weight*mwidth**2*mmes0**2 + weight=weight*(1d0+dtan(al1)**2)*(1d0+dtan(al2)**2) + weight=weight/4d0/mmes1/mmes2 + +ccccccccccccccccc + + mwidth1=mwidth*((1d0-4d0*0.13957018d0**2/mmes1**2)/ + & (1d0-4d0*0.13957018d0**2/mmes0**2))**1.5d0 + mwidth2=mwidth*((1d0-4d0*0.13957018d0**2/mmes2**2)/ + & (1d0-4d0*0.13957018d0**2/mmes0**2))**1.5d0 + + if(mmes1.lt.2d0*0.13957018d0) goto 677 + if(mmes2.lt.2d0*0.13957018d0) goto 677 + + weight=weight*2d0*mmes0*mmes1*mwidth1/pi + weight=weight*2d0*mmes0*mmes2*mwidth2/pi + weight=weight*mmes1**2*mmes2**2/mmes0**4 + weight=weight/((mmes0**2-mmes1**2)**2+mmes1**4* + & mwidth1**2/mmes0**2) + weight=weight/((mmes0**2-mmes2**2)**2+mmes2**4* + & mwidth2**2/mmes0**2) + + endif + + +cccccccccccccccccccccccccccccccccccccccccccccccccccc + +c pi rapidities + + call r2455(ranx1) + call r2455(ranx2) + call r2455(ranx3) + call r2455(ranx4) + + phix1=2d0*pi*ranx1 + +c r2max=1d0 +c r2=r2max*ranx1 +c ptxsq1=-dlog(r2)/bjac1 ! generate exp p_t^2 (can be more efficient) + + ptxsqmin=0d0 + ptxsqmax=10d0 ! generate linear p_t^2 + + ptxsq1=ptxsqmin+(ptxsqmax-ptxsqmin)*ranx2 + + q(1,5)=dsqrt(ptxsq1)*dcos(phix1) + q(2,5)=dsqrt(ptxsq1)*dsin(phix1) + q(1,6)=-pt1(1)-pt2(1)-q(1,5) + q(2,6)=-pt1(2)-pt2(2)-q(2,5) + ptxsq2=q(1,6)**2+q(2,6)**2 + + rmx1=dsqrt(mmes1**2+ptxsq1) + rmx2=dsqrt(mmes2**2+ptxsq2) + + ymax1=dlog(rts/rmx1) + + if(cuts.eq.'true')then + if(rmax.lt.ymax1)then + ymax1=rmax + endif + endif + + ymin1=-ymax1 + + yx1=ymin1+(ymax1-ymin1)*ranx3 + + ymax2=(dlog((rts-rmx1*dexp(yx1))/rmx2)) + ymin2=-(dlog((rts-rmx1*dexp(-yx1))/rmx2)) + + if(cuts.eq.'true')then + if(ymax2.gt.rmax)then + ymax2=rmax + endif + if(ymin2.lt.-rmax)then + ymin2=-rmax + endif + endif + + if(ymax2.lt.ymin2)then + weight=0d0 + goto 700 + endif + + yx2=ymin2+(ymax2-ymin2)*ranx4 + x1=(rmx2*dexp(yx2)+rmx1*dexp(yx1))/rts + x2=(rmx2*dexp(-yx2)+rmx1*dexp(-yx1))/rts + + weight=weight*(ptxsqmax-ptxsqmin) +c weight=weight*dexp(bjac1*ptxsq1)*r2max/bjac1 + weight=weight*(ymax1-ymin1)*(ymax2-ymin2) + + + q(3,5)=rmx1*(dexp(yx1)-dexp(-yx1))/2d0 + q(4,5)=rmx1*(dexp(yx1)+dexp(-yx1))/2d0 + + q(3,6)=rmx2*(dexp(yx2)-dexp(-yx2))/2d0 + q(4,6)=rmx2*(dexp(yx2)+dexp(-yx2))/2d0 + +ccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c impose massive on-shell condition by solving +c p1+ + cc1/p2- = aa1 +c p2- + cc2/p1+ = aa2 + + aa1=bp*(1d0-x1) + aa2=bp*(1d0-x2) + cc1=0.5d0*(pt2sq+mp**2) + cc2=0.5d0*(pt1sq+mp**2) + + root1sq=(cc1-cc2-aa1*aa2)**2-4d0*cc2*aa1*aa2 + root2sq=(cc2-cc1-aa1*aa2)**2-4d0*cc1*aa1*aa2 + if(root1sq.le.0d0.or.root2sq.le.0d0)then + weight=0d0 + goto 700 + endif + p1p=(cc2-cc1+aa1*aa2+dsqrt(root1sq))/(2d0*aa2) + p2m=(cc1-cc2+aa1*aa2+dsqrt(root2sq))/(2d0*aa1) + p1m=(pt1sq+mp**2)/(2d0*p1p) + p2p=(pt2sq+mp**2)/(2d0*p2m) + + if(p1p.lt.0d0.or.p1m.lt.0d0.or.p2p.lt.0d0.or.p2m.lt.0d0)then + weight=0d0 + goto 700 + endif + + t1=-rts*p1m*rt2 + t2=-rts*p2p*rt2 + + q(1,3)=pt1(1) + q(2,3)=pt1(2) + q(3,3)=(p1p-p1m)/rt2 + q(4,3)=(p1p+p1m)/rt2 + + q(1,4)=pt2(1) + q(2,4)=pt2(2) + q(3,4)=(p2p-p2m)/rt2 + q(4,4)=(p2p+p2m)/rt2 + +ccccccccccccccccccccccccccccccccccccccccccccccccc + + do k=1,4 + svec(k)=q(k,5)+q(k,6) + enddo + + call dot(svec,svec,sh) + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c rho0 decays + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + if(pflag.eq.'rho')then + + kk=6 + + do k=5,6 + + if(k.eq.5)then + mmesp=mmes1 + else + mmesp=mmes2 + endif + call rambo(npart,mmesp,am,pout,wt) ! call RAMBO + + do j=1,4 + pboo(j)=q(j,k) + enddo + + do h=1,2 + do j=1,4 + pcm(j)=pout(j,h) + enddo + call boost(mmesp,pboo,pcm,plb) + kk=kk+1 + do j=1,4 + q(j,kk)=plb(j) + enddo + enddo + + enddo + + endif + + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c place cuts + + if(cuts.eq.'true')then + + call cut(icut) + + if(icut.eq.0)then + weight=0d0 + weight0=0d0 + weight1=0d0 + weight2=0d0 + ncut=ncut+1 + goto 700 + endif + + endif + +c print*,i + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c Write 4-momenta to event record array + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + if(output.eq.'lhe')then + + do k=3,4 + do j=1,4 + pup(j,k)=q(j,k) + enddo + pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2) + enddo + + if(pflag.eq.'rho')then + + do k=5,10 + do j=1,4 + pup(j,k)=q(j,k) + enddo + pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2) + enddo + + else + + do k=5,6 + do j=1,4 + pup(j,k)=q(j,k) + enddo + pup(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2) + enddo + + endif + + else + + do k=3,4 + do j=1,4 + phep(j,k)=q(j,k) + enddo + phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2) + enddo + + if(pflag.eq.'rho')then + + do k=5,10 + do j=1,4 + phep(j,k)=q(j,k) + enddo + phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2) + enddo + + else + + do k=5,6 + do j=1,4 + phep(j,k)=q(j,k) + enddo + phep(5,k)=dsqrt(q(4,k)**2-q(3,k)**2-q(2,k)**2-q(1,k)**2) + enddo + + endif + + endif + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c Weight evaluation + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c pp-->pp(M_1 M_2) matrix element + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + + call schimc(pt1x,pt1y,pt2x,pt2y,zmat) + + weight=weight*cdabs(zmat)**2 + weight=weight/norm**2 + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + weight=weight/(s**2*(16d0)**3*pi**5)*389d3 + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +ccc Probability of no additional particles produced in production subprocess + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + if(fsi.eq.'true')then + + if(pflag.eq.'pipm'.or.pflag.eq.'pi0')then + + if(dsqrt(sh).gt.mf127)then + m0=mf127 + c=0.7d0 + cont=1d0-(m0-2d0*mmes)**2/m0**2 + exn=c*dlog(((dsqrt(sh)-2d0*mmes)**2)/m0**2+cont) + snp=dexp(-exn) + else + snp=1d0 + endif + + elseif(pflag.eq.'kpkm'.or.pflag.eq.'ks')then + + if(dsqrt(sh).gt.mf1525)then + m0=mf1525 + c=0.7d0 + cont=1d0-(m0-2d0*mmes)**2/m0**2 + exn=c*dlog(((dsqrt(sh)-2d0*mmes)**2)/m0**2+cont) + snp=dexp(-exn) + else + snp=1d0 + endif + + else + + snp=1d0 + + endif + + else + + snp=1d0 + + endif + + if(snp.gt.1d0)then + snp=1d0 + endif + + weight=weight*snp + +ccccccccccccccc + + if(pflag.eq.'pi0'.or.pflag.eq.'ks'.or.pflag.eq.'rho')then + weight=weight/2d0 ! symmetry factor + endif + +cccccccccccccccccccccccccccccccccccccccccccccccccccccc + + 666 weight=weight + + if(weightm.lt.weight)then + weightm=weight + endif + + + if(ll.eq.1)then + + xwgtup=weight + + wthist=weight/dfloat(ntotal) + + call binit(wthist) + + else + + call r2455(ranhist) + + if(ranhist.lt.weight/weightm)then + + xwgtup=1d0 + + num=num+1 + + +c call binit(sumt/nev) + + write(35,302)num,nup + + if(output.eq.'lhe')then + + do j=1,nup + write(35,300)j,istup(j),idup(j),mothup(1,j),mothup(2,j), + & icolup(1,j),icolup(2,j),pup(1,j),pup(2,j), + & pup(3,j),pup(4,j),pup(5,j),vtimup(j),spinup(j) + enddo + + else + + do j=1,nup + write(35,301)j,idhep(j),isthep(j),jmohep(1,j), + & jmohep(2,j),jdahep(1,j),jdahep(2,j), + & phep(1,j),phep(2,j),phep(3,j),phep(4,j),phep(5,j) + & ,vhep(1,j),vhep(2,j),vhep(3,j),vhep(4,j) + enddo + + endif + + + endif + + endif + + if(ll.eq.2)then + + if(num.gt.nev-1)then ! exit loop once required no. of unweighted events generated + + ntotal=i +! goto 888 + call terminate() + + endif + endif + + 700 sum=sum+weight + sum1=sum1+weight**2 + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c c +c End of event loop c +c c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + call terminate() + 300 format(i4,1x,i4,1x,i8,1x,i4,1x,i4,1x,i4,1x,i4,1x,E13.6,1x, + &E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6) + 301 format(i4,1x,i8,1x,i4,1x,i4,1x,i4,1x,i4,1x,i4,1x,E13.6,1x,E13.6 + &,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x,E13.6,1x + &,E13.6) + 302 format(1x,i8,1x,i5,1x,E13.6) + + return + end + + subroutine terminate + character prefix*50,fsp*10,order*10,pflag*10,fsi*10,formf*10 + &,ppbar*10,output*10,mregge*10,cuts*10,unw*10 + + 888 continue + sum=sum/dfloat(ntotal) + sum1=sum1/dfloat(ntotal) + var=dsqrt((sum1-sum**2)/dfloat(ntotal)) + EFF=1d0*(ntotal-ncut)/ntotal + + if(ll.eq.1)then + print*,'...done!' + write(6,*) + sumt=sum + else + print*,'...done!' + write(6,*) + endif + + + +ccccc create histograms + + if(ll.eq.1)then + + do i=1,nhist + call histo2(i,0) + enddo + + endif + + if(ll.eq.1)then + if(pflag.eq.'pipm')then + write(6,*) + write(6,*)'pi+pi- production' + elseif(pflag.eq.'pi0')then + write(6,*) + write(6,*)'pi0pi0 production' + elseif(pflag.eq.'kpkm')then + write(6,*) + write(6,*)'K+K- production' + elseif(pflag.eq.'ks')then + write(6,*) + write(6,*)'K0K0 production' + elseif(pflag.eq.'rho')then + write(6,*) + write(6,*)'rho0rho0 production' + endif + endif + + if(ll.eq.1)then + if(pflag.eq.'rho')then + write(6,221)rts,ntotal,sum,var + else + write(6,222)rts,mmes,ntotal,sum,var + endif + else + if(output.eq.'lhe')then + print*,'LHE output' + else + print*,'HEPEVT output' + endif + write(6,223)nev,sum,var + endif + +! enddo + + close(35) + + 221 format(/, + . 3x,'collider energy (GeV) : ',f10.3,/, + . 3x,'number of events : ',i9,/, + . 3x,'sigma (nb) : ',f15.5,' +- ',f15.5, + . /) + + 222 format(/, + . 3x,'collider energy (GeV) : ',f10.3,/, + . 3x,'meson mass (GeV) : ',f10.5,/, + . 3x,'number of events : ',i9,/, + . 3x,'sigma (nb) : ',f15.5,' +- ',f15.5, + . /) + + 223 format(3x,'number of events : ',i9,/, + . 3x,'sigma (nb) : ',f15.5,' +- ',f15.5, + . /) + + call cpu_time(t2) + print*,'time elapsed = ', t2, ' s' + + return + end + +ccccc Pomeron -- (off-shell) meson form factor + + function fpi(x) + implicit double precision(a-y) + character formf*10 + common/vars/s,rts,mmes,yx + common/ff/formf + common/ffpars/bexp,ao,bo,aoo,boo + + if(formf.eq.'exp')then + fpi=exp((x-mmes**2)*bexp) + elseif(formf.eq.'orexp')then + fpi=exp(-(dsqrt(-x+mmes**2+ao**2)-ao)*bo) + elseif(formf.eq.'power')then + fpi=1d0/(1d0-(x-mmes**2)/aoo) + endif + + return + end + +ccccccc Pom Pom --> meson pair amplitude + + subroutine wev(matt,v1,v2) + implicit double precision(a-y) + implicit complex*16(z) + double precision q(4,20),svec(4),tvec(4),uvec(4),v1(4),v2(4) + integer k + character*10 mregge,pflag + complex*16 matt + common/mom/q + common/vars/s,rts,mmes,yx + common/wvars/sig0,bb,t1,t2 + common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m + common/sec/cpom,cf,crho,aff,ar + common/hvars/sh,th,uh + common/regge/mregge + common/flags/ iin, pflag, fsi, ppbar, output, cuts, unw + + + zi=(0d0,1d0) + pi=dacos(-1d0) + + do k=1,4 + q(k,10)=v1(k) + q(k,11)=v2(k) + enddo + + do k=1,4 + tvec(k)=q(k,10)-q(k,5) + uvec(k)=q(k,10)-q(k,6) + enddo + + call dot(tvec,tvec,th) + call dot(uvec,uvec,uh) + + sht1=2d0*(q(4,5)*q(4,3)-q(3,5)*q(3,3)-q(2,5)*q(2,3)- + &q(1,5)*q(1,3))+mmes**2 + + sht2=2d0*(q(4,6)*q(4,4)-q(3,6)*q(3,4)-q(2,6)*q(2,4)- + &q(1,6)*q(1,4))+mmes**2 + + shu1=2d0*(q(4,6)*q(4,3)-q(3,6)*q(3,3)-q(2,6)*q(2,3)- + &q(1,6)*q(1,3))+mmes**2 + + shu2=2d0*(q(4,5)*q(4,4)-q(3,5)*q(3,4)-q(2,5)*q(2,4)- + &q(1,5)*q(1,4))+mmes**2 + + bt1=(2d0*alphap*dlog(sht1)) + bt2=(2d0*alphap*dlog(sht2)) + bu1=(2d0*alphap*dlog(shu1)) + bu2=(2d0*alphap*dlog(shu2)) + + bt1r=(2d0*alphapr*dlog(sht1)) + bt2r=(2d0*alphapr*dlog(sht2)) + bu1r=(2d0*alphapr*dlog(shu1)) + bu2r=(2d0*alphapr*dlog(shu2)) + + sig0t1=sig0*sht1**0.0808d0/0.389d0 + sig0t2=sig0*sht2**0.0808d0/0.389d0 + sig0u1=sig0*shu1**0.0808d0/0.389d0 + sig0u2=sig0*shu2**0.0808d0/0.389d0 + + t11=v1(1)**2+v1(2)**2 + t22=v2(1)**2+v2(2)**2 + + zmt=(zi*dexp(-bt1*t11/2d0)*cpom*sht1**alpha0 + & +((aff+zi)*cf*sht1**alpha0r+(ar-zi)*crho*sht1**alpha0r) + & *dexp(-bt1r*t11/2d0)) + zmt=zmt*(zi*dexp(-bt2*t22/2d0)*cpom*sht2**alpha0 + & +((aff+zi)*cf*sht2**alpha0r-(ar-zi)*crho*sht2**alpha0r) + & *dexp(-bt2r*t11/2d0)) + zmt=zmt*fpi(th)**2/(mmes**2-th) + + zmu=(zi*dexp(-bu1*t11/2d0)*cpom*shu1**alpha0 + & +((aff+zi)*cf*shu1**alpha0r-(ar-zi)*crho*shu1**alpha0r) + & *dexp(-bu1r*t11/2d0)) + zmu=zmu*(zi*dexp(-bu2*t22/2d0)*cpom*shu2**alpha0 + & +((aff+zi)*cf*shu2**alpha0r+(ar-zi)*crho*shu2**alpha0r) + & *dexp(-bu2r*t22/2d0)) + zmu=zmu*fpi(uh)**2/(mmes**2-uh) + + matt=zmu+zmt + + return + end + +c binning subroutine + + subroutine binit(wt) + implicit double precision(a-y) + double precision q(4,20),pt1(2),pt2(2),ptx(2) + common/mom/q + common/mompt/pt1,pt2,ptx + common/vars/s,rts,mmes,yx + common/vars1/ptgam,etagam,ptel2,ptel1,etael1,etael2 + common/vars2/ptpi1,etapi1,ptpi2,etapi2 + common/ang/cost,costa,costb + common/hvars/sh,th,uh + + ypisum=0.5d0*dlog((q(4,5)+q(4,6)+q(3,5)+q(3,6)) + & /(q(4,5)+q(4,6)-q(3,5)-q(3,6))) + + ypi1=0.5d0*dlog((q(4,5)+q(3,5)) + & /(q(4,5)-q(3,5))) + ypi2=0.5d0*dlog((q(4,6)+q(3,6)) + & /(q(4,6)-q(3,6))) + + call histo1(1,30,0d0,4.5d0,dsqrt(sh),wt) + + return + end + +ccccc Cutting subroutine +ccccc +ccccc User can define their own cuts on any particle momenta, stored in array q(i,j) +ccccc + + subroutine cut(icut) + implicit double precision(a-y) + double precision q(4,20) + integer icut + common/mom/q + common/vars/s,rts,mmes,yx + common/vars1/ptgam,etagam,ptel2,ptel1,etael1,etael2 + common/vars2/ptpi1,etapi1,ptpi2,etapi2 + common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ecut,rmax,rmin,mcut + common/hvars/sh,th,uh + integer iii + common/ivar/iii + + icut=0 + +c -- meson 1 variables + ptpi1=dsqrt(q(1,5)**2+q(2,5)**2) + pmodpi1=dsqrt(q(1,5)**2+q(2,5)**2+q(3,5)**2) + etapi1=.5d0*dlog((pmodpi1+q(3,5)) + . /(pmodpi1-q(3,5))) + + ypi1=.5d0*dlog((q(4,5)+q(3,5)) + . /(q(4,5)-q(3,5))) + + etpi1=q(4,5)*dsqrt(pmodpi1**2-q(3,5)**2)/pmodpi1 + +c -- meson 2 variables + ptpi2=dsqrt(q(1,6)**2+q(2,6)**2) + pmodpi2=dsqrt(q(1,6)**2+q(2,6)**2+q(3,6)**2) + etapi2=.5d0*dlog((pmodpi2+q(3,6)) + . /(pmodpi2-q(3,6))) + + ypi2=.5d0*dlog((q(4,6)+q(3,6)) + . /(q(4,6)-q(3,6))) + + etpi2=q(4,6)*dsqrt(pmodpi2**2-q(3,6)**2)/pmodpi2 + +c print*,ptpi1,ptpi2 + +ccccc example cuts.... + +c if(dsqrt(sh).lt.0.8d0)return + + xf1=2d0*q(3,3)/rts + xf2=2d0*q(3,4)/rts + +c if(dabs(xf1).lt.0.9d0)return +c if(dabs(xf2).lt.0.9d0)return + +ccccccccc Basic cuts described at start of code + + if(ptpi1.lt.ecut)return + if(ptpi2.lt.ecut)return + if(ypi1.gt.rmax)return + if(ypi2.gt.rmax)return + if(ypi1.lt.rmin)return + if(ypi2.lt.rmin)return + + icut=1 + + return + end + +c prints histograms + + subroutine histo1(ih,ib,x0,x1,x,w) + implicit real*8(a-h,o-z) + character*1 regel(30),blank,star + dimension h(20,100),hx(20),io(20),iu(20),ii(20) + dimension y0(20),y1(20),ic(20) + data regel / 30*' ' /,blank /' ' /,star /'*'/ + save + open(10,file="output.dat") + y0(ih)=x0 + y1(ih)=x1 + ic(ih)=ib + if(x.lt.x0) goto 11 + if(x.gt.x1) goto 12 + ix=idint((x-x0)/(x1-x0)*dble(ib))+1 + h(ih,ix)=h(ih,ix)+w + if(h(ih,ix).gt.hx(ih)) hx(ih)=h(ih,ix) + ii(ih)=ii(ih)+1 + return + 11 iu(ih)=iu(ih)+1 + return + 12 io(ih)=io(ih)+1 + return + entry histo2(ih,il) + ib1=ic(ih) + x01=y0(ih) + x11=y1(ih) + bsize=(x11-x01)/dble(ib1) + hx(ih)=hx(ih)*(1.d0+1.d-06) + if(il.eq.0) write(6,21) ih,ii(ih),iu(ih),io(ih) + if(il.eq.1) write(6,22) ih,ii(ih),iu(ih),io(ih) + 21 format(' no.',i3,' lin : inside,under,over ',3i6) + 22 format(' no.',i3,' log : inside,under,over ',3i6) + if(ii(ih).eq.0) goto 28 + write(6,23) + 23 format(35(1h ),3(10h----+----i)) + do 27 iv=1,ib1 + z=(dble(iv)-0.5d0)/dble(ib1)*(x11-x01)+x01 + if(il.eq.1) goto 24 + iz=idint(h(ih,iv)/hx(ih)*30.)+1 + goto 25 + 24 iz=-1 + if(h(ih,iv).gt.0.d0) + .iz=idint(dlog(h(ih,iv))/dlog(hx(ih))*30.)+1 + 25 if(iz.gt.0.and.iz.le.30) regel(iz)=star + write(6,26) z,h(ih,iv)/bsize,(regel(i),i=1,30) +c write(10,*)z-0.125d0/2d0,h(ih,iv)/bsize ! Print histogram to file + write(10,*)z,h(ih,iv)/bsize ! Print histogram to file +c write(10,*)z,h(ih,iv)/bsize ! Print histogram to file + 26 format(1h ,2g15.6,4h i,30a1,1hi) + 36 format(1h ,2g15.6) + if(iz.gt.0.and.iz.le.30) regel(iz)=blank + 27 continue + write(6,23) + return + 28 write(6,29) + 29 format(' no entries inside histogram') + return + entry histo3(ih) + do 31 i=1,100 + 31 h(ih,i)=0. + hx(ih)=0. + io(ih)=0 + iu(ih)=0 + ii(ih)=0 + close(10) + return + end + + +cccc Initializes soft model parameters + + subroutine initpars(in) + implicit double precision(a-y) + integer in,i1,i2 + integer nch + double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5) + common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm + common/ipars/nch + common/vars/s,rts,mmes,yx + + pi=dacos(-1d0) + + if(in.eq.1)then + + ep=0.13d0 + asp=0.08d0 + ep1=0d0 + sigo=23d0/0.39d0 + gaa(3)=0.4d0 + gd2=0.3d0 + nch=2 + ntf=0 + + cc0(1)=0.45d0 + bm(1)=3d0 + bb0(1)=0.1d0 + pp0(1)=0.92d0 + bex(1)=8.5d0 + + cc0(2)=0.45d0 + bm(2)=1.5d0 + bb0(2)=0.5d0 + pp0(2)=0.1d0 + bex(2)=4.5d0 + + cc0(3)=1d0 + bm(3)=0d0 + bb0(3)=0.8d0 + pp0(3)=0.5d0 + bex(3)=0.5d0 + + elseif(in.eq.2)then + + ep=0.115d0 + asp=0.11d0 + ep1=0d0 + sigo=33d0/0.39d0 + gaa(3)=0.6d0 + gd2=0.16d0 + nch=2d0 + ntf=0d0 + + cc0(1)=0.63d0 + bm(1)=3d0 + bb0(1)=0.1d0 + pp0(1)=0.5d0 + bex(1)=8d0 + + cc0(2)=0.47d0 + bm(2)=1.5d0 + bb0(2)=0.5d0 + pp0(2)=0.1d0 + bex(2)=6d0 + + cc0(3)=1d0 + bm(3)=0d0 + bb0(3)=0.8d0 + pp0(3)=0.5d0 + bex(3)=0.5d0 + + elseif(in.eq.3)then + + ep=0.093d0 + asp=0.075d0 + ep1=0d0 + sigo=60d0/0.39d0 + gaa3=4.8d0 + gd2=1.03d0 + nch=2 + nga=1 + cc0(1)=0.55d0 + bm(1)=3d0 + bb0(1)=0.27d0 + pp0(1)=0.48d0 + bex(1)=5.3d0 + + cc0(2)=0.48d0 + bm(2)=1.5d0 + bb0(2)=0.1d0 + pp0(2)=1d0 + bex(2)=3.8d0 + + cc0(3)=0.24d0 + + elseif(in.eq.4)then + + ep=0.11d0 !!! Capital delta + asp=0.06d0 !!! alpha' + ep1=0d0 !!! Zero in all models, matters (?) + sigo=50d0/0.39d0 !!! sigma_0 (GeV^-2) + gaa3=6d0 !!! k2/k(1.8 TeV) + gd2=1.3d0 !!! k1/k(1.8 TeV) + nch=2 + nga=1 !!! A flag (?) + cc0(1)=0.6d0 !!! d1 + bm(1)=3d0 !!! Doesn't matter + bb0(1)=0.45d0 !!! c1-0.08 (added back later) + pp0(1)=0.5d0 !!! 2*|a_1|^2 + bex(1)=7.2d0 !!! b1 + + cc0(2)=0.48d0 !!! d2 + bm(2)=1.5d0 !!! Doesn't matter + bb0(2)=0.16d0 !!! c2-0.08 (added back later) + pp0(2)=1d0 !!! |a_2|^2 is set later + bex(2)=4.2d0 !!! b2 + + cc0(3)=0.12d0 !!! Beta : k^2_min ~ s^Beta + + endif + + if(nch.eq.3) pp0(3)=3d0-pp0(2)-pp0(1) + if(nch.eq.2) pp0(2)=2d0-pp0(1) !!! Set |a_2|^2 + if(nch.eq.1) pp0(1)=1d0 + + if(in.eq.3.or.in.eq.4)then + + gamm=(1800d0/rts)**cc0(3) + ga1=1d0/(1d0+gamm*gd2) + ga2=1d0/(1d0+gamm*gaa3) + gaa(1)=2d0*ga1/(ga1+ga2) !!! gamma_1 (+) Eq (15) + gaa(2)=2d0*ga2/(ga1+ga2) !!! gamma_2 (-) Eq (15) + + elseif(in.eq.1.or.in.eq.2)then + + if(nch.eq.2)then + gaa(1)=1d0+dsqrt(gd2) + gaa(2)=1d0-dsqrt(gd2) + gaa(3)=0. + elseif(nch.eq.1)then + gaa(1)=1d0 + gaa(2)=0d0 + gaa(3)=0d0 + endif + + endif + + sum=0d0 + + do i1=1,nch + do i2=1,nch + sum=sum+gaa(i1)*gaa(i2)*pp0(i1)*pp0(i2)/dble(nch)**2 + enddo + enddo + + norm=sum + + return + end + +ccccc Calculates proton opacity + + subroutine calcop + implicit double precision(a-y) + integer nb,i1,i2,nch,ib + common/ipars/nch + double precision op(5,5,10000,2),oph(5,5,10000,2) + common/opac/op,oph + + nb=900 + hb=100d0/dble(nb) + + print*,'Calculating opacity' + + + do ib=1,nb+1 + bt=dble(ib-1)*hb + + do i1=1,nch + do i2=1,nch + + call opacity(i1,i2,bt,fr,fr1) + + + op(i1,i2,ib,1)=bt + op(i1,i2,ib,2)=fr + oph(i1,i2,ib,1)=bt + oph(i1,i2,ib,2)=fr1 + + write(40,*)bt,fr,fr1 + + enddo + enddo + enddo + + return + end + +cccc Calculates screening amplitude (in k_t space) + + subroutine calcscreen + implicit double precision(a-y) + integer ns,i1,i2,nch,ib + common/ipars/nch + double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2) + common/spac/sca,sca1 + + ns=900 + ksqma=8.2d0 + inck=ksqma/dble(ns) + ksqmin=0.001d0 + lginck=(dlog(ksqma/ksqmin))/dble(ns) + + print*,'Calculating screening amplitude' + + do ib=0,ns+1 + + ksq=dble(ib-1)*inck + lgksq=ksq + + if(ib.eq.0)then + ksq=0d0 + lgksq=0d0 + else + lgksq=dble(ib-1)*lginck+dlog(ksqmin) + ksq=dexp(lgksq) + endif + + do i1=1,nch + do i2=1,nch + + call screening(i1,i2,ksq,sc,sc1) + + sca(i1,i2,ib,1)=lgksq + sca(i1,i2,ib,2)=sc + sca1(i1,i2,ib,1)=lgksq + sca1(i1,i2,ib,2)=sc1 + + enddo + enddo + + enddo + + return + end + +ccccc Interpolation.... + + subroutine opacityint(i,j,bt,out,out1) + implicit double precision(a-y) + double precision op(5,5,10000,2),oph(5,5,10000,2) + integer i,j,it + common/opac/op,oph + + incbt=op(1,1,2,1)-op(1,1,1,1) + it=nint(bt/incbt) + if(dble(it).gt.(bt/incbt))then + it=it-1 + endif + + m=(op(i,j,it+2,2)-op(i,j,it+1,2)) + &/(op(i,j,it+2,1)-op(i,j,it+1,1)) + del=bt-op(1,1,it+1,1) + mh=(oph(i,j,it+2,2)-oph(i,j,it+1,2)) + &/(oph(i,j,it+2,1)-oph(i,j,it+1,1)) + delh=bt-oph(1,1,it+1,1) + + out=m*del+op(i,j,it+1,2) + out1=mh*delh+oph(i,j,it+1,2) + + return + end + + subroutine screeningint(i,j,ktsq,out,out1) + implicit double precision(a-y) + double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2) + integer i,j,it,ns + common/spac/sca,sca1 + + if(ktsq.lt.0.001d0)then + + it=0 + + m=(sca(i,j,it+1,2)-sca(i,j,it,2)) + & /(dexp(sca(i,j,it+1,1))-sca(i,j,it,1)) + del=ktsq-sca(1,1,it,1) + m1=sca1(i,j,it+1,2)-sca1(i,j,it,2) + & /(dexp(sca1(i,j,it+1,1))-sca1(i,j,it,1)) + del1=ktsq-sca1(1,1,it,1) + + out=m*del+sca(i,j,it,2) + out1=m1*del1+sca1(i,j,it,2) + + elseif(ktsq.lt.8d0)then + + ktmin=sca(1,1,1,1) + inckt=sca(1,1,2,1)-sca(1,1,1,1) + it=nint((dlog(ktsq)-ktmin)/inckt) + if(dble(it).gt.((dlog(ktsq)-ktmin)/inckt))then + it=it-1 + endif + + m=(sca(i,j,it+2,2)-sca(i,j,it+1,2)) + & /(sca(i,j,it+2,1)-sca(i,j,it+1,1)) + del=dlog(ktsq)-sca(1,1,it+1,1) + m1=sca1(i,j,it+2,2)-sca1(i,j,it+1,2) + & /(sca1(i,j,it+2,1)-sca1(i,j,it+1,1)) + del1=dlog(ktsq)-sca1(1,1,it+1,1) + + out=m*del+sca(i,j,it+1,2) + out1=m1*del1+sca1(i,j,it+1,2) + + else + + out=0d0 + out1=0d0 + + endif + + return + end + +cccc Integrates round Pomeron loop (to calculate screened amplitude) + + subroutine schimc(p1x,p1y,p2x,p2y,out) + implicit double precision(a-y) + integer nx,jx,jy,nch,i1,i2,it,k + double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5) + double precision sca(5,5,0:40000,2),sca1(5,5,0:40000,2) + double precision a1(4),a2(4),q(4,20) + complex*16 out,x0,x00 + common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm + common/ipars/nch + common/spac/sca,sca1 + common/mom/q + common/wvars/sig0,bb,t1,t2 + common/alphas/alphap,alpha0,alphapr,alpha0r,alphapm,alpha0m + + do k=1,4 + a1(k)=q(k,1)-q(k,3) + a2(k)=q(k,2)-q(k,4) + enddo + + nx=10 + hx=2d0/dble(nx) + + out=0d0 + + do jx=-nx,nx + do jy=-nx,nx + + tpx=dble(jx)*hx + tpy=dble(jy)*hx + tp2=tpx**2+tpy**2 + wt=hx*hx + + p1xp=p1x-tpx + p1yp=p1y-tpy + t12=p1xp**2+p1yp**2 + p2xp=tpx+p2x + p2yp=tpy+p2y + t22=p2xp**2+p2yp**2 + + a1(1)=-p1xp + a1(2)=-p1yp + a2(1)=-p2xp + a2(2)=-p2yp + + do i1=1,nch + do i2=1,nch + + call screeningint(i1,i2,tp2,sc,sc1) + call wev(x0,a1,a2) + + x0=x0*pp0(i1)*pp0(i2)/dble(nch)**2 + x0=x0*dexp(-((t12+0.08d0+bb0(i1))*bex(i1))**cc0(i1)+ + & (bex(i1)*(bb0(i1)+0.08d0))**cc0(i1)) + x0=x0*dexp(-((t22+0.08d0+bb0(i2))*bex(i2))**cc0(i2)+ + & (bex(i2)*(bb0(i2)+0.08d0))**cc0(i2)) + + out=out+x0*wt*sc + + enddo + enddo + + + enddo + enddo + + a1(1)=-p1x + a1(2)=-p1y + a2(1)=-p2x + a2(2)=-p2y + + t11=p1x**2+p1y**2 + t22=p2x**2+p2y**2 + + call formfac(t11,t22,x00p) + + call wev(x00,a1,a2) + + + + out=out+x00*x00p + + return + end + +ccccc Pomeron -- diffrative eignstate i form factor + + subroutine formfac(t1,t2,out) + implicit double precision(a-y) + double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5) + common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm + integer nch,i1,i2 + common/ipars/nch + + out=0d0 + + delta1=dexp(-t1/2d0) + delta2=dexp(-t2/2d0) + + do i1=1,nch + do i2=1,nch + + wt=gaa(i1)*gaa(i2)*pp0(i1)*pp0(i2)/dble(nch)**2 + wt=wt*dexp(-((t1+0.08d0+bb0(i1))*bex(i1))**cc0(i1)+ + & (bex(i1)*(bb0(i1)+0.08d0))**cc0(i1)) + wt=wt*dexp(-((t2+0.08d0+bb0(i2))*bex(i2))**cc0(i2)+ + & (bex(i2)*(bb0(i2)+0.08d0))**cc0(i2)) + + out=out+wt + + enddo + enddo + + return + end + + subroutine screening(i,j,ktsq,out,out1) + implicit double precision(a-y) + integer i,j,ib,nb,it,nt,nch + double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5) + common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm + common/ipars/nch + common/vars/s,rts,mmes,yx + + pi=dacos(-1d0) + + nb=5000 + hb=99d0/dble(nb) + + out=0d0 + out1=0d0 + + do ib=1,nb + bt=dble(ib)*hb + wt=-bt/2d0/pi*hb + + call opacityint(i,j,bt,fr,fr1) + + sige=sigo*dexp(dlog(rts)*2d0*ep) + fr=fr*gaa(i)*gaa(j)*sige + fr1=fr1*gaa(i)*gaa(j)*sige + + out=out+wt*(1d0-dexp(-fr/2d0))*besj0(bt*dsqrt(ktsq)) + & *gaa(i)*gaa(j) + out1=out1+wt*(1d0-dexp(-fr1/2d0))*besj0(bt*dsqrt(ktsq)) + & *gaa(i)*gaa(j) + + enddo + + return + end + + + subroutine opacity(i,j,bt,out,out1) + implicit double precision(a-y) + integer i,j,ib,nb,it,nt,nch + double precision cc0(5),bm(5),bb0(5),pp0(5),bex(5),gaa(5) + common/pars/cc0,bm,bb0,pp0,bex,asp,sigo,gaa,ep,norm + common/ipars/nch + common/vars/s,rts,mmes,yx + + pi=dacos(-1d0) + + ampi=0.02d0 + amro=1d0 + a4=4d0*ampi + alo=dlog(amro/ampi) + bpol=2.4d0 + + nt=6000 + htt=6d0/dble(nt)**2 + + out=0d0 + out1=0d0 + + do it=0,nt + t=dble(it)**2*htt + if(it.eq.0) t=1d-8 + wt=htt*2d0*dble(it)/4d0/pi + if(it.eq.0) wt=wt/2d0 + bes0=besj0(bt*dsqrt(t)) + + ffi=dexp(-((t+0.08d0+bb0(i))*bex(i))**cc0(i)+ + & (bex(i)*(bb0(i)+0.08d0))**cc0(i)) + ffj=dexp(-((t+0.08d0+bb0(j))*bex(j))**cc0(j)+ + & (bex(j)*(bb0(j)+0.08d0))**cc0(j)) + + asp1=asp + form1=dlog(ffi*ffj)-2d0*t*asp1*dlog(rts) +cccccccccccccccc + ah=dsqrt(1d0+a4/t) + h1pi=2d0*a4+t*(alo-ah*ah*ah*dlog((1d0+ah)/(ah-1d0))) !!! Pion loop insertion + h1pi=h1pi*sigo/(72d0*pi**3*(1d0+t/bpol)**2) +ccccccccccccccc + ww=bes0*dexp(form1-2d0*h1pi*dlog(rts)) + aspt=t*asp+h1pi + + out=out+ww*wt + out1=out1+bes0*wt*ffi*ffj + + enddo + + + return + end + + function rgauss(r1,r2) + implicit double precision (a-y) + + pi=dacos(-1d0) + rgauss=dsqrt(-2d0*dlog(r1))*dsin(2d0*pi*r2) + + return + end + +c boosting subroutine + + subroutine boost(p,pboo,pcm,plb) + real*8 pboo(4),pcm(4),plb(4),p,fact + plb(4)=(pboo(4)*pcm(4)+pboo(3)*pcm(3) + & +pboo(2)*pcm(2)+pboo(1)*pcm(1))/p + fact=(plb(4)+pcm(4))/(p+pboo(4)) + do 10 j=1,3 + 10 plb(j)=pcm(j)+fact*pboo(j) + return + end + +c calculates lorentzian dot product for real 4-vectors + + subroutine dot(v1,v2,dt) + double precision v1(4),v2(4),dt + + dt=-(v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3) + &-v1(4)*v2(4)) + + return + end +c + +* +* subtractive mitchell-moore generator +* ronald kleiss - october 2, 1987 +* +* the algorithm is N(i)=[ N(i-24) - N(i-55) ]mod M, +* implemented in a cirucular array with identifcation +* of NR(i+55) and nr(i), such that effectively: +* N(1) <--- N(32) - N(1) +* N(2) <--- N(33) - N(2) .... +* .... N(24) <--- N(55) - N(24) +* N(25) <--- N(1) - N(25) .... +* .... N(54) <--- N(30) - N(54) +* N(55) <--- N(31) - N(55) +* +* in this version M =2**30 and RM=1/M=2.D0**(-30.D0) +* +* the array NR has been initialized by putting NR(i)=i +* and subsequently running the algorithm 100,000 times. +* + + subroutine R2455(Ran) + IMPLICIT REAL*8(A-H,O-Z) + DIMENSION N(55) + DATA N/ + . 980629335, 889272121, 422278310,1042669295, 531256381, + . 335028099, 47160432, 788808135, 660624592, 793263632, + . 998900570, 470796980, 327436767, 287473989, 119515078, + . 575143087, 922274831, 21914605, 923291707, 753782759, + . 254480986, 816423843, 931542684, 993691006, 343157264, + . 272972469, 733687879, 468941742, 444207473, 896089285, + . 629371118, 892845902, 163581912, 861580190, 85601059, + . 899226806, 438711780, 921057966, 794646776, 417139730, + . 343610085, 737162282,1024718389, 65196680, 954338580, + . 642649958, 240238978, 722544540, 281483031,1024570269, + . 602730138, 915220349, 651571385, 405259519, 145115737/ + DATA M/1073741824/ + DATA RM/0.9313225746154785D-09/ + DATA K/55/,L/31/ + IF(K.EQ.55) THEN + K=1 + ELSE + K=K+1 + ENDIF + IF(L.EQ.55) THEN + L=1 + ELSE + L=L+1 + ENDIF + J=N(L)-N(K) + IF(J.LT.0) J=J+M + N(K)=J + RAN=J*RM + END + + + + SUBROUTINE RAMBO(N,ET,XM,P,WT) +C------------------------------------------------------ +C +C RAMBO +C +C RA(NDOM) M(OMENTA) B(EAUTIFULLY) O(RGANIZED) +C +C A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR +C AUTHORS: S.D. ELLIS, R. KLEISS, W.J. STIRLING +C THIS IS VERSION 1.0 - WRITTEN BY R. KLEISS +C +C N = NUMBER OF PARTICLES (>1, IN THIS VERSION <101) +C ET = TOTAL CENTRE-OF-MASS ENERGY +C XM = PARTICLE MASSES ( DIM=100 ) +C P = PARTICLE MOMENTA ( DIM=(4,100) ) +C WT = WEIGHT OF THE EVENT +C +C------------------------------------------------------ + implicit none +! boost arguments + integer n + double precision et,xm(100),p(4,100),wt +! function + double precision rn +! internal variables + double precision q(4,100),z(100),r(4), + & b(3),p2(100),xm2(100),e(100),v(100) + integer iwarn(5) + double precision acc + integer itmax,ibegin + data acc/1.0d-14/,itmax/6/,ibegin/0/,iwarn/5*0/ + + integer i,k,iter + double precision twopi,po2log,xmt,c,s,f,rmas,g,a,x,bq, + & accu,f0,g0,wt2,wt3,wtm,x2,xmax + integer nm + save + + +C +C INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT + IF(IBEGIN.NE.0) GOTO 103 + IBEGIN=1 + TWOPI=8.*DATAN(1.D0) + PO2LOG=DLOG(TWOPI/4.) + Z(2)=PO2LOG + DO 101 K=3,100 + 101 Z(K)=Z(K-1)+PO2LOG-2.*DLOG(DFLOAT(K-2)) + DO 102 K=3,100 + 102 Z(K)=(Z(K)-DLOG(DFLOAT(K-1))) +C +C CHECK ON THE NUMBER OF PARTICLES + 103 IF(N.GT.1.AND.N.LT.101) GOTO 104 + PRINT 1001,N + STOP +C +C CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES + 104 XMT=0. + NM=0 + DO 105 I=1,N + IF(XM(I).NE.0.D0) NM=NM+1 + 105 XMT=XMT+DABS(XM(I)) + IF(XMT.LE.ET) GOTO 201 + PRINT 1002,XMT,ET + STOP +C +C THE PARAMETER VALUES ARE NOW ACCEPTED +C +C GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE + 201 DO 202 I=1,N + C=2.*RN(1)-1. + S=DSQRT(1.-C*C) + F=TWOPI*RN(2) + Q(4,I)=-DLOG(RN(3)*RN(4)) + Q(3,I)=Q(4,I)*C + Q(2,I)=Q(4,I)*S*DCOS(F) + 202 Q(1,I)=Q(4,I)*S*DSIN(F) + +C +C CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION + DO 203 I=1,4 + 203 R(I)=0. + DO 204 I=1,N + DO 204 K=1,4 + 204 R(K)=R(K)+Q(K,I) + RMAS=DSQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2) + DO 205 K=1,3 + 205 B(K)=-R(K)/RMAS + G=R(4)/RMAS + A=1./(1.+G) + X=ET/RMAS +C +C TRANSFORM THE Q'S CONFORMALLY INTO THE P'S + DO 207 I=1,N + BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I) + DO 206 K=1,3 + 206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ)) + 207 P(4,I)=X*(G*Q(4,I)+BQ) +C +C CALCULATE WEIGHT AND POSSIBLE WARNINGS + WT=PO2LOG + IF(N.NE.2) WT=(2.*N-4.)*DLOG(ET)+Z(N) + IF(WT.GE.-180.D0) GOTO 208 + IF(IWARN(1).LE.5) PRINT 1004,WT + IWARN(1)=IWARN(1)+1 + 208 IF(WT.LE. 174.D0) GOTO 209 + IF(IWARN(2).LE.5) PRINT 1005,WT + IWARN(2)=IWARN(2)+1 +C +C RETURN FOR WEIGHTED MASSLESS MOMENTA + 209 IF(NM.NE.0) GOTO 210 + WT=DEXP(WT) + RETURN +C +C MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X + 210 XMAX=DSQRT(1.-(XMT/ET)**2) + DO 301 I=1,N + XM2(I)=XM(I)**2 + 301 P2(I)=P(4,I)**2 + ITER=0 + X=XMAX + ACCU=ET*ACC + 302 F0=-ET + G0=0. + X2=X*X + DO 303 I=1,N + E(I)=DSQRT(XM2(I)+X2*P2(I)) + F0=F0+E(I) + 303 G0=G0+P2(I)/E(I) + IF(DABS(F0).LE.ACCU) GOTO 305 + ITER=ITER+1 + IF(ITER.LE.ITMAX) GOTO 304 +C PRINT 1006,ITMAX,ACCU,DABS(F0) + WRITE(99,1006)ITMAX,ACCU,DABS(F0) + IF (DABS(F0).GT.1.0E-6) THEN + WRITE(*,1007)DABS(F0) + ENDIF + GOTO 305 + 304 X=X-F0/(X*G0) + GOTO 302 + 305 DO 307 I=1,N + V(I)=X*P(4,I) + DO 306 K=1,3 + 306 P(K,I)=X*P(K,I) + 307 P(4,I)=E(I) +C +C CALCULATE THE MASS-EFFECT WEIGHT FACTOR + WT2=1. + WT3=0. + DO 308 I=1,N + WT2=WT2*V(I)/E(I) + 308 WT3=WT3+V(I)**2/E(I) + WTM=(2.*N-3.)*DLOG(X)+DLOG(WT2/WT3*ET) +C +C RETURN FOR WEIGHTED MASSIVE MOMENTA + WT=WT+WTM + IF(WT.GE.-180.D0) GOTO 309 + IF(IWARN(3).LE.5) PRINT 1004,WT + IWARN(3)=IWARN(3)+1 + 309 IF(WT.LE. 174.D0) GOTO 310 + IF(IWARN(4).LE.5) PRINT 1005,WT + IWARN(4)=IWARN(4)+1 + 310 WT=DEXP(WT) + RETURN +C + 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED') + 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT', + . ' SMALLER THAN TOTAL ENERGY =',D15.6) + 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW') + 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY OVERFLOW') + 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE', + . ' DESIRED ACCURACY =',D10.3,' (STOPS AT',D10.3,')') + 1007 FORMAT(' RAMBO WARNS: END OF ITERATION ACCURACY TOO LOW =',D10.3) + END +C + + function rn(idum) +* +* SUBTRACTIVE MITCHELL-MOORE GENERATOR +* RONALD KLEISS - OCTOBER 2, 1987 +* +* THE ALGORITHM IS N(I)=[ N(I-24) - N(I-55) ]MOD M, +* IMPLEMENTED IN A CIRUCULAR ARRAY WITH IDENTIFCATION +* OF NR(I+55) AND NR(I), SUCH THAT EFFECTIVELY: +* N(1) <--- N(32) - N(1) +* N(2) <--- N(33) - N(2) .... +* .... N(24) <--- N(55) - N(24) +* N(25) <--- N(1) - N(25) .... +* .... N(54) <--- N(30) - N(54) +* N(55) <--- N(31) - N(55) +* +* IN THIS VERSION M =2**30 AND RM=1/M=2.D0**(-30D0) +* +* THE ARRAY NR HAS BEEN INITIALIZED BY PUTTING NR(I)=I +* AND SUBSEQUENTLY RUNNING THE ALGORITHM 100,000 TIMES. +* + + implicit none + double precision rn + integer idum + integer n(55) + data n/ + . 980629335, 889272121, 422278310,1042669295, 531256381, + . 335028099, 47160432, 788808135, 660624592, 793263632, + . 998900570, 470796980, 327436767, 287473989, 119515078, + . 575143087, 922274831, 21914605, 923291707, 753782759, + . 254480986, 816423843, 931542684, 993691006, 343157264, + . 272972469, 733687879, 468941742, 444207473, 896089285, + . 629371118, 892845902, 163581912, 861580190, 85601059, + . 899226806, 438711780, 921057966, 794646776, 417139730, + . 343610085, 737162282,1024718389, 65196680, 954338580, + . 642649958, 240238978, 722544540, 281483031,1024570269, + . 602730138, 915220349, 651571385, 405259519, 145115737/ + double precision eps + double precision rm + integer j,k,l,m + + data eps/1D-9/ + data m/1073741824/ + data rm/0.9313225746154785D-09/ + data k/55/,l/31/ + + + 1 CONTINUE + IF(K.EQ.55) THEN + K=1 + ELSE + K=K+1 + ENDIF + IF(L.EQ.55) THEN + L=1 + ELSE + L=L+1 + ENDIF + J=N(L)-N(K) + IF(J.LT.0) J=J+M + N(K)=J + RN=J*RM + IF(RN.LT.EPS) GOTO 1 + IF(RN.GT.1D0-EPS) GOTO 1 + RETURN + END -- 2.43.0