diff --git a/CHANGELOG b/CHANGELOG new file mode 100644 index 0000000..e96d3ea --- /dev/null +++ b/CHANGELOG @@ -0,0 +1,105 @@ +2016-10-07 Cristian Baldenegro + + * Added EFT for AA->AZ, according to Sylvain and Gero's computations. Datacard: dataQED_AZanom + * In 2017-02-15 corrected the ME to include two operators which induce the Zgamma process. + +2016-09-23 Cristian Baldenegro + + * Correction for AA->Spin2->AA (-1 factor when calling Mxxxx matrix elements) + * Couldn't run dataQED_ZZanom anymore; now it's corrected. + +2016-06-20 Cristian Baldenegro + + * Exclusive AA-> HH (IPROC=16071) and AA->Dijet (IPROC=16070) were included + * Matrix elements calculated by Gero von Gersdorff and Sylvain Fichet + * Correspond to Datacards: dataQED_HHexoticSPin0EvenResonances and dataQED_GluGluexoticSpin0EvenResonances respectively + * New parameters (F0H and F0G) + +2016-05-20 Cristian Baldenegro + + * Exclusive AA->ZZ (IPROC=16066) and AA->AZ (IPROC=16068) + * Exclusive AA->WW (IPROC=16067) and AA->AA (IPROC=16069) processes included + * Helicity amplitudes from Gero von Gersdorff and Sylvain Fichet + * Corresponding datacards: dataQED_ZZexoticSpin0EvenResonances and dataQED_AZexoticSpin0EvenResonances +dataQED_WWexoticSpin0EvenResonances and dataQED_AAexoticSpin2Resonances + * New parameter F0Z (1/coupling for Phi->ZZ) AF0ZG (1/coupling for Phi->gamma Z) + * New parameter F0W (1/coupling for phi->WW) + * Such exclusive productions are included in the same subroutine calling section as Matthias' (AAANOM=3, AAEXOTIC=1) + +2015-07-24 Matthias Saimpert + * IMPORTANT BUG CORRECTIONS: symmetry factor was omitted for + gamgam->gamgam, now added. As a result, gamgam->gamgam xsec twice smaller + +2015-06-25 Matthias Saimpert + * exclusive gamma gamma->gamma gamma studies, exotic spin 0-even resonance added + * IPROC=16065, associated datacard, dataQED_AAexoticSpin0EvenResonances + * new parameters AAF0 (1/coupling), AAW (width const), AAA2 (dumping factor + for sqrt(s) dependence of the width) + +2015-01-20 Christophe Royon , Murilo Rangel + * correct for normalisation for reggeon (applied in PDFs and not in flux) + +2014-11-26 Christophe Royon + * new code for heavy ion flux (nflux=11) + * new cards to transmit BMIN + * new code and fluxes for proton ion, ion proton via Pom/Photon exchange + * bug fix for photon pom + +2014-09-11 Matthias Saimpert + * Fix for gamma gamma->gamma gamma, high energy limit implemented (sqrt(s)>>m) + * Final release for gamma gamma->gamma gamma :) + +2014-08-18 Matthias Saimpert + * Anomalous aaaa with compHEP routine found to be wrong (around a factor 2.1 + too small) + * Code from Gero and Sylvain replaced the compHEP routine, everything is now + in agreement + * Possibility to switch back to compHEP routine by commenting/uncommenting + fpmc.f l. 3435 + * The new EFT routine takes interferences with SM aa production via W loop + into account + +2014-03-20 Matthias Saimpert + * General update/upgrade of the QED Diphoton exclusive process + * IPROC=16060, SM QED Exc diphoton production (W + massive fermions loops included) + * (old routine is IPROC=19800 kept for comparison purpose but is most likely + to be WRONG) + * IPROC=16061, SM QED Exc diphoton production at high mass (massive fermion + loops EXCLUDED) + * IPROC=16062, SM QED Exc diphoton production at low mass (W loop EXCLUDED) + * IPROC=16063, General exotic vector contributions with SM interf (SM fermion + loops EXCLUDED + * IPROC=16064, General exotic fermion contributions with SM interf (SM + fermion EXCLUDED) + * Corresponding datacards added + +2014-01-21 Christophe Royon + * with Goergiy Stelmakh, + * add Reggeon (nflux=10), Pomeron-Reggeon (nflux=19), Reggeon Pomeron (nflux=21) + * Photon-Pomeron (nflux=20), Pomeron-Photon (nflux=22) + +2013-12-12 Matthias Saimpert + * Adding aa->aa anomalous couplings (IPROC=16016) based on arXiv:1311.6815 [hep-ph] + * Adding ~/External/comphep_interface/sqme_aaaa folder + * Modifications of ~/Examples/ffcard., module., module_reco., fpmc_main. and ~/Fpmc/fpmc. + +2013-03-26 Oldrich Kepka + * Adding setup_lxplus.sh, to be used on lxplus to configure CERNLIB + +2012-06-05 Rafal Staszewski + * New Examples/fpmc_main.f created to supercede module.f and module_reco.f. + module and module_reco are disabled in the Makefile. The new exacutable is + called fpmc. + * The information from jet alg. was not stored into ntuple. A change in + External/ntuple.f was needed to fix it. + * New OUTPUT option added to steering, determining the generator output. + Possible values: + 0 : Only the text output with the cross section at the end. + 1 : Ntuple with few histograms and a list of particles for each event. + 2 : Jet CONE algorithm is performed and its output is also stored in the + ntuple. + +2011-05-09 Oldrich Kepka + * Tagging version fpmc_v01-05 - cleanup of nutples, version for testing + + diff --git a/CMakeLists.txt b/CMakeLists.txt index af2e663..3f06b72 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,6 +3,7 @@ project(FPMC) set(PROJECT_VERSION 1) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wno-long-long -Wno-deprecated-declarations -pedantic -lgfortran -std=c++11 -g") +set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -g -O1 -fno-automatic -fPIC") set(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake") set(MODULES Fpmc Herwig) diff --git a/Examples/FPMC.INC b/Examples/FPMC.INC deleted file mode 100644 index afc31f3..0000000 --- a/Examples/FPMC.INC +++ /dev/null @@ -1 +0,0 @@ - INCLUDE '../Fpmc/fpmc.inc' \ No newline at end of file diff --git a/Examples/HERWIG65.INC b/Examples/HERWIG65.INC deleted file mode 100644 index 0bbb0b7..0000000 --- a/Examples/HERWIG65.INC +++ /dev/null @@ -1,2 +0,0 @@ -* INCLUDE 'herwig6510.inc' - INCLUDE 'herwig6500.inc' diff --git a/Examples/fpmc-bare.cc b/Examples/fpmc-bare.cc new file mode 100644 index 0000000..94dae0c --- /dev/null +++ b/Examples/fpmc-bare.cc @@ -0,0 +1,19 @@ +#include "Fpmc/Fpmc.h" + +#include + +using namespace std; + +int main( int argc, char* argv[] ) +{ + fpmc::Fpmc generator( argv[1] ); + + generator.parameters().dump(); + + hwevnt_t ev; + for( unsigned int evt = 0; evt < 10; ++evt ) { + cout << "[FPMC] Processing event " << ( evt+1 ) << endl; + generator.next( ev ); + } + return 0; +} diff --git a/Examples/fpmc-wrapper.cc b/Examples/fpmc-wrapper.cc new file mode 100644 index 0000000..075ed64 --- /dev/null +++ b/Examples/fpmc-wrapper.cc @@ -0,0 +1,57 @@ +#include "Fpmc/HepMCWrapper.h" +#include "Fpmc/ArgsParser.h" + +#include +#include + +#ifdef HEPMC_VERSION2 +#include "HepMC/IO_GenEvent.h" +#else +#include "HepMC/WriterAscii.h" +#endif + +using namespace std; + +int main( int argc, char* argv[] ) +{ + fpmc::ArgsParser args( argc, argv, { "cfg", "nevents", "comenergy" }, { "fileout" } ); + const auto required_params = args.required_parameters(), optional_params = args.optional_parameters(); + + //---------------------------- + // Required parameters + string datacard_ = required_params.at( "cfg" ); + unsigned int maxEvents_ = stoi( required_params.at( "nevents" ) ); + double comEnergy_ = stof( required_params.at( "comenergy" ) ); + + // Optional parameters + string outputFileName_ = "fpmc.hepmc"; + if ( optional_params.count("fileout") != 0 ) outputFileName_ = optional_params.at( "fileout" ); + + stringstream oss; + oss << "=========================================================" << endl + << "FPMC (Wrapper) will initialize with parameters: " << endl + << " Datacard: " << datacard_ << endl + << " N events: " << maxEvents_ << endl + << " COM energy: " << comEnergy_ << endl + << " Output file: " << outputFileName_ << endl + << "=========================================================" << endl; + cout << oss.str(); + + fpmc::HepMCWrapper generator( comEnergy_, datacard_.c_str() ); + + +#ifdef HEPMC_VERSION2 + HepMC::IO_GenEvent output( outputFileName_, ios::out ); +#else + HepMC::WriterAscii output( outputFileName_ ); +#endif + for( unsigned int evt = 0; evt < maxEvents_; ++evt ) { + cout << "[FPMC Wrapper] Processing event " << ( evt+1 ) << endl; +#ifdef HEPMC_VERSION2 + output.write_event( generator.event() ); +#else + output.write_event( *generator.event() ); +#endif + } + return 0; +} diff --git a/Examples/herwig6500.inc b/Examples/herwig6500.inc deleted file mode 100644 index 70ff5d7..0000000 --- a/Examples/herwig6500.inc +++ /dev/null @@ -1,344 +0,0 @@ -C ****COMMON BLOCK FILE FOR HERWIG VERSION 6.4**** -C -C ALTERATIONS: Layout completely overhauled for 5.9 -C -C -C New common blocks added for version 6.1: -C HWCLUS,HWSUSY,HWRPAR,HWMINB -C -C New variables added for version 6.1: -C OMHMIX,ET2MIX,PH3MIX,IOP4JT,NPRFMT, -C PRNDEF,PRNTEX,PRNWEB,EFFMIN,GCUTME, -C IOP4JT,NPRFMT see HWPRAM -C Y4JT,DURHAM see HWHARD -C QORQQB,QBORQQ see HWPROP -C NRECO see HWUCLU -C TXNAME see HWUNAM -C PPCL,NCL,IDCL see HWCLUS -C TANB,ALPHAH,COSBPA,SINBPA,COSBMA, -C SINBMA,COSA,SINA,COSB,SINB,COTB, -C ZMIXSS,ZMXNSS,ZSGNSS,LFCH,RFCH, -C SLFCH,SRFCH,WMXUSS,WMXVSS,WSGNSS, -C QMIXSS,LMIXSS,THETAT,THETAB,THETAL, -C ATSS,ABSS,ALSS,MUSS,FACTSS,GHWWSS, -C GHZZSS,GHDDSS,GHUUSS,GHWHSS,GHSQSS, -C XLMNSS,RMMNSS,DMSSM,SENHNC, -C SSPARITY,SUSYIN see HWSUSY -C LAMDA1,LAMDA2,LAMDA3,HRDCOL,RPARTY, -C COLUPD see HWRPAR -C PMBN1,PMBN2,PMBN3,PMBK1,PMBK2, -C PMBM1,PMBM2,PMBP1,PMBP2,PMBP3 see HWMINB -C -C New parameters added for version 6.1: -C NMXCL -C -C Parameter NMXRES raised to 500 -C -C Scalar variables changed to arrays of size 2: -C CLSMR,PSPLT,CLDIR see HWPRAM -C -C NEW for HERWIG6.200 -C -C New common blocks added for version 6.2: -C HWGRAV see HWHGRV -C -C NEW for HERWIG6.202 -C -C New common block added for version 6.202: -C HW6202 -C which contains all other new variables, parameters and -C control flags introduced since version 6.1, so that -C other common blocks become identical to version 6.1 -C -C New parameters added for version 6.2: -C VIPWID,DXRCYL,DXZMAX,DXRSPH,LRSUSY see HWIGIN -C GRVLAM,EMGRV,GAMGRV see HWHGRV -C -C New control flags added for version 6.2: -C WZRFR see HWBJCO -C FIX4JT see HWIGIN -C IMSSM,IHIGGS see HWUINC -C -C New variable added for version 6.2: -C PARITY see HWUINC -C -C Version 6.203: -C -C NMXHEP raised to 4000 for LHC studies -C -C - IMPLICIT NONE - DOUBLE PRECISION ZERO,ONE,TWO,THREE,FOUR,HALF - PARAMETER (ZERO =0.D0, ONE =1.D0, TWO =2.D0, - & THREE=3.D0, FOUR=4.D0, HALF=0.5D0) -C - DOUBLE PRECISION - & ACCUR,AFCH,ALPFAC,ALPHEM,ANOMSC,ASFIXD,AVWGT,B1LIM,BETAF,BRFRAC, - & BRHIG,BTCLM,CAFAC,CFFAC,CLDKWT,CLMAX,CLPOW,CLQ,CLSMR,CMMOM,COSS, - & COSTH,CSPEED,CTHRPW,CTMAX,DECPAR,DECWT,DISF,DKLTM,EBEAM1,EBEAM2, - & EMLST,EMMAX,EMMIN,EMPOW,EMSCA,ENHANC,ENSOF,EPOLN,ETAMIX,EVWGT, - & EXAG,F0MIX,F1MIX,F2MIX,GAMH,GAMMAX,GAMW,GAMWT,GAMZ,GAMZP,GCOEF, - & GEV2NB,GEV2MM,GPOLN,H1MIX,HBAR,HARDST,OMEGA0,PBEAM1,PBEAM2,PDIQK, - & PGSMX,PGSPL,PHEP,PHIMIX,PHIPAR,PHOMAS,PIFAC,PLTCUT,PPAR,PPOLN, - & PRECO,PRSOF,PSPLT,PTINT,PTMAX,PTMIN,PTPOW,PTRMS,PXRMS,PWT,Q2MAX, - & Q2MIN,Q2POW,Q2WWMN,Q2WWMX,QCDL3,QCDL5,QCDLAM,QDIQK,QEV,QFCH,QG, - & QLIM,QSPAC,QV,QWT,REPWT,RESN,RHOHEP,RHOPAR,RLTIM,RMASS,RMIN, - & RSPIN,SCABI,SINS,SNGWT,SWEIN,SWTEF,SUD,THMAX,TLOUT,TMTOP,TMNISR, - & TQWT,VCKM,VFCH,VGCUT,VHEP,VMIN2,VPAR,VPCUT,VQCUT,VTXPIP,VTXQDK, - & WBIGST,WGTMAX,WGTSUM,WHMIN,WSQSUM,XFACT,XLMIN,XMIX,XMRCT,XX, - & XXMIN,YBMAX,YBMIN,YJMAX,YJMIN,YMIX,YMRCT,YWWMAX,YWWMIN,ZBINM, - & ZJMAX,ZMXISR,Y4JT,EFFMIN,PPCL, - & TANB,ALPHAH,COSBPA,SINBPA,COSBMA,SINBMA,COSA,SINA,COSB,SINB,COTB, - & ZMIXSS,ZMXNSS,ZSGNSS,LFCH,RFCH,SLFCH,SRFCH, WMXUSS,WMXVSS,WSGNSS, - & QMIXSS,LMIXSS,THETAT,THETAB,THETAL,ATSS,ABSS,ALSS,MUSS,FACTSS, - & GHWWSS,GHZZSS,GHDDSS,GHUUSS,GHWHSS,GHSQSS -C--fix by PR 12/7/02 to avoid problems with nag compiler - DOUBLE PRECISION - & XLMNSS,RMMNSS,DMSSM,SENHNC,SSPARITY,LAMDA1,LAMDA2,LAMDA3, - & PMBN1,PMBN2,PMBN3,PMBK1,PMBK2,PMBM1,PMBM2,PMBP1,PMBP2,PMBP3, - & OMHMIX,ET2MIX,PH3MIX,GCUTME -C - INTEGER - & CLDIR,IAPHIG,IBRN,IBSH,ICHRG,ICO,IDCMF,IDHEP,IDHW,IDK,IDKPRD,IDN, - & IDPAR,IDPDG,IERROR,IFLAV,IFLMAX,IFLMIN,IHPRO,IMQDK,INHAD,INTER, - & IOPDKL,IOPHIG,IOPREM,IPART1,IPART2,IPRINT,IPRO,IPROC,ISLENT, - & ISPAC,ISTAT,ISTHEP,ISTPAR,JCOPAR,JDAHEP,JDAPAR,JMOHEP,JMOPAR, - & JNHAD,LNEXT,LOCN,LOCQ,LRSUD,LSTRT,LWEVT,LWSUD,MAPQ,MAXER,MAXEV, - & MAXFL,MAXPR,MODBOS,MODMAX,MODPDF,NBTRY,NCLDK,NCOLO,NCTRY,NDKYS, - & NDTRY,NETRY,NEVHEP,NEVPAR,NFLAV,NGSPL,NHEP,NME,NMODES,NMXCDK, - & NMXDKS,NMXHEP,NMXJET,NMXMOD,NMXPAR,NMXQDK,NMXRES,NMXSUD,NPAR, - & NPRODS,NQDK,NQEV,NRES,NRN,NSPAC,NSTRU,NSTRY,NSUD,NUMER,NUMERU, - & NWGTS,NZBIN,SUDORD,IOP4JT,HRDCOL,NMXCL,NCL,IDCL,NPRFMT,NRECO -C - LOGICAL - & AZSOFT,AZSPIN,BGSHAT,BREIT,CLRECO,COLISR,DKPSET,FROST,FSTEVT, - & FSTWGT,GENEV,GENSOF,HARDME,HVFCEN,MAXDKL,MIXING,NOSPAC,NOWGT, - & PRNDEC,PIPSMR,PRVTX,RSTAB,SOFTME,TMPAR,TPOL,USECMF,VTOCDK,VTORDK, - & ZPRIME,RPARTY,COLUPD,PRNDEF,PRNTEX,PRNWEB,DURHAM,SUSYIN, - & QORQQB,QBORQQ -C - CHARACTER*4 - & BDECAY - CHARACTER*8 - & PART1,PART2,RNAME - CHARACTER*20 - & AUTPDF - CHARACTER*37 - & TXNAME -C -C New standard event common - PARAMETER (NMXHEP=4000) - COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP), - & JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP) -C -C Beams, process and number of events - COMMON/HWBEAM/IPART1,IPART2 - COMMON/HWBMCH/PART1,PART2 - COMMON/HWPROC/EBEAM1,EBEAM2,PBEAM1,PBEAM2,IPROC,MAXEV -C -C Basic parameters (and quantities derived from them) - COMMON/HWPRAM/AFCH(16,2),ALPHEM,B1LIM,BETAF,BTCLM,CAFAC,CFFAC, - & CLMAX,CLPOW,CLSMR(2),CSPEED,ENSOF,ETAMIX,F0MIX,F1MIX,F2MIX,GAMH, - & GAMW,GAMZ,GAMZP,GEV2NB,H1MIX,PDIQK,PGSMX,PGSPL(4),PHIMIX,PIFAC, - & PRSOF,PSPLT(2),PTRMS,PXRMS,QCDL3,QCDL5,QCDLAM,QDIQK,QFCH(16),QG, - & QSPAC,QV,SCABI,SWEIN,TMTOP,VFCH(16,2),VCKM(3,3),VGCUT,VQCUT, - & VPCUT,ZBINM,EFFMIN,OMHMIX,ET2MIX,PH3MIX,GCUTME, - & IOPREM,IPRINT,ISPAC,LRSUD,LWSUD,MODPDF(2),NBTRY,NCOLO,NCTRY, - & NDTRY,NETRY,NFLAV,NGSPL,NSTRU,NSTRY,NZBIN,IOP4JT(2),NPRFMT, - & AZSOFT,AZSPIN,CLDIR(2),HARDME,NOSPAC,PRNDEC,PRVTX,SOFTME,ZPRIME, - & PRNDEF,PRNTEX,PRNWEB -C - COMMON/HWPRCH/AUTPDF(2),BDECAY -C -C Parton shower common (same format as /HEPEVT/) - PARAMETER (NMXPAR=500) - COMMON/HWPART/NEVPAR,NPAR,ISTPAR(NMXPAR),IDPAR(NMXPAR), - & JMOPAR(2,NMXPAR),JDAPAR(2,NMXPAR),PPAR(5,NMXPAR),VPAR(4,NMXPAR) -C -C Parton polarization common - COMMON/HWPARP/DECPAR(2,NMXPAR),PHIPAR(2,NMXPAR),RHOPAR(2,NMXPAR), - & TMPAR(NMXPAR) -C -C Electroweak boson common - PARAMETER (MODMAX=5) - COMMON/HWBOSC/ALPFAC,BRHIG(12),ENHANC(12),GAMMAX,RHOHEP(3,NMXHEP), - & IOPHIG,MODBOS(MODMAX) -C -C Parton colour common - COMMON/HWPARC/JCOPAR(4,NMXPAR) -C -C other HERWIG branching, event and hard subprocess common blocks - COMMON/HWBRCH/ANOMSC(2,2),HARDST,PTINT(3,2),XFACT,INHAD,JNHAD, - & NSPAC(7),ISLENT,BREIT,FROST,USECMF -C - COMMON/HWEVNT/AVWGT,EVWGT,GAMWT,TLOUT,WBIGST,WGTMAX,WGTSUM,WSQSUM, - & IDHW(NMXHEP),IERROR,ISTAT,LWEVT,MAXER,MAXPR,NOWGT,NRN(2),NUMER, - & NUMERU,NWGTS,GENSOF -C - COMMON/HWHARD/ASFIXD,CLQ(7,6),COSS,COSTH,CTMAX,DISF(13,2),EMLST, - & EMMAX,EMMIN,EMPOW,EMSCA,EPOLN(3),GCOEF(7),GPOLN,OMEGA0,PHOMAS, - & PPOLN(3),PTMAX,PTMIN,PTPOW,Q2MAX,Q2MIN,Q2POW,Q2WWMN,Q2WWMX,QLIM, - & SINS,THMAX,Y4JT,TMNISR,TQWT,XX(2),XLMIN,XXMIN,YBMAX,YBMIN,YJMAX, - & YJMIN,YWWMAX,YWWMIN,WHMIN,ZJMAX,ZMXISR,IAPHIG,IBRN(2),IBSH, - & ICO(10),IDCMF,IDN(10),IFLMAX,IFLMIN,IHPRO,IPRO,MAPQ(6),MAXFL, - & BGSHAT,COLISR,FSTEVT,FSTWGT,GENEV,HVFCEN,TPOL,DURHAM -C -C Arrays for particle properties (NMXRES = max no of particles defined) - PARAMETER(NMXRES=500) - COMMON/HWPROP/RLTIM(0:NMXRES),RMASS(0:NMXRES),RSPIN(0:NMXRES), - & ICHRG(0:NMXRES),IDPDG(0:NMXRES),IFLAV(0:NMXRES),NRES, - & VTOCDK(0:NMXRES),VTORDK(0:NMXRES), - & QORQQB(0:NMXRES),QBORQQ(0:NMXRES) -C - COMMON/HWUNAM/RNAME(0:NMXRES),TXNAME(2,0:NMXRES) -C -C Arrays for particle decays (NMXDKS = max total no of decays, -C NMXMOD = max no of modes for a particle) - PARAMETER(NMXDKS=4000,NMXMOD=200) - COMMON/HWUPDT/BRFRAC(NMXDKS),CMMOM(NMXDKS),DKLTM(NMXRES), - & IDK(NMXDKS),IDKPRD(5,NMXDKS),LNEXT(NMXDKS),LSTRT(NMXRES),NDKYS, - & NME(NMXDKS),NMODES(NMXRES),NPRODS(NMXDKS),DKPSET,RSTAB(0:NMXRES) -C -C Weights used in cluster decays - COMMON/HWUWTS/REPWT(0:3,0:4,0:4),SNGWT,DECWT,QWT(3),PWT(12), - & SWTEF(NMXRES) -C -C Parameters for cluster decays (NMXCDK = max total no of cluster -C decay channels) - PARAMETER(NMXCDK=4000) - COMMON/HWUCLU/CLDKWT(NMXCDK),CTHRPW(12,12),PRECO,RESN(12,12), - & RMIN(12,12),LOCN(12,12),NCLDK(NMXCDK),NRECO,CLRECO -C -C Variables controling mixing and vertex information -C--VTXPIP should have been a 5-vector, problems with NAG compiler - COMMON/HWDIST/EXAG,GEV2MM,HBAR,PLTCUT,VMIN2,VTXPIP(5),XMIX(2), - & XMRCT(2),YMIX(2),YMRCT(2),IOPDKL,MAXDKL,MIXING,PIPSMR -C -C Arrays for temporarily storing heavy-b,c-hadrons decaying partonicaly -C (NMXQDK = max no such decays in an event) - PARAMETER (NMXQDK=20) - COMMON/HWQDKS/VTXQDK(4,NMXQDK),IMQDK(NMXQDK),LOCQ(NMXQDK),NQDK -C -C Parameters for Sudakov form factors -C (NMXSUD= max no of entries in lookup table) - PARAMETER (NMXSUD=1024) - COMMON/HWUSUD/ACCUR,QEV(NMXSUD,6),SUD(NMXSUD,6),INTER,NQEV,NSUD, - & SUDORD -C - PARAMETER (NMXJET=200) -C -C SUSY parameters - COMMON/HWSUSY/ - & TANB,ALPHAH,COSBPA,SINBPA,COSBMA,SINBMA,COSA,SINA,COSB,SINB,COTB, - & ZMIXSS(4,4),ZMXNSS(4,4),ZSGNSS(4), LFCH(16),RFCH(16), - & SLFCH(16,4),SRFCH(16,4), WMXUSS(2,2),WMXVSS(2,2), WSGNSS(2), - & QMIXSS(6,2,2),LMIXSS(6,2,2), - & THETAT,THETAB,THETAL,ATSS,ABSS,ALSS,MUSS,FACTSS, - & GHWWSS(3),GHZZSS(3),GHDDSS(4),GHUUSS(4),GHWHSS(3), - & GHSQSS(4,6,2,2),XLMNSS,RMMNSS,DMSSM,SENHNC(24),SSPARITY,SUSYIN -C -C R-Parity violating parameters and colours - COMMON /HWRPAR/ LAMDA1(3,3,3),LAMDA2(3,3,3), - & LAMDA3(3,3,3),HRDCOL(2,5),RPARTY,COLUPD -C -C Parameters for minimum bias/soft underlying event - COMMON/HWMINB/ - & PMBN1,PMBN2,PMBN3,PMBK1,PMBK2,PMBM1,PMBM2,PMBP1,PMBP2,PMBP3 -C -C Cluster common used by soft event routines - PARAMETER (NMXCL=500) - COMMON/HWCLUS/PPCL(5,NMXCL),IDCL(NMXCL),NCL -C -C Parameters for resonant graviton production - DOUBLE PRECISION GRVLAM,EMGRV,GAMGRV - COMMON/HWGRAV/GRVLAM,EMGRV,GAMGRV -C -C Other new parameters for version 6.2 - DOUBLE PRECISION VIPWID,DXRCYL,DXZMAX,DXRSPH - LOGICAL WZRFR,FIX4JT - INTEGER IMSSM,IHIGGS,PARITY,LRSUSY - COMMON/HW6202/VIPWID(3),DXRCYL,DXZMAX,DXRSPH,WZRFR,FIX4JT, - & IMSSM,IHIGGS,PARITY,LRSUSY -C -C New parameters for version 6.203 - DOUBLE PRECISION ABWGT,ABWSUM,AVABW - INTEGER NNEGWT,NNEGEV - LOGICAL NEGWTS - COMMON/HW6203/ABWGT,ABWSUM,AVABW,NNEGWT,NNEGEV,NEGWTS -C -C New parameters for version 6.3 - INTEGER IMAXCH,IMAXOP - PARAMETER (IMAXCH=20,IMAXOP=40) - DOUBLE PRECISION MJJMIN,CHNPRB(IMAXCH) - INTEGER IOPSTP,IOPSH - LOGICAL OPTM,CHON(IMAXCH) - COMMON/HW6300/MJJMIN,CHNPRB,IOPSTP,IOPSH,OPTM,CHON -C New PDF's for version 6.3 - INTEGER NXMRS,NQMRS,NPMRS - PARAMETER(NXMRS=49,NQMRS=37,NPMRS=8) - DOUBLE PRECISION FMRS(3,NPMRS,NXMRS,NQMRS+1) - COMMON /HWPMRS/FMRS -C Circe interface for version 6.3 - INTEGER CIRCOP,CIRCAC,CIRCVR,CIRCRV,CIRCCH - COMMON /HWCIRC/CIRCOP,CIRCAC,CIRCVR,CIRCRV,CIRCCH -C New parameters and commons for spin correlations -C--constants for the arrays - INTEGER NMXSPN,NCFMAX - PARAMETER(NMXSPN=50,NCFMAX=3) - INTEGER NMODE2,NMODE3,NDIAGR,NMODEB,NMODE4 - PARAMETER(NMODE2=500,NMODE3=500,NDIAGR=8,NMODEB=50,NMODE4=4) -C--common block for X --> X gauge boson - DOUBLE PRECISION ABMODE(2,NMODEB),BBMODE(2,12,NMODEB), - & PBMODE(12,NMODEB),WTBMAX(12,NMODEB) - INTEGER IDBPRT(NMODEB),IBMODE(NMODEB),IBDRTP(NMODEB),NBMODE - COMMON /HWDSPB/ABMODE,BBMODE,PBMODE,WTBMAX,IDBPRT,IBDRTP,IBMODE, - & NBMODE -C--common block for two body decays - DOUBLE PRECISION A2MODE(2,NMODE2),P2MODE(NMODE2),WT2MAX(NMODE2) - INTEGER ID2PRT(NMODE2),I2DRTP(NMODE2),N2MODE - COMMON /HWDSP2/A2MODE,P2MODE,WT2MAX,ID2PRT,I2DRTP,N2MODE -C--common block for three body decays - DOUBLE PRECISION A3MODE(2,NDIAGR,NMODE3),B3MODE(2,NDIAGR,NMODE3), - & P3MODE(NMODE3),WT3MAX(NMODE3),SPN3CF(NCFMAX,NCFMAX,NMODE3) - INTEGER ID3PRT(NMODE3),I3MODE(NDIAGR,NMODE3), - & I3DRTP(NDIAGR,NMODE3),N3MODE,NDI3BY(NMODE3),N3NCFL(NMODE3), - & I3DRCF(NDIAGR,NMODE3) - COMMON /HWDSP3/A3MODE,B3MODE,P3MODE,WT3MAX,SPN3CF,ID3PRT,I3MODE, - & I3DRTP,N3MODE,NDI3BY,N3NCFL,I3DRCF -C--common block for four body decays - DOUBLE PRECISION A4MODE(2,12,NMODE4),B4MODE(2,12,NMODE4), - & P4MODE(12,12,NMODE4),WT4MAX(12,12,NMODE4) - INTEGER ID4PRT(NMODE4),I4MODE(2,NMODE4),N4MODE - COMMON /HWDSP4/A4MODE,B4MODE,P4MODE,WT4MAX,ID4PRT,I4MODE,N4MODE -C--common block for spin correlations in event - INTEGER NDECSY,NSEARCH,LRDEC,LWDEC - LOGICAL SYSPIN,THREEB,FOURB - CHARACTER *6 TAUDEC - COMMON /HWDSPN/NDECSY,NSEARCH,LRDEC,LWDEC,SYSPIN,THREEB, - & FOURB,TAUDEC - - INTEGER IDSPN(NMXSPN),JMOSPN(NMXSPN),JDASPN(2,NMXSPN),NSPN, - & ISNHEP(NMXHEP),NSNTRY,NCFL(NMXSPN),SPCOPT - DOUBLE COMPLEX MESPN(2,2,2,2,NCFMAX,NMXSPN),RHOSPN(2,2,NMXSPN) - DOUBLE PRECISION SPNCFC(NCFMAX,NCFMAX,NMXSPN) - LOGICAL DECSPN(NMXSPN) - COMMON /HWSPIN/MESPN,RHOSPN,SPNCFC,IDSPN,JMOSPN,JDASPN, - & NSPN,ISNHEP,NSNTRY,DECSPN,NCFL,SPCOPT - INTEGER JAK1,JAK2,ITDKRC,IFPHOT - COMMON /HWSTAU/ JAK1,JAK2,ITDKRC,IFPHOT -C -C--common block for Les Houches interface to store information we need -C - INTEGER MAXHRP - PARAMETER (MAXHRP=100) - DOUBLE PRECISION LHWGT(MAXHRP),LHWGTS(MAXHRP),LHMXSM, - & LHXSCT(MAXHRP),LHXERR(MAXHRP),LHXMAX(MAXHRP) - INTEGER LHIWGT(MAXHRP),ITYPLH,LHNEVT(MAXHRP) - LOGICAL LHSOFT,LHGLSF - COMMON /HWGUPR/LHWGT,LHWGTS,LHXSCT,LHXERR,LHXMAX,LHMXSM,LHIWGT, - & LHNEVT,ITYPLH,LHSOFT,LHGLSF -C -C--common block for HERWIG6.5 -C - LOGICAL PRESPL - COMMON /HW6500/ PRESPL diff --git a/Fpmc/External/excl_aaaa/CMakeLists.txt b/Fpmc/External/excl_aaaa/CMakeLists.txt index 71f67cb..f3bd07f 100644 --- a/Fpmc/External/excl_aaaa/CMakeLists.txt +++ b/Fpmc/External/excl_aaaa/CMakeLists.txt @@ -1,4 +1,5 @@ file(GLOB_RECURSE excl_aaaa_sources *.cpp) add_library(excl_aaaa OBJECT excl_aaaa_wraper.cpp ${excl_aaaa_sources}) + include_directories(commons) set_target_properties(excl_aaaa PROPERTIES LINKER_LANGUAGE CXX POSITION_INDEPENDENT_CODE ON) diff --git a/Fpmc/External/ntuple.f b/Fpmc/External/ntuple.f new file mode 100644 index 0000000..2b94bb8 --- /dev/null +++ b/Fpmc/External/ntuple.f @@ -0,0 +1,91 @@ +c//////////////////////////////////////////////// +c// @purpose: Simple ntuple +c// +c// @date: 20/02/2011 +c// @author: O. Kepka (kepkao@fzu.cz) +c//////////////////////////////////////////////// + + +c______________________________________________________________________________ +c___ create ntuple + + SUBROUTINE NTINIT +c______________________________________________________________________________ + + implicit none + INCLUDE '../Examples/ffcard.inc' + REAL hmemor + COMMON/pawc/hmemor(1000000) + +***** common block - generator level particles + integer ngenmax + parameter(ngenmax=1000) + integer ngen + real px(ngenmax),py(ngenmax),pz(ngenmax) + real e(ngenmax),m(ngenmax) + integer id(ngenmax) + integer ii(ngenmax) + integer ist(ngenmax) + + common /gener/ngen, + & px,py,pz,e,m,id,ii,ist + +***** common block - bjorken-x of the parton from pomeron +c common /remnant/xg1b,xg2b + + +c***** jets + integer njetmax + parameter(njetmax=30) + integer njets + real + & pxjet(njetmax),pyjet(njetmax), + & pzjet(njetmax),ejet(njetmax) + common/jets/njets, + & pxjet,pyjet, + & pzjet,ejet + + + INTEGER MODE, ICYCLE, ISTAT + + +* + print '(A)',' hropen ...' + + + call hbnt(777,'ntuple',' ') + + print *,'Creating the ntuple 777' + + +c call hbname(777,'remnant',xg1b,'xg1,'// +c & 'xg2') + + + + call hbname(777,'gener',ngen,'ngen[0,1000]:I,'// + & 'px(ngen),py(ngen),pz(ngen),e(ngen),rm(ngen),id(ngen),'// + & 'ii(ngen), ist(ngen)') + + if(UOUTPUT.ge.2) + & call hbname(777,'jets',njets,'njets[0,30]:I,'// + & 'pxjet(njets),pyjet(njets),'// + & 'pzjet(njets),pejet(njets)') + + + + END + +c______________________________________________________________________________ + + SUBROUTINE NTEND +c______________________________________________________________________________ + +*---Close the ntuple: + INTEGER idbg,igener,irad,ifrad + COMMON /CONST1/ idbg,igener,irad,ifrad + + call hrout(777,icycle,' ') + + END + diff --git a/Fpmc/External/ntuple_ww.f b/Fpmc/External/ntuple_ww.f new file mode 100644 index 0000000..028019c --- /dev/null +++ b/Fpmc/External/ntuple_ww.f @@ -0,0 +1,334 @@ +c rutine to initialize ntuples for detector simulation + + SUBROUTINE NTINIT + + + REAL hmemor + COMMON/pawc/hmemor(8000000) + +c______________________________________________________________________________ +c common block for ntuples +c exactly the same structure as simul_interf_cms.f + + +***** gener common block + integer ngenmax + parameter(ngenmax=1000) + integer ngen + real px(ngenmax),py(ngenmax),pz(ngenmax) + real e(ngenmax),rm(ngenmax) + integer id(ngenmax) + + common /gener/ngen, + & px,py,pz,e,rm,id + +***** remnant common block + common /remnant/xg1b,xg2b + + +***** renostruction common block +C CHR 25/05/99 modif pour ntuple reconstruit + integer nrecmax + parameter(nrecmax=100) + integer nrec + real typrec(nrecmax),pxrec(nrecmax),pyrec(nrecmax) + real pzrec(nrecmax) + real erec(nrecmax),qrec(nrecmax),eemrec(nrecmax), + & ehadrec(nrecmax) + real etrarec(nrecmax),widrec(nrecmax),var6(nrecmax) + real ctagrec(nrecmax),btag1rec(nrecmax),btag2rec(nrecmax) + integer ntrarec(nrecmax) + real sum15ec(nrecmax),sum15hc(nrecmax), sum40(nrecmax) + real sum40ec(nrecmax) + + common /recons1/ nrec, + & typrec,pxrec,pyrec,pzrec,erec,qrec,eemrec,ehadrec, + & etrarec,widrec,var6, + & ctagrec,btag1rec,btag2rec, + & ntrarec, + & sum15ec,sum15hc, sum40, sum40ec + + +***** global common block + real circul,transen,transmass,ptmisstx,ptmissty + real ptmissx,ptmissy,pznu1,pznu2 + + common /global/ circul,transen,transmass,ptmisstx,ptmissty, + . ptmissx,ptmissy,pznu1,pznu2 + + +******** track common block +* integer ntramax +* parameter(ntramax=150) +* integer ntra +* real pxtra(ntramax),pytra(ntramax),pztra(ntramax) +* real qtra(ntramax) + +* common /track/ ntra,pxtra,pytra,pztra,qtra + + +***** cluster common block +c integer nclmax +c parameter(nclmax=500) +c integer nclu +c real pxrclu(nclmax),pyrclu(nclmax),pzrclu(nclmax) +c real eclu(nclmax) +c real eemclu(nclmax) +c real ehadclu(nclmax),widtclu(nclmax),rmisset,phimet +c integer multclu(nclmax) + +c common /cluster/ nclu,pxrclu,pyrclu,pzrclu,eclu,eemclu, +c & ehadclu,widtclu,multclu,rmisset,phimet + + +***** trigger common block +c integer ntrimax +c parameter(ntrimax=50) +c integer ntri +c real typtri(ntrimax),etatri(ntrimax),phitri(ntrimax) +c real ettri(ntrimax) +c real clutri(ntrimax) +c real pttri(ntrimax),tratri(ntrimax) + +c common /trigger/ ntri,typtri,etatri,phitri,ettri,clutri,pttri, +c & tratri + +***** cell common block +c integer ncellmax +c parameter(ncellmax=5000) +c integer ncell +c real pxcell(ncellmax),pycell(ncellmax) +c real etsheta(ncellmax) +c real ecell(ncellmax),etcell(ncellmax) +c real etacell(ncellmax) +c real phicell(ncellmax),rmrkcell(ncellmax) + +c common /cells/ ncell, pxcell, pycell, etsheta, ecell, etcell, +c + etacell,phicell,rmrkcell + + +c OK:10/06: following is not used also. +***** jets : first jets from LUCELL, then LUCLUS + integer njetmax + parameter(njetmax=30) + integer njetscel + real + & pxjet(njetmax),pyjet(njetmax), + & pzjet(njetmax),pejet(njetmax), + & pmjet(njetmax) + common/celljets/ + & njetscel, + & pxjet,pyjet, + & pzjet,pejet, + & pmjet + + + integer nclumax + parameter(nclumax=30) + integer njetsclu + real + & pxclu(nclumax),pyclu(nclumax), + & pzclu(nclumax),peclu(nclumax), + & pmclu(nclumax) + common/clujets/ + & njetsclu, + & pxclu,pyclu, + & pzclu,peclu, + & pmclu + + INTEGER i +*...Max. number of preclusters/ clusters/jets + INTEGER NCLUS + PARAMETER (NCLUS = 50) +*KEEP,PTMIS. +*...Et missing true (neutrinos), Et miss. measured, reconstructed Z-compon. of +*...neutrino from W -> l nu decay (if any). + REAL PTMT, PTM, PZNU + COMMON /PTMIS/ PTMT(2), PTM(2), PZNU(2) +*KEEP,TRANCH. +*...Circularity, transverse energy and mass values calculated in sub. TRAPAR + REAL CIRC, TREN, TRMA + COMMON /TRANCH/ CIRC, TREN, TRMA +*KEEP,LEPTS. +*...Desired leptons found by user (PROLEP, then after CISOL check and, +*...possibly, after GLOBJF), +*...PL - 4-VECTORS, KL - code : +-1(2) for e(mu)-+ ; +*...INDLEP - number of each lepton in /PLIST/, LHIT - in /HITPAR/, +*...MRL - origin (code of parent, =999 for SUSY particles), +*...LISOL - isolation mark (0 - isolated, 1 - non-isolated in tracker, +*... 2 - non-isolated in calo, 3 - non-isolated in both tracker and calo, +*... 4 - isolated in jet,), NLE - number of leptons. + INTEGER NLMAX, KL, INDLEP, LHIT, MRL, LISOL, NLE + PARAMETER (NLMAX = 100) + REAL PL + COMMON /LEPTS/ PL(4,NLMAX), KL(NLMAX), INDLEP(NLMAX), + & LHIT(NLMAX), MRL(NLMAX), LISOL(NLMAX), NLE +*KEEP,GAMMAS. +*...Gammas found (PROGAM, then CISOL and possibly GLOBJF), +*...PG - 4-VECTORS, +*...INDGAM - number of each lepton in /PLIST/, +*...IGHIT - number in /HITPAR/, +*...LISGAM - isolation mark : 0 - isolated, +*... 2 - non-isolated in calo, 4 - isolated in jet; +*...NGAM - number of gammas. + INTEGER NGMAX, INDGAM, IGHIT, LISGAM, NGAM + REAL PG + PARAMETER (NGMAX = 100) + COMMON /GAMMAS/ PG(4,NGMAX), INDGAM(NGMAX), + & IGHIT(NGMAX), LISGAM(NGMAX), NGAM +*KEEP,JETOUT. +*...Jets found by subroutine GLOBJF + INTEGER NJG + REAL PJG + COMMON /JETOUT/ PJG(4,NCLUS), NJG +*KEEP,JEMARK. +*...Jet marks filled by sub. QJMAT : 0 - quark/gluon not assigned, +*... > 0 - PYTHIA`s code of the quark/gluon (without sign). +*...CALL TAUTAG can add 1000 to this mark of jet if recognized as 1-prong tau. +*...BTAG contains a likelihood of jet b-tagging (if IBTAG =/= 0) + INTEGER MRKJET + REAL BTAG + COMMON /JEMARK/ MRKJET(NCLUS), BTAG(NCLUS) +*KEEP,DNAMES. +*...Names and LUNs of the I/O files. + INTEGER LDAT1, LDAT2, LDAT3, LDAT4, LDAT5, LDAT6, + & LRNDM1, LRNDM2, IDNTIN, IDNTOUT + CHARACTER*80 DATNA1, DATNA2, DATNA3, DATNA4, DATNA5, DATNA6, + & RNDMN1, RNDMN2, HISTNA, HISTNB + COMMON /DNAMES/ LDAT1, LDAT2, LDAT3, LDAT4, LDAT5, LDAT6, + & LRNDM1, LRNDM2, IDNTIN, IDNTOUT, + & DATNA1, DATNA2, DATNA3, DATNA4, DATNA5, DATNA6, + & RNDMN1, RNDMN2, HISTNA, HISTNB +*KEND. +*............ + + INTEGER NLE_CWN, KL_CWN, MRL_CWN, LISOL_CWN + REAL PL_CWN + COMMON /LEP_CWN/ NLE_CWN, PL_CWN(4,NLMAX), KL_CWN(NLMAX), + & MRL_CWN(NLMAX), LISOL_CWN(NLMAX) + + INTEGER NGAM_CWN, LISGAM_CWN + REAL PG_CWN + COMMON /GAM_CWN/ NGAM_CWN, PG_CWN(4,NGMAX), LISGAM_CWN(NGMAX) + + INTEGER NJG_CWN, MRKJET_CWN + REAL PJG_CWN, BTAG_CWN + COMMON /JET_CWN/ NJG_CWN, PJG_CWN(4,NCLUS), MRKJET_CWN(NCLUS), + & BTAG_CWN(NCLUS) + + + INTEGER MODE, ICYCLE, ISTAT + + + +c______________________________________________________________________________ +c create ntuple +* + print '(A)',' hropen ...' + +* call hropen(31,'pythia','higgs.ntp','n',1024,istat) + +c IF (ISTAT.NE. 0) THEN +c PRINT * ,'*** ERROR ',ISTAT,' DURING OPEN OF FILE 31' +c CLOSE(31) +c RETURN +c ENDIF + +c print '(A)',' hbnt ...' + +c call hmdir('//pythia','s') +c call hcdir('//pythia',' ') + call hbnt(777,'ntuple',' ') + + print *,'on vient de creer le ntuple 777' + + + +c call hbname(777,'LUCELL',njetscel,'njetscel[0,30]:U,'// +c &'pxjet(njetscel),pyjet(njetscel),'// +c &'pzjet(njetscel),pejet(njetscel),pmjet(njetscel)') + +c call hbname(777,'LUCLUS',njetsclu,'njetsclu[0,30]:U,'// +c &'pxclu(njetsclu),pyclu(njetsclu),'// +c &'pzclu(njetsclu),peclu(njetsclu),pmclu(njetsclu)') + + call hbname(777,'recons1',nrec,'nrec[0,100]:I,'// + & 'typrec(nrec),'// + & 'pxrec(nrec),pyrec(nrec),pzrec(nrec),'// + & 'erec(nrec),'// + & 'qrec(nrec),eemrec(nrec),ehadrec(nrec),'// + & 'etrarec(nrec),widrec(nrec),var6(nrec),'// + & 'ctagrec(nrec),btag1rec(nrec),btag2rec(nrec),'// + & 'ntrarec(nrec),'// + & 'sum15ec(nrec),sum15hc(nrec), sum40(nrec), sum40ec(nrec)') + +c call hbname(777,'track',ntra,'ntra[0,100]:I,'// +c & 'pxtra(ntra),pytra(ntra),pztra(ntra),qtra(ntra)') + + + call hbname(777,'global',circul,'circul,transen,transmass,'// + & 'ptmisstx,ptmissty,'// + & 'ptmissx,ptmissy,pznu1,pznu2') + + +* do not put cells in root tuple +c call hbname(777,'cells',ncell,'ncell[0,5000]:I,'// +c & 'pxcell(ncell),'// +c & 'pycell(ncell), etsheta(ncell), ecell(ncell), etcell(ncell),'// +c & 'etacell(ncell),phicell(ncell),rmrkcell(ncell)') + + + + + call hbname(777,'remnant',xg1b,'xg1,'// + & 'xg2') + + +c call hbname(777,'cluster',nclu,'nclu[0,500]:I,'// +c & 'pxrclu(nclu),pyrclu(nclu),pzrclu(nclu),eclu(nclu),'// +c & 'eemclu(nclu),ehadclu(nclu),'// +c & 'widtclu(nclu),multclu(nclu),rmisset,phimet') + +c call hbname(777,'trigger',ntri,'ntri[0,50]:I,'// +c & 'typtri(ntri),etatri(ntri),phitri(ntri),ettri(ntri),'// +c & 'clutri(ntri),pttri(ntri),tratri(ntri)') + + call hbname(777,'gener',ngen,'ngen[0,1000]:I,'// + & 'px(ngen),py(ngen),pz(ngen),e(ngen),rm(ngen),id(ngen)') + +c CALL HBNAME(777,'TRANCH',CIRC,'CIRC,'// +c & 'TREN,TRMA') +c CALL HBNAME(777,'PTMIS', PTMT,'PTMT(2),'// +c & 'PTM(2),PZNU(2)') +c CALL HBNAME(777,'LEPTONS',NLE_CWN,'NLE_CWN[0,100],'// +c & 'PL_CWN(4,NLE_CWN),KL_CWN(NLE_CWN)[-2,2],'// +c & 'MRL_CWN(NLE_CWN),LISOL_CWN(NLE_CWN)[0,4]') +c CALL HBNAME(777,'GAMMAS',NGAM_CWN,'NGAM_CWN[0,100],'// +c & 'PG_CWN(4,NGAM_CWN),'// +c & 'LISGAM_CWN(NGAM_CWN)[0,4]') +c CALL HBNAME(777,'JETS',NJG_CWN,'NJG_CWN[0,50],'// +c & 'PJG_CWN(4,NJG_CWN),'// +c & 'MRKJET_NJG(NJG_CWN)[0,2000],'// +c & 'BTAG_NJG(NJG_CWN)') + + return + end + + + + SUBROUTINE NTEND + +*KEND. +* close the ntuple: + +*KEEP,CONST. + INTEGER idbg,igener,irad,ifrad + COMMON /CONST1/ idbg,igener,irad,ifrad + + + + call hrout(777,icycle,' ') + + + END + diff --git a/Fpmc/External/reco_ok/HERWIG65.INC b/Fpmc/External/reco_ok/HERWIG65.INC new file mode 100644 index 0000000..7a080d5 --- /dev/null +++ b/Fpmc/External/reco_ok/HERWIG65.INC @@ -0,0 +1,2 @@ +* INCLUDE 'herwig6510.inc' + INCLUDE '../../Herwig/herwig6500.inc' diff --git a/Fpmc/External/reco_ok/pxcone_mod_new.f b/Fpmc/External/reco_ok/pxcone_mod_new.f new file mode 100644 index 0000000..466aae2 --- /dev/null +++ b/Fpmc/External/reco_ok/pxcone_mod_new.f @@ -0,0 +1,988 @@ + SUBROUTINE PXCONE(MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM, + + MXJET,NJET,PJET,IPASS,IJMUL,IERR) +*.********************************************************* +*. ------ +*. PXCONE +*. ------ +* +* MW 07/02/2003: huge change - implement mode=3 (E-scheme) +* +*. +*.********** Pre Release Version 26.2.93 +*. +*. Driver for the Cone Jet finding algorithm of L.A. del Pozo. +*. Based on algorithm from D.E. Soper. +*. Finds jets inside cone of half angle CONER with energy > EPSLON. +*. Jets which receive more than a fraction OVLIM of their energy from +*. overlaps with other jets are excluded. +*. Output jets are ordered in energy. +*. If MODE.EQ.2 momenta are stored as (eta,phi,,pt) +*. Usage : +*. +*. INTEGER ITKDM,MXTRK +*. PARAMETER (ITKDM=4.or.more,MXTRK=1.or.more) +*. INTEGER MXJET, MXTRAK, MXPROT +*. PARAMETER (MXJET=10,MXTRAK=900,MXPROT=500) +*. INTEGER IPASS (MXTRAK),IJMUL (MXJET) +*. INTEGER NTRAK,NJET,IERR,MODE +*. DOUBLE PRECISION PTRAK (ITKDM,MXTRK),PJET (5,MXJET) +*. DOUBLE PRECISION CONER, EPSLON, OVLIM +*. NTRAK = 1.to.MXTRAK +*. CONER = ... +*. EPSLON = ... +*. OVLIM = ... +*. CALL PXCONE (MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,MXJET, +*. + NJET,PJET,IPASS,IJMUL,IERR) +*. +*. INPUT : MODE 1=>e+e-, 2=>hadron-hadron 3 RunII: Escheme +*. INPUT : NTRAK Number of particles +*. INPUT : ITKDM First dimension of PTRAK array +*. INPUT : PTRAK Array of particle 4-momenta (Px,Py,Pz,E) +*. INPUT : CONER Cone size (half angle) in radians +*. INPUT : EPSLON Minimum Jet energy (GeV) +*. INPUT : OVLIM Maximum fraction of overlap energy in a jet +*. INPUT : MXJET Maximum possible number of jets +*. OUTPUT : NJET Number of jets found +*. OUTPUT : PJET 5-vectors of jets +*. OUTPUT : IPASS(k) Particle k belongs to jet number IPASS(k) +*. IPASS = -1 if not assosciated to a jet +*. OUTPUT : IJMUL(i) Jet i contains IJMUL(i) particles +*. OUTPUT : IERR = 0 if all is OK ; = -1 otherwise +*. +*. CALLS : PXSEAR, PXSAME, PXNEW, PXTRY, PXORD, PXUVEC, PXOLAP +*. CALLED : User +*. +*. AUTHOR : L.A. del Pozo +*. CREATED : 26-Feb-93 +*. LAST MOD : 2-Mar-93 +*. +*. Modification Log. +*. 2-Jan-97: M Wobisch - fix bug concerning COS2R in eta phi mode +*. 4-Apr-93: M H Seymour - Change 2d arrays to 1d in PXTRY & PXNEW +*. 2-Apr-93: M H Seymour - Major changes to add boost-invariant mode +*. 1-Apr-93: M H Seymour - Increase all array sizes +*. 30-Mar-93: M H Seymour - Change all REAL variables to DOUBLE PRECISION +*. 30-Mar-93: M H Seymour - Change OVLIM into an input parameter +*. 2-Mar-93: L A del Pozo - Fix Bugs in PXOLAP +*. 1-Mar-93: L A del Pozo - Remove Cern library routine calls +*. 1-Mar-93: L A del Pozo - Add Print out of welcome and R and Epsilon +*. +*.********************************************************* +C+SEQ,DECLARE. +*** External Arrays + INTEGER ITKDM,MXJET,NTRAK,NJET,IERR,MODE + INTEGER IPASS (*),IJMUL (MXJET) + DOUBLE PRECISION PTRAK (ITKDM,*),PJET (5,*), CONER, EPSLON, OVLIM +*** Internal Arrays + INTEGER MXPROT, MXTRAK + PARAMETER (MXPROT=500, MXTRAK=900) + DOUBLE PRECISION PP(4,MXTRAK), PU(3,MXTRAK), PJ(4,MXPROT) + LOGICAL JETLIS(MXPROT,MXTRAK) +*** Used in the routine. + DOUBLE PRECISION COSR,COS2R, VSEED(3), VEC1(3), VEC2(3),PTSQ,PPSQ, + + COSVAL,PXMDPI +cMWobisch + DOUBLE PRECISION RSEP, MWDIST +cMWobisch + LOGICAL UNSTBL + INTEGER I,J,N,MU,N1,N2, ITERR + INTEGER NCALL, NPRINT + DOUBLE PRECISION ROLD, EPSOLD, OVOLD, HMW(3,MXTRAK) + COMMON /HMWCOMM/ HMW + SAVE NCALL,NPRINT,ROLD, EPSOLD, OVOLD + DATA NCALL,NPRINT /0,0/ + DATA ROLD,EPSOLD,OVOLD/0.,0.,0./ + +cMWobisch +c*************************************** + RSEP = 2D0 +c RSEP = 1D0 ! only meaningful in NLO calculations!!!! +c*************************************** +cMWobisch + IERR=0 +* +*** INITIALIZE + IF(NCALL.LE.0) THEN + ROLD = 0. + EPSOLD = 0. + OVOLD = 0. + ENDIF + NCALL = NCALL + 1 +* +*** Print welcome and Jetfinder parameters + IF((CONER.NE.ROLD .OR. EPSLON.NE.EPSOLD .OR. OVLIM.NE.OVOLD) + + .AND. NPRINT.LE.10) THEN + WRITE (6,*) + WRITE (6,*) ' *********** PXCONE: Cone Jet-finder ***********' + WRITE (6,*) ' Written by Luis Del Pozo of OPAL' + WRITE (6,*) ' Modified for eta-phi by Mike Seymour' + WRITE (6,*) ' E-scheme for eta-phi by Markus Wobisch' + WRITE (6,*) ' mode=',mode + WRITE(6,1000)' Cone Size R = ',CONER,' Radians' + WRITE(6,1001)' Min Jet energy Epsilon = ',EPSLON,' GeV' + WRITE(6,1002)' Overlap fraction parameter = ',OVLIM + WRITE (6,*) ' ***********************************************' +cMWobisch + IF (RSEP .lt. 1.999) THEN + WRITE(6,*) ' ' + WRITE (6,*) ' ******************************************' + WRITE (6,*) ' ******************************************' + WRITE(6,*) ' M Wobisch: private change !!!!!!!!!!!! ' + WRITE(6,*) ' Rsep is set to ',RSEP + WRITE(6,*) ' this is ONLY meaningful in a NLO calculation' + WRITE(6,*) ' ------------------------ ' + WRITE(6,*) ' please check what you''re doing!!' + WRITE(6,*) ' or ask: wobisch@fnal.gov --' + WRITE (6,*) ' ******************************************' + WRITE (6,*) ' ******************************************' + WRITE (6,*) ' ******************************************' + WRITE(6,*) ' ' + WRITE(6,*) ' ' + ENDIF +cMWobisch + + WRITE (6,*) +1000 FORMAT(A18,F5.2,A10) +1001 FORMAT(A29,F5.2,A5) +1002 FORMAT(A33,F5.2) + NPRINT = NPRINT + 1 + ROLD=CONER + EPSOLD=EPSLON + OVOLD=OVLIM + ENDIF +* +*** Copy calling array PTRAK to internal array PP(4,NTRAK) +* + IF (NTRAK .GT. MXTRAK) THEN + WRITE (6,*) ' PXCONE: Ntrak too large' + IERR=-1 + RETURN + ENDIF + IF (MODE.NE.2) THEN + DO 100 I=1, NTRAK + DO 101 J=1,4 + PP(J,I)=PTRAK(J,I) +101 CONTINUE + +cMW - create 'help' array with particle eta/phi values - maybe later pT + if (mode.eq.3) then + PTSQ=PTRAK(1,I)**2+PTRAK(2,I)**2 + PPSQ=(SQRT(PTSQ+PTRAK(3,I)**2)+ABS(PTRAK(3,I)))**2 + IF (PTSQ.LE.4.25d-18*PPSQ) THEN + HMW(1,i)=20d0 + ELSE + HMW(1,i)=0.5d0*LOG(PPSQ/PTSQ) + ENDIF + HMW(1,i)=SIGN(HMW(1,i),PTRAK(3,I)) + IF (PTSQ.EQ.0d0) THEN + HMW(2,i)=0d0 + ELSE + HMW(2,i)=ATAN2(PTRAK(2,I),PTRAK(1,I)) + ENDIF + endif + +100 CONTINUE + ELSE +*** Converting to eta,phi,pt if necessary + DO 104 I=1,NTRAK + PTSQ=PTRAK(1,I)**2+PTRAK(2,I)**2 + PPSQ=(SQRT(PTSQ+PTRAK(3,I)**2)+ABS(PTRAK(3,I)))**2 + IF (PTSQ.LE.4.25E-18*PPSQ) THEN + PP(1,I)=20 + ELSE + PP(1,I)=0.5*LOG(PPSQ/PTSQ) + ENDIF + PP(1,I)=SIGN(PP(1,I),PTRAK(3,I)) + IF (PTSQ.EQ.0) THEN + PP(2,I)=0 + ELSE + PP(2,I)=ATAN2(PTRAK(2,I),PTRAK(1,I)) + ENDIF + PP(3,I)=0 + PP(4,I)=SQRT(PTSQ) + PU(1,I)=PP(1,I) + PU(2,I)=PP(2,I) + PU(3,I)=PP(3,I) +104 CONTINUE + ENDIF +* +*** Zero output variables +* + NJET=0 + DO 102 I = 1, NTRAK + DO 103 J = 1, MXPROT + JETLIS(J,I) = .FALSE. +103 CONTINUE +102 CONTINUE + CALL PXZERV(4*MXPROT,PJ) + CALL PXZERI(MXJET,IJMUL) +* + IF (MODE.EQ.3) THEN ! MW -> new mode + COSR = 1-CONER**2 + COS2R = 1-(RSEP*CONER)**2 + ELSEIF (MODE.NE.2) THEN + COSR = COS(CONER) +cMW COS2R = COS(CONER) ! probably a bug?? + COS2R = COS(2*CONER) + ELSE +*** Purely for convenience, work in terms of 1-R**2 + COSR = 1-CONER**2 +cMW -- select Rsep: 1-(Rsep*CONER)**2 + COS2R = 1-(RSEP*CONER)**2 +cORIGINAL COS2R = 1-(2*CONER)**2 + ENDIF + UNSTBL = .FALSE. + IF (MODE.NE.2) THEN + CALL PXUVEC(NTRAK,PP,PU,IERR) + IF (IERR .NE. 0) RETURN + ENDIF +*** Look for jets using particle diretions as seed axes +* + DO 110 N = 1,NTRAK + DO 120 MU = 1,3 + VSEED(MU) = PU(MU,N) +120 CONTINUE + CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED, + & NJET,JETLIS,PJ,UNSTBL,IERR) +c IF (IERR .NE. 0) RETURN + IF (IERR .NE. 0) then + write(*,*) '1: ',vseed + RETURN + endif +110 CONTINUE + +cMW - for Rsep=1 goto 145 +c GOTO 145 + +*** Now look between all pairs of jets as seed axes. + DO 140 N1 = 1,NJET-1 + VEC1(1)=PJ(1,N1) + VEC1(2)=PJ(2,N1) + VEC1(3)=PJ(3,N1) + IF (MODE.NE.2) CALL PXNORV(3,VEC1,VEC1,ITERR) + DO 150 N2 = N1+1,NJET + VEC2(1)=PJ(1,N2) + VEC2(2)=PJ(2,N2) + VEC2(3)=PJ(3,N2) + IF (MODE.NE.2) CALL PXNORV(3,VEC2,VEC2,ITERR) + CALL PXADDV(3,VEC1,VEC2,VSEED,ITERR) + IF (MODE.NE.2) THEN + CALL PXNORV(3,VSEED,VSEED,ITERR) + ELSE + VSEED(1)=VSEED(1)/2 + VSEED(2)=VSEED(2)/2 + ENDIF +C---ONLY BOTHER IF THEY ARE BETWEEN 1 AND 2 CONE RADII APART + IF (MODE.eq.3) THEN + COSVAL = 1-MWDIST(VEC1,VEC2) + ELSEIF (MODE.NE.2) THEN + COSVAL=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3) + ELSE + IF (ABS(VEC1(1)).GE.20.OR.ABS(VEC2(1)).GE.20) THEN + COSVAL=-1000 + ELSE + COSVAL=1- + + ((VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2) + ENDIF + ENDIF + + IF (COSVAL.LE.COSR.AND.COSVAL.GE.COS2R) + + CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET, + + JETLIS,PJ,UNSTBL,IERR) +c CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET, +c + JETLIS,PJ,UNSTBL,IERR) +c IF (IERR .NE. 0) RETURN + IF (IERR .NE. 0) then + write(*,*) '2: ',vseed + RETURN + endif +150 CONTINUE +140 CONTINUE + IF (UNSTBL) THEN + IERR=-1 + WRITE (6,*) ' PXCONE: Too many iterations to find a proto-jet' + RETURN + ENDIF + + 145 CONTINUE +*** Now put the jet list into order by jet energy, eliminating jets +*** with energy less than EPSLON. + CALL PXORD(MODE,EPSLON,NJET,NTRAK,JETLIS,PJ) +* +*** Take care of jet overlaps + CALL PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM) +* +*** Order jets again as some have been eliminated, or lost energy. + CALL PXORD(MODE,EPSLON,NJET,NTRAK,JETLIS,PJ) +* +*** All done!, Copy output into output arrays + IF (NJET .GT. MXJET) THEN + WRITE (6,*) ' PXCONE: Found more than MXJET jets' + IERR=-1 + GOTO 99 + ENDIF + IF (MODE.NE.2) THEN + DO 300 I=1, NJET + DO 310 J=1,4 + PJET(J,I)=PJ(J,I) +310 CONTINUE +300 CONTINUE + ELSE + DO 315 I=1, NJET + PJET(1,I)=PJ(4,I)*COS(PJ(2,I)) + PJET(2,I)=PJ(4,I)*SIN(PJ(2,I)) + PJET(3,I)=PJ(4,I)*SINH(PJ(1,I)) + PJET(4,I)=PJ(4,I)*COSH(PJ(1,I)) + 315 CONTINUE + ENDIF + DO 320 I=1, NTRAK + IPASS(I)=-1 + DO 330 J=1, NJET + IF (JETLIS(J,I)) THEN + IJMUL(J)=IJMUL(J)+1 + IPASS(I)=J + ENDIF +330 CONTINUE +320 CONTINUE +99 RETURN + END +*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper +*-- Author : +C----------------------------------------------------------------------- + SUBROUTINE PXNORV(N,A,B,ITERR) + INTEGER I,N,ITERR + DOUBLE PRECISION A(N),B(N),C + C=0 + DO 10 I=1,N + C=C+A(I)**2 + 10 CONTINUE + IF (C.LE.0) RETURN + C=1/SQRT(C) + DO 20 I=1,N + B(I)=A(I)*C + 20 CONTINUE + END +*CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper +*CMZ : 1.06/00 15/03/94 12.17.46 by P. Schleper +*-- Author : +* +C+DECK,PXOLAP. + SUBROUTINE PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM) +* +*** Looks for particles assigned to more than 1 jet, and reassigns them +*** If more than a fraction OVLIM of a jet's energy is contained in +*** higher energy jets, that jet is neglected. +*** Particles assigned to the jet closest in angle (a la CDF, Snowmass). +C+SEQ,DECLARE. + INTEGER MXTRAK, MXPROT + PARAMETER (MXTRAK=900,MXPROT=500) + INTEGER NJET, NTRAK, MODE + LOGICAL JETLIS(MXPROT,MXTRAK) + DOUBLE PRECISION PJ(4,MXPROT),PP(4,MXTRAK),PXMDPI + INTEGER I,J,N,MU + LOGICAL OVELAP + DOUBLE PRECISION EOVER, OVERJET, MWDIST + DOUBLE PRECISION OVLIM + INTEGER ITERR, IJMIN, IJET(MXPROT), NJ + DOUBLE PRECISION VEC1(3), VEC2(3), COST, THET, THMIN + DATA IJMIN/0/ +* + IF (NJET.LE.1) RETURN +*** Look for jets with large overlaps with higher energy jets. + DO 100 I = 2,NJET +*** Find overlap energy between jets I and all higher energy jets. + EOVER = 0.0 + DO 110 N = 1,NTRAK + OVELAP = .FALSE. + DO 120 J= 1,I-1 + IF (JETLIS(I,N).AND.JETLIS(J,N)) THEN + OVELAP = .TRUE. + ENDIF +120 CONTINUE + IF (OVELAP) THEN + if (MODE.ne.3) then + EOVER = EOVER + PP(4,N) + else + EOVER = EOVER + sqrt(PP(1,N)**2+PP(2,N)**2) + endif + ENDIF +110 CONTINUE +*** Is the fraction of energy shared larger than OVLIM? + if (MODE.ne.3) then + OVERJET = PJ(4,I) + else + OVERJET = sqrt(PJ(1,I)**2+PJ(2,I)**2) + endif + IF (EOVER.GT.OVLIM*OVERJET) THEN +*** De-assign all particles from Jet I + DO 130 N = 1,NTRAK + JETLIS(I,N) = .FALSE. +130 CONTINUE + ENDIF +100 CONTINUE +*** Now there are no big overlaps, assign every particle in +*** more than 1 jet to the closet jet. +*** Any particles now in more than 1 jet are assigned to the CLOSET +*** jet (in angle). + DO 140 I=1,NTRAK + NJ=0 + DO 150 J=1, NJET + IF(JETLIS(J,I)) THEN + NJ=NJ+1 + IJET(NJ)=J + ENDIF +150 CONTINUE + IF (NJ .GT. 1) THEN +*** Particle in > 1 jet - calc angles... + VEC1(1)=PP(1,I) + VEC1(2)=PP(2,I) + VEC1(3)=PP(3,I) + THMIN=0. + DO 160 J=1,NJ + VEC2(1)=PJ(1,IJET(J)) + VEC2(2)=PJ(2,IJET(J)) + VEC2(3)=PJ(3,IJET(J)) + IF (MODE.EQ.3) THEN + THET = MWDIST(VEC1,VEC2) + ELSEIF (MODE.NE.2) THEN + CALL PXANG3(VEC1,VEC2,COST,THET,ITERR) + ELSE + THET=(VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2 + ENDIF + IF (J .EQ. 1) THEN + THMIN=THET + IJMIN=IJET(J) + ELSEIF (THET .LT. THMIN) THEN + THMIN=THET + IJMIN=IJET(J) + ENDIF +160 CONTINUE +*** Assign track to IJMIN + DO 170 J=1,NJET + JETLIS(J,I) = .FALSE. +170 CONTINUE + JETLIS(IJMIN,I)=.TRUE. + ENDIF +140 CONTINUE +*** Recompute PJ + DO 200 I = 1,NJET + DO 210 MU = 1,4 + PJ(MU,I) = 0.0 +210 CONTINUE + DO 220 N = 1,NTRAK + IF( JETLIS(I,N) ) THEN + IF (MODE.NE.2) THEN + DO 230 MU = 1,4 + PJ(MU,I) = PJ(MU,I) + PP(MU,N) +230 CONTINUE + ELSE + PJ(1,I)=PJ(1,I) + + + PP(4,N)/(PP(4,N)+PJ(4,I))*(PP(1,N)-PJ(1,I)) + PJ(2,I)=PJ(2,I) + + + PP(4,N)/(PP(4,N)+PJ(4,I))*PXMDPI(PP(2,N)-PJ(2,I)) + PJ(4,I)=PJ(4,I)+PP(4,N) + ENDIF + ENDIF +220 CONTINUE +200 CONTINUE + RETURN + END +*CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper +*CMZ : 1.06/00 14/03/94 15.37.45 by P. Schleper +*-- Author : +* +C+DECK,PXORD. + SUBROUTINE PXORD(MODE,EPSLON,NJET,NTRAK,JETLIS,PJ) +* +*** Routine to put jets into order and eliminate tose less than EPSLON +C+SEQ,DECLARE. + INTEGER MXTRAK,MXPROT, MODE + PARAMETER (MXTRAK=900,MXPROT=500) + INTEGER I, J, INDEX(MXPROT) + DOUBLE PRECISION PTEMP(4,MXPROT), ELIST(MXPROT) + INTEGER NJET,NTRAK + LOGICAL JETLIS(MXPROT,MXTRAK) + LOGICAL LOGTMP(MXPROT,MXTRAK) + DOUBLE PRECISION EPSLON,PJ(4,MXPROT), ECUT +*** Puts jets in order of energy: 1 = highest energy etc. +*** Then Eliminate jets with energy below EPSLON +* +*** Copy input arrays. + DO 100 I=1,NJET + DO 110 J=1,4 + PTEMP(J,I)=PJ(J,I) +110 CONTINUE + DO 120 J=1,NTRAK + LOGTMP(I,J)=JETLIS(I,J) +120 CONTINUE +100 CONTINUE + DO 150 I=1,NJET + if (MODE.ne.3) then + ELIST(I)=PJ(4,I) + else +cMW - new Escheme -> sort in pT + ELIST(I)=sqrt(PJ(1,I)**2+PJ(2,I)**2) + endif +150 CONTINUE +*** Sort the energies... + CALL PXSORV(NJET,ELIST,INDEX,'I') +*** Fill PJ and JETLIS according to sort ( sort is in ascending order!!) + DO 200 I=1, NJET + DO 210 J=1,4 + PJ(J,I)=PTEMP(J,INDEX(NJET+1-I)) +210 CONTINUE + DO 220 J=1,NTRAK + JETLIS(I,J)=LOGTMP(INDEX(NJET+1-I),J) +220 CONTINUE +200 CONTINUE +** Jets are now in order +*** Now eliminate jets with less than Epsilon energy + DO 300, I=1, NJET + if (MODE.eq.3) then + ECUT = sqrt(PJ(1,I)**2+PJ(2,I)**2) + else + ecut = PJ(4,I) + endif + IF (ECUT .LT. EPSLON) THEN + NJET=NJET-1 + PJ(4,I)=0. + ENDIF +300 CONTINUE + RETURN + END + +******************************************************************** +*CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper +*CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper +*-- Author : +C+DECK,PXSEAR. + SUBROUTINE PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET, + + JETLIS,PJ,UNSTBL,IERR) +* +C+SEQ,DECLARE. + INTEGER MXTRAK, MXPROT + PARAMETER (MXTRAK=900,MXPROT=500) + INTEGER NTRAK, IERR, MODE + DOUBLE PRECISION COSR,PU(3,MXTRAK),PP(4,MXTRAK),VSEED(3) + LOGICAL UNSTBL + LOGICAL JETLIS(MXPROT,MXTRAK) + INTEGER NJET + DOUBLE PRECISION PJ(4,MXPROT) +*** Using VSEED as a trial axis , look for a stable jet. +*** Check stable jets against those already found and add to PJ. +*** Will try up to MXITER iterations to get a stable set of particles +*** in the cone. + INTEGER MU,N,ITER + LOGICAL PXSAME,PXNEW,OK + LOGICAL NEWLIS(MXTRAK),OLDLIS(MXTRAK) + DOUBLE PRECISION OAXIS(3),NAXIS(3),PNEW(4) + INTEGER MXITER + PARAMETER(MXITER = 30) +* + DO 100 MU=1,3 + OAXIS(MU) = VSEED(MU) +100 CONTINUE + DO 110 N = 1,NTRAK + OLDLIS(N) = .FALSE. +110 CONTINUE + DO 120 ITER = 1,MXITER + CALL PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,PNEW,NEWLIS,OK) +*** Return immediately if there were no particles in the cone. + IF (.NOT.OK) THEN + RETURN + ENDIF + IF(PXSAME(NEWLIS,OLDLIS,NTRAK)) THEN +*** We have a stable jet. + IF (PXNEW(NEWLIS,JETLIS,NTRAK,NJET)) THEN +*** And the jet is a new one. So add it to our arrays. +*** Check arrays are big anough... + IF (NJET .EQ. MXPROT) THEN + WRITE (6,*) ' PXCONE: Found more than MXPROT proto-jets' + IERR = -1 + RETURN + ENDIF + NJET = NJET + 1 + DO 130 N = 1,NTRAK + JETLIS(NJET,N) = NEWLIS(N) +130 CONTINUE + DO 140 MU=1,4 + PJ(MU,NJET)=PNEW(MU) +140 CONTINUE + ENDIF + RETURN + ENDIF +*** The jet was not stable, so we iterate again + DO 150 N=1,NTRAK + OLDLIS(N)=NEWLIS(N) +150 CONTINUE + DO 160 MU=1,3 + OAXIS(MU)=NAXIS(MU) +160 CONTINUE +120 CONTINUE + UNSTBL = .TRUE. + write(*,*) 'MW unstable: ',vseed + RETURN + END +*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper +*-- Author : +C----------------------------------------------------------------------- + SUBROUTINE PXSORV(N,A,K,OPT) +C Sort A(N) into ascending order +C OPT = 'I' : return index array K only +C OTHERWISE : return sorted A and index array K +C----------------------------------------------------------------------- + INTEGER NMAX + PARAMETER (NMAX=500) +* +* INTEGER N,I,J,K(N),IL(NMAX),IR(NMAX) +*LUND + INTEGER N,I,J,K(NMAX),IL(NMAX),IR(NMAX) + CHARACTER OPT +* +* DOUBLE PRECISION A(N),B(NMAX) + DOUBLE PRECISION A(NMAX),B(NMAX) +*LUND + IF (N.GT.NMAX) STOP 'Sorry, not enough room in Mike''s PXSORV' + IL(1)=0 + IR(1)=0 + DO 10 I=2,N + IL(I)=0 + IR(I)=0 + J=1 + 2 IF(A(I).GT.A(J)) GO TO 5 + 3 IF(IL(J).EQ.0) GO TO 4 + J=IL(J) + GO TO 2 + 4 IR(I)=-J + IL(J)=I + GO TO 10 + 5 IF(IR(J).LE.0) GO TO 6 + J=IR(J) + GO TO 2 + 6 IR(I)=IR(J) + IR(J)=I + 10 CONTINUE + I=1 + J=1 + GO TO 8 + 20 J=IL(J) + 8 IF(IL(J).GT.0) GO TO 20 + 9 K(I)=J + B(I)=A(J) + I=I+1 + IF(IR(J)) 12,30,13 + 13 J=IR(J) + GO TO 8 + 12 J=-IR(J) + GO TO 9 + 30 IF(OPT.EQ.'I') RETURN + DO 31 I=1,N + 31 A(I)=B(I) + 999 END + +********************************************************************* +*CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper +*CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper +*-- Author : +* +C+DECK,PXTRY. + SUBROUTINE PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS, + + PNEW,NEWLIS,OK) +* +C+SEQ,DECLARE. + INTEGER MXTRAK + PARAMETER (MXTRAK=900) + INTEGER NTRAK,MODE +*** Note that although PU and PP are assumed to be 2d arrays, they +*** are used as 1d in this routine for efficiency + DOUBLE PRECISION COSR,PU(3*MXTRAK),PP(4*MXTRAK),OAXIS(3),PXMDPI + LOGICAL OK + LOGICAL NEWLIS(MXTRAK) + DOUBLE PRECISION NAXIS(3),PNEW(4), + + PART(3) +*** Finds all particles in cone of size COSR about OAXIS direction. +*** Calculates 4-momentum sum of all particles in cone (PNEW) , and +*** returns this as new jet axis NAXIS (Both unit Vectors) + INTEGER N,MU,NPU,NPP + DOUBLE PRECISION COSVAL,NORMSQ,NORM + DOUBLE PRECISION PHIAXIS,RPAXIS + DOUBLE PRECISION HMW(3*MXTRAK) + COMMON /HMWCOMM/ HMW +* + OK = .FALSE. + DO 100 MU=1,4 + PNEW(MU)=0.0 +100 CONTINUE + NPU=-3 + NPP=-4 + + IF (MODE.eq.3) THEN + P = sqrt(oaxis(1)**2+oaxis(2)**2+oaxis(3)**2) + IF ((P-abs(oaxis(3))).GE. 1D-08) then + RPAXIS = 0.5* log((P+oaxis(3))/(P-oaxis(3))) + PHIAXIS = ATAN2(oaxis(2),oaxis(1)) + ELSE + RPAXIS = 20d0 + PHIAXIS = 0d0 + ENDIF + ENDIF + + DO 110 N=1,NTRAK + NPU=NPU+3 + NPP=NPP+4 + IF (MODE.eq.3) THEN +c PART(1) = PU(1+NPU) +c PART(2) = PU(2+NPU) +c PART(3) = PU(3+NPU) +c COSVAL = 1-MWDIST2(PART,RPAXIS,PHIAXIS) +c write(*,*) '>>> MWcosval1 ',cosval + IF (ABS(HMW(1+NPU)).GE.20.OR.RPAXIS.GE.20) THEN + COSVAL=-1000 + ELSE + COSVAL=1- + + ((RPAXIS-HMW(1+NPU))**2+PXMDPI(PHIAXIS-HMW(2+NPU))**2) + ENDIF +c write(*,*) ' MWcosval2 ',cosval + ELSEIF (MODE.NE.2) THEN + COSVAL=0.0 + DO 120 MU=1,3 + COSVAL=COSVAL+OAXIS(MU)*PU(MU+NPU) +120 CONTINUE + ELSE + IF (ABS(PU(1+NPU)).GE.20.OR.ABS(OAXIS(1)).GE.20) THEN + COSVAL=-1000 + ELSE + COSVAL=1- + + ((OAXIS(1)-PU(1+NPU))**2+PXMDPI(OAXIS(2)-PU(2+NPU))**2) + ENDIF + ENDIF +c - check if particle is inside cone + IF (COSVAL.GE.COSR)THEN + NEWLIS(N) = .TRUE. + OK = .TRUE. + IF (MODE.NE.2) THEN + DO 130 MU=1,4 + PNEW(MU) = PNEW(MU) + PP(MU+NPP) +130 CONTINUE + ELSE + PNEW(1)=PNEW(1) + + + PP(4+NPP)/(PP(4+NPP)+PNEW(4))*(PP(1+NPP)-PNEW(1)) + PNEW(2)=PNEW(2) + + + PP(4+NPP)/(PP(4+NPP)+PNEW(4)) + + *PXMDPI(PP(2+NPP)-PNEW(2)) + PNEW(4)=PNEW(4)+PP(4+NPP) + ENDIF + ELSE + NEWLIS(N)=.FALSE. + ENDIF +110 CONTINUE +*** If there are particles in the cone, calc new jet axis + IF (OK) THEN + IF (MODE.NE.2) THEN + NORMSQ = 0.0 + DO 140 MU = 1,3 + NORMSQ = NORMSQ + PNEW(MU)**2 +140 CONTINUE + NORM = SQRT(NORMSQ) + ELSE + NORM = 1 + ENDIF + DO 150 MU=1,3 + NAXIS(MU) = PNEW(MU)/NORM +150 CONTINUE + ENDIF + RETURN + END + +********************************************************************* +*CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper +*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper +*-- Author : +C+DECK,PXUVEC. +* + SUBROUTINE PXUVEC(NTRAK,PP,PU,IERR) +* +*** Routine to calculate unit vectors PU of all particles PP +C+SEQ,DECLARE. + INTEGER MXTRAK + PARAMETER (MXTRAK=900) + INTEGER NTRAK, IERR + DOUBLE PRECISION PP(4,MXTRAK) + DOUBLE PRECISION PU(3,MXTRAK) + INTEGER N,MU + DOUBLE PRECISION MAG + DO 100 N=1,NTRAK + MAG=0.0 + DO 110 MU=1,3 + MAG=MAG+PP(MU,N)**2 +110 CONTINUE + MAG=SQRT(MAG) + IF (MAG.EQ.0.0) THEN + WRITE(6,*)' PXCONE: An input particle has zero mod(p)' + IERR=-1 + RETURN + ENDIF + DO 120 MU=1,3 + PU(MU,N)=PP(MU,N)/MAG +120 CONTINUE +100 CONTINUE + RETURN + END +*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper +*-- Author : +C----------------------------------------------------------------------- + SUBROUTINE PXZERI(N,A) + INTEGER I,N,A(N) + DO 10 I=1,N + A(I)=0 + 10 CONTINUE + END +*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper +*-- Author : +C----------------------------------------------------------------------- +C This is a set of routines written by Mike Seymour to provide the +C services presumably normally provided by standard OPAL routines +C PXZERV zeroes a vector +C PXZERI zeroes a vector of integers +C PXNORV normalizes a vector +C PXADDV adds two vectors +C PXSORV sorts a vector (copied from HERWIG) +C PXANG3 finds the angle (and its cosine) between two vectors +C PXMDPI moves its argument onto the range [-pi,pi) +C----------------------------------------------------------------------- + SUBROUTINE PXZERV(N,A) + INTEGER I,N + DOUBLE PRECISION A(N) + DO 10 I=1,N + A(I)=0 + 10 CONTINUE + END +*-- Author : +C----------------------------------------------------------------------- + SUBROUTINE PXADDV(N,A,B,C,ITERR) + INTEGER I,N,ITERR + DOUBLE PRECISION A(N),B(N),C(N) + DO 10 I=1,N + C(I)=A(I)+B(I) + 10 CONTINUE + END +*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper +*-- Author : +C----------------------------------------------------------------------- + SUBROUTINE PXANG3(A,B,COST,THET,ITERR) + INTEGER ITERR + DOUBLE PRECISION A(3),B(3),C,COST,THET + C=(A(1)**2+A(2)**2+A(3)**2)*(B(1)**2+B(2)**2+B(3)**2) + IF (C.LE.0) RETURN + C=1/SQRT(C) + COST=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*C + THET=ACOS(COST) + END +*CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper +*-- Author : P. Schleper 28/02/94 + LOGICAL FUNCTION PXNEW(TSTLIS,JETLIS,NTRAK,NJET) +* + INTEGER MXTRAK,MXPROT + PARAMETER (MXTRAK=900,MXPROT=500) + INTEGER NTRAK,NJET +*** Note that although JETLIS is assumed to be a 2d array, it +*** it is used as 1d in this routine for efficiency + LOGICAL TSTLIS(MXTRAK),JETLIS(MXPROT*MXTRAK) +*** Checks to see if TSTLIS entries correspond to a jet already found +*** and entered in JETLIS + INTEGER N, I, IN + LOGICAL MATCH +* + PXNEW = .TRUE. + DO 100 I = 1,NJET + MATCH = .TRUE. + IN=I-MXPROT + DO 110 N = 1,NTRAK + IN=IN+MXPROT + IF(TSTLIS(N).NEQV.JETLIS(IN)) THEN + MATCH = .FALSE. + GO TO 100 + ENDIF +110 CONTINUE + IF (MATCH) THEN + PXNEW = .FALSE. + RETURN + ENDIF +100 CONTINUE + RETURN + END +*CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper +*-- Author : P. Schleper 28/02/94 + LOGICAL FUNCTION PXSAME(LIST1,LIST2,N) +* + LOGICAL LIST1(*),LIST2(*) + INTEGER N +*** Returns T if the first N elements of LIST1 are the same as the +*** first N elements of LIST2. + INTEGER I +* + PXSAME = .TRUE. + DO 100 I = 1,N + IF ( LIST1(I).NEQV.LIST2(I) ) THEN + PXSAME = .FALSE. + RETURN + ENDIF +100 CONTINUE + RETURN + END +*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper +*-- Author : +C----------------------------------------------------------------------- + FUNCTION PXMDPI(PHI) + IMPLICIT NONE +C---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI) + DOUBLE PRECISION PXMDPI,PHI,PI,TWOPI,THRPI,EPS + PARAMETER (PI=3.141592654,TWOPI=6.283185307,THRPI=9.424777961) + PARAMETER (EPS=1E-15) + PXMDPI=PHI + IF (PXMDPI.LE.PI) THEN + IF (PXMDPI.GT.-PI) THEN + GOTO 100 + ELSEIF (PXMDPI.GT.-THRPI) THEN + PXMDPI=PXMDPI+TWOPI + ELSE + PXMDPI=-MOD(PI-PXMDPI,TWOPI)+PI + ENDIF + ELSEIF (PXMDPI.LE.THRPI) THEN + PXMDPI=PXMDPI-TWOPI + ELSE + PXMDPI=MOD(PI+PXMDPI,TWOPI)-PI + ENDIF + 100 IF (ABS(PXMDPI).LT.EPS) PXMDPI=0 + END +C*********************************************************************** +c + FUNCTION MWDIST(vec1,vec2) + IMPLICIT NONE +C--- computes Eta-Phi distance **2 of two four-vectors + DOUBLE PRECISION mwdist, vec1(3),vec2(3), p1,p2, dph, + + PXMDPI + + + P1 = sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2) + IF ((P1-abs(VEC1(3))).LE. 1D-08) then + MWDIST = 1000d0 + return + endif + + P2 = sqrt(vec2(1)**2+vec2(2)**2+vec2(3)**2) + IF ((P2-abs(VEC2(3))).LE. 1D-08) then + MWDIST = 1000d0 + return + endif + + mwdist = + + PXMDPI(ATAN2(VEC1(2),VEC1(1))-ATAN2(VEC2(2),VEC2(1)))**2 + + + 0.25d0* + + log( (P1+VEC1(3))*(P2-VEC2(3))/(P1-VEC1(3))/(P2+VEC2(3)) )**2 + +c write (*,*) 'MWDIST = ',mwdist ,' ',ph1,ph2 + + return + end +C*********************************************************************** diff --git a/Fpmc/External/reco_ok/reco_ok.f b/Fpmc/External/reco_ok/reco_ok.f new file mode 100644 index 0000000..529bfb5 --- /dev/null +++ b/Fpmc/External/reco_ok/reco_ok.f @@ -0,0 +1,340 @@ +* OK: 1/12/2006 interface to the pxcone jet algo. +* 1) convert hep to structure of pxcone(pythia typ) +* 2) do the clustering +* 3) fill the appropriate structure in ntuples + SUBROUTINE MYRECO(IERR) + !IERR = 0 OK + !IERR!= 0 error + + include 'HERWIG65.INC' + integer nfoundjet + + +C ... PXCONE variables, two pxcone algorithms for jet reconstructions + INTEGER ITKDM,MXTRK + ! ... ctyrvektor, max. pocet castic v jednom eventu + PARAMETER (ITKDM=4,MXTRK=4000) + INTEGER MXJET, MXTRAK, MXPROT + PARAMETER (MXJET=20,MXTRAK=4000,MXPROT=500) + INTEGER IPASS (MXTRAK),IJMUL (MXJET) + INTEGER NTRAK,MODE, NJET, IERR + DOUBLE PRECISION PTRAK(ITKDM,MXTRK),PJET(5,MXJET) + DOUBLE PRECISION CONER, EPSLON, OVLIM + + +c ... particles on generator level + integer ngenmax + parameter(ngenmax=1000) + integer ngen + real px(ngenmax),py(ngenmax),pz(ngenmax) + real e(ngenmax),m(ngenmax) + integer id(ngenmax) + + common /gener/ngen, + & px,py,pz,e,m,id + +c ... jets + integer njetmax + parameter(njetmax=30) + integer njets + real + & pxjet(njetmax),pyjet(njetmax), + & pzjet(njetmax),ejet(njetmax) + common/jets/njets, + & pxjet,pyjet, + & pzjet,ejet + + + + data nevhep /0/ + +c ... variables for selecting leptons + real*8 pt, pttresh + integer absid + integer nleptons + INTEGER MXLEP + PARAMETER (MXLEP=20) + REAL*8 PLEP(5,MXLEP) +c ... local variables for misset + real*8 misspx, misspy + + +c ... loop variables + INTEGER I, J, N, IPART + DOUBLE PRECISION RAP + + integer idiv + + +c------------------------------------------------------------------- + +C ... PXCONE parameters ... + MODE = 2 !Snow mass scheme + CONER = 0.7 + EPSLON = 4 +c EPSLON = 8 + OVLIM = 0.5d0 + +c------------------------------------------------------------------- +c particles on the generator level + call vzero(px,ngenmax) + call vzero(py,ngenmax) + call vzero(pz,ngenmax) + call vzero(e,ngenmax) + call vzero(m,ngenmax) + call vzero(id,ngenmax) + + + N=NHEP + IPART=0 + do 1515 I=1,N + IF(ISTHEP(I).EQ.1) THEN + IPART=IPART+1 + px(IPART)=sngl(PHEP(1,I)) + py(IPART)=sngl(PHEP(2,I)) + pz(IPART)=sngl(PHEP(3,I)) + e(IPART) =sngl(PHEP(4,I)) + m(IPART)=sngl(PHEP(5,I)) + id(IPART)=IDHEP(I) +c rm(i)=sngl(PTRAK(5,i)) +c print '(A,4F8.2)','gen:',px(i),py(i),pz(i),e(i) +c print '(A,F8.2,I6)','rm id :',p(i,5),k(i,2) + ENDIF +1515 continue + ngen=IPART + + + idiv = MOD(nevhep, 1000) + if(nevhep.lt.100) then + print *,' *** cone_interface nevhep :', nevhep + else if ( idiv.EQ.0 ) then + print *,' *** cone_interface nevhep :', nevhep + ENDIF + + +c------------------------------------------------------------------- + NTRAK=0 + N=NHEP +c print '(A,I6)', 'NHEP', NHEP + DO 180 I=1,N + !if particle in the final state + IF(ISTHEP(I).EQ.1) THEN + RAP = DABS( 0.5*DLOG( (PHEP(4,I)+PHEP(3,I))/ + . (PHEP(4,I)-PHEP(3,I)) ) ) +c print*,rap + IF(RAP.LT.4) THEN + NTRAK=NTRAK+1 + DO 160 J=1,5 + PTRAK(J,NTRAK)=PHEP(J,I) + 160 CONTINUE + ENDIF + +c print '(A,4F8.2)','phep:',PHEP(1,i),PHEP(2,i), +c . PHEP(3,i),PHEP(4,i) + + ENDIF + 180 CONTINUE + + +c------------------------------------------------------------------- + !clustering + + NJET=0 !midpoint alg. + CALL PXCONE(MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,MXJET, + + NJET,PJET,IPASS,IJMUL,IERR) + + IF(IERR.NE.0) THEN + print *, 'PXCONE did not converge' + RETURN + ENDIF + + + + + +c fill jets + do 1550 i=1, njet + pxjet(i) = PJET(1,i) + pyjet(i) = PJET(2,i) + pzjet(i) = PJET(3,i) + ejet(i) = PJET(4,i) +c print '(A, 4F8.2)', 'pjet x y z e :', PJET(1,i),PJET(2,i), +c . PJET(3,i),PJET(4,i) +c print *,'' + 1550 continue + njets = njet + + + RETURN + END + + + + +********************************************** +* R.S. 05.06.2012 * +* MYRECO copy with additional parameter that * +* decides whether to run jet cone algorithm. * +********************************************** + + SUBROUTINE FPMC_RECO(output,IERR) + !IERR = 0 OK + !IERR!= 0 error + + include 'HERWIG65.INC' + integer output + integer nfoundjet + +C ... PXCONE variables, two pxcone algorithms for jet reconstructions + INTEGER ITKDM,MXTRK + ! ... ctyrvektor, max. pocet castic v jednom eventu + PARAMETER (ITKDM=4,MXTRK=4000) + INTEGER MXJET, MXTRAK, MXPROT + PARAMETER (MXJET=20,MXTRAK=4000,MXPROT=500) + INTEGER IPASS (MXTRAK),IJMUL (MXJET) + INTEGER NTRAK,MODE, NJET, IERR + DOUBLE PRECISION PTRAK(ITKDM,MXTRK),PJET(5,MXJET) + DOUBLE PRECISION CONER, EPSLON, OVLIM + +c ... particles on generator level + integer ngenmax + parameter(ngenmax=1000) + integer ngen + real px(ngenmax),py(ngenmax),pz(ngenmax) + real e(ngenmax),m(ngenmax) + integer id(ngenmax) + integer ii(ngenmax) + integer ist(ngenmax) + + common /gener/ngen, + & px,py,pz,e,m,id,ii,ist + +c ... jets + integer njetmax + parameter(njetmax=30) + integer njets + real + & pxjet(njetmax),pyjet(njetmax), + & pzjet(njetmax),ejet(njetmax) + common/jets/njets, + & pxjet,pyjet, + & pzjet,ejet + + data nevhep /0/ + +c ... variables for selecting leptons + real*8 pt, pttresh + integer absid + integer nleptons + INTEGER MXLEP + PARAMETER (MXLEP=20) + REAL*8 PLEP(5,MXLEP) +c ... local variables for misset + real*8 misspx, misspy + +c ... loop variables + INTEGER I, J, N, IPART + DOUBLE PRECISION RAP + + integer idiv + +c------------------------------------------------------------------- + +C ... PXCONE parameters ... + MODE = 2 !Snow mass scheme + CONER = 0.7 + EPSLON = 4 +c EPSLON = 8 + OVLIM = 0.5d0 + +c------------------------------------------------------------------- +c particles on the generator level + call vzero(px,ngenmax) + call vzero(py,ngenmax) + call vzero(pz,ngenmax) + call vzero(e,ngenmax) + call vzero(m,ngenmax) + call vzero(id,ngenmax) + + N=NHEP + IPART=0 + do 1515 I=1,N + IF((ISTHEP(I).EQ.1.OR.IDHEP(I).EQ.5.OR.IDHEP(I).EQ.-5)) THEN + IPART=IPART+1 + px(IPART)=sngl(PHEP(1,I)) + py(IPART)=sngl(PHEP(2,I)) + pz(IPART)=sngl(PHEP(3,I)) + e(IPART) =sngl(PHEP(4,I)) + m(IPART)=sngl(PHEP(5,I)) + id(IPART)=IDHEP(I) + ii(IPART) = I + ist(IPART) = ISTHEP(I) +c rm(i)=sngl(PTRAK(5,i)) +c print '(A,4F8.2)','gen:',px(i),py(i),pz(i),e(i) +c print '(A,F8.2,I6)','rm id :',p(i,5),k(i,2) + ENDIF +1515 continue + ngen=IPART + + +c idiv = MOD(nevhep, 1000) +c if(nevhep.lt.100) then +c print *,' *** cone_interface nevhep :', nevhep +c else if ( idiv.EQ.0 ) then +c print *,' *** cone_interface nevhep :', nevhep +c ENDIF + + +c------------------------------------------------------------------- + NTRAK=0 + N=NHEP +c print '(A,I6)', 'NHEP', NHEP + DO 180 I=1,N + !if particle in the final state + IF(ISTHEP(I).EQ.1 .and. ABS(IDPDG(I)).GT.100 ) THEN + RAP = DABS( 0.5*DLOG( (PHEP(4,I)+PHEP(3,I))/ + . (PHEP(4,I)-PHEP(3,I)) ) ) + IF(RAP.LT.5) THEN + NTRAK=NTRAK+1 + DO 160 J=1,5 + PTRAK(J,NTRAK)=PHEP(J,I) + 160 CONTINUE + ENDIF + +c print '(A,4F8.2)','phep:',PHEP(1,i),PHEP(2,i), +c . PHEP(3,i),PHEP(4,i) + + ENDIF + 180 CONTINUE + + +c------------------------------------------------------------------- + !clustering + if(output.ge.2) then + + NJET=0 !midpoint alg. + CALL PXCONE(MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,MXJET, + + NJET,PJET,IPASS,IJMUL,IERR) + + IF(IERR.NE.0) THEN + print *, 'PXCONE did not converge' + RETURN + ENDIF + +c fill jets + do 1550 i=1, njet + pxjet(i) = PJET(1,i) + pyjet(i) = PJET(2,i) + pzjet(i) = PJET(3,i) + ejet(i) = PJET(4,i) + 1550 continue + njets = njet + + endif + + + RETURN + END + + + diff --git a/Herwig/herwig.h b/Fpmc/interface/herwig.h similarity index 78% rename from Herwig/herwig.h rename to Fpmc/interface/herwig.h index 70e3aa0..d230edb 100644 --- a/Herwig/herwig.h +++ b/Fpmc/interface/herwig.h @@ -2,16 +2,8 @@ #define HERWIG_INC //-------------------------HERWIG common block ------------------------------------- -// COMMON block stuff, that doesn't come with the HerwigWrapper6_4.h .... -/*C Arrays for particle properties (NMXRES = max no of particles defined) - PARAMETER(NMXRES=500) - COMMON/HWPROP/RLTIM(0:NMXRES),RMASS(0:NMXRES),RSPIN(0:NMXRES), - & ICHRG(0:NMXRES),IDPDG(0:NMXRES),IFLAV(0:NMXRES),NRES, - & VTOCDK(0:NMXRES),VTORDK(0:NMXRES), - & QORQQB(0:NMXRES),QBORQQ(0:NMXRES) */ - -static const int nmxres = 500+1; // we need NMXRES+1 entries ... +static const int nmxres = 500; static const int nmxsud = 1024; // max number of entries for Sudakov lookup table static const int nmxcdk = 4000; static const int nmxhep = 4000; @@ -20,11 +12,6 @@ static const int modmax = 50; static const int maxhrp = 100; static const int imaxch = 20; -static const int NPROC = 117; -static const int MAXMS = 100; -static const int NPSIMP = 16; -static const double SMALL = 1.e-20; - #ifdef __cplusplus extern "C" { @@ -55,7 +42,9 @@ extern "C" } hwevnt_t; extern hwevnt_t hwevnt_; - typedef struct { char PART1[8], PART2[8]; } hwbmch_t; + typedef struct { + char PART1[8], PART2[8]; + } hwbmch_t; extern hwbmch_t hwbmch_; typedef struct { @@ -77,14 +66,19 @@ extern "C" } hwprch_t; extern hwprch_t hwprch_; - + //--- arrays for particle properties (NMXRES = max no of particles defined) typedef struct { - double RLTIM[nmxres], RMASS[nmxres], RSPIN[nmxres]; - int ICHRG[nmxres], IDPDG[nmxres],IFLAV[nmxres], NRES; - int VTOCDK[nmxres], VTORDK[nmxres], QORQQB[nmxres], QBORQQ[nmxres]; + double RLTIM[nmxres+1], RMASS[nmxres+1], RSPIN[nmxres+1]; + int ICHRG[nmxres+1], IDPDG[nmxres+1],IFLAV[nmxres+1], NRES; + int VTOCDK[nmxres+1], VTORDK[nmxres+1], QORQQB[nmxres+1], QBORQQ[nmxres+1]; // starting from 0... } hwprop_t; extern hwprop_t hwprop_; + typedef struct { + unsigned char RNAME[8][nmxres],TXNAME[37][2][nmxres]; + } hwunam_t; + extern hwunam_t hwunam_; + //--- parameters for Sudakov form factors typedef struct { double ACCUR, QEV[6][nmxsud],SUD[6][nmxsud]; @@ -202,41 +196,23 @@ extern "C" extern hw6300_t hw6300_; //-------------------------- JIMMY COMMON BLOCK ------------------------------- -/* - DOUBLE PRECISION YGAMMA, JMZMIN, JMRAD, PTJIM - DOUBLE PRECISION PHAD, JMU2, JMV2, SMALL, JMARRY -c JMARRY is the array storing gamma-p xsec at various z, & -c max weight for each z - DOUBLE PRECISION TOTSCAT, NLOST - - INTEGER MAXMS, NPSIMP, MSFLAG, JMPTYP, JCMVAR, NPROC - LOGICAL ANOMOFF - - PARAMETER( NPROC = 117 ) - PARAMETER( MAXMS = 100 ) ! Maximum multiple scatters - PARAMETER( NPSIMP = 16 ) ! No. of Simpson rule (YBJ) -C intervals (must be even) - PARAMETER( SMALL = 1.0D-20 ) - INTEGER JMOUT, JMBUG, FN_TYPE, NSCAT, JMUEO, MAXMSTRY - PARAMETER(JMOUT = 6) - COMMON / JMPARM / PTJIM, YGAMMA, JMZMIN, JMRAD(264) - & ,PHAD, JMU2, JMV2, JMARRY( 6+MAXMS,0:NPSIMP ) - & ,NLOST, TOTSCAT, ANOMOFF, JCMVAR, JMUEO - & ,JMPTYP(NPROC), JMBUG, FN_TYPE, MSFLAG, MAXMSTRY - DOUBLE PRECISION JMPROC, JMVETO - COMMON / JMEVNT/ JMPROC(NPROC) - &, JMVETO(2,13), NSCAT -*/ + static const int nproc = 117; + static const int maxms = 100; + static const int npsimp = 16; + //static const double small = 1.e-20; + //static const int jmout = 6; typedef struct { - double PTJIM,YGAMMA,JMZMIN,JMRAD[264],PHAD,JMU2,JMV2,JMARRY[NPSIMP+1][6+MAXMS], - NLOST,TOTSCAT; - int ANAMOFF,JCMVAR,JMUEO,JMPTYP[NPROC],JMBUG,FN_TYPE,MSFLAG,MAXMSTRY; + double PTJIM,YGAMMA,JMZMIN,JMRAD[264],PHAD,JMU2,JMV2; + double JMARRY[npsimp+1][6+maxms]; // array storing gamma-p xsec at various z, and max weight for each z + double NLOST,TOTSCAT; + int ANOMOFF; // logical + int JCMVAR,JMUEO,JMPTYP[nproc],JMBUG,FN_TYPE,MSFLAG,MAXMSTRY; } jmparm_t; extern jmparm_t jmparm_; typedef struct { - double JMPROC[NPROC],JMVETO[13][2]; + double JMPROC[nproc],JMVETO[13][2]; int NSCAT; } jmevnt_t; extern jmevnt_t jmevnt_; @@ -266,7 +242,7 @@ C intervals (must be even) #define hwhdecay hwhdecay_ //--------------------------------------------------------------- - void hwuidt_( int* iopt, int* ipdg, int* iwig, char nwig[8] ); + void hwuidt_( int* iopt, int* ipdg, int* iwig, unsigned char nwig[8] ); #define hwuidt hwuidt_ double hwualf_( int *mode, double* scale ); #define hwualf hwualf_ @@ -274,6 +250,8 @@ C intervals (must be even) #define hwuaem hwuaem_ //--------------------------------------------------------------- + void hwudat_(); +#define hwudat hwudat_ void hweini_(); // initialise elementary process #define hweini hweini_ void hwuine_(); // initialise event @@ -282,7 +260,7 @@ C intervals (must be even) #define hwepro hwepro_ void hwigin_(); // initialise other common blocks #define hwigin hwigin_ - void hwusta_( const char*, int ); // make any particle stable + void hwusta_( const unsigned char*, int ); // make any particle stable #define hwusta hwusta_ void hwuinc_(); // compute parameter-dependent constants #define hwuinc hwuinc_ @@ -302,6 +280,8 @@ C intervals (must be even) #define hwmevt hwmevt_ void hwufne_(); // event generation completed, wrap up event, ... #define hwufne hwufne_ + void hwabeg_(); +#define hwabeg hwabeg_ #ifdef __cplusplus } diff --git a/Fpmc/src/Fpmc.cc b/Fpmc/src/Fpmc.cc index 45bff9d..79a7b49 100644 --- a/Fpmc/src/Fpmc.cc +++ b/Fpmc/src/Fpmc.cc @@ -70,7 +70,6 @@ namespace fpmc /*for ( unsigned int i = 0; i < 500; ++i ) { std::cout << "hwprop for particle " << i << ": " << hwprop_.RLTIM[i] << "\t" << hwprop_.RMASS[i] << "\t" << hwprop_.RSPIN[i] << "\t" << hwprop_.ICHRG[i] << "\t" << hwprop_.IDPDG[i] << "\t" << hwprop_.IFLAV[i] << std::endl; }*/ -// exit(0); params_.fetchHWPRAM( hwpram_ ); hwpram_.IPRINT = herwigVerbosity_; diff --git a/Fpmc/src/FpmcParameters.cc b/Fpmc/src/FpmcParameters.cc index 52075ef..7348e75 100644 --- a/Fpmc/src/FpmcParameters.cc +++ b/Fpmc/src/FpmcParameters.cc @@ -7,7 +7,7 @@ namespace fpmc FpmcParameters::FpmcParameters() : std::map( { // particles masses - { "rmass", "0.0" }, { "wmass", "80.425" }, { "hmass", "125.0" }, { "tmass", "174.3" }, { "mst1", "250.0" }, { "msb1", "250.0" }, + { "wmass", "80.425" }, { "hmass", "125.0" }, { "tmass", "174.3" }, { "mst1", "250.0" }, { "msb1", "250.0" }, // kinematics cuts { "ecms", "14000.0" }, { "yjmin", "-6.0" }, { "yjmax", "6.0" }, { "ptmin", "0.0" }, { "ptmax", "1.e8" }, { "emmin", "10.0" }, { "emmax", "1.e8" }, { "ywwmin", "0.0" }, { "ywwmax", "0.1" }, { "q2wwmn", "0.0" }, { "q2wwmx", "4.0" }, diff --git a/Herwig/CMakeLists.txt b/Herwig/CMakeLists.txt index 418482f..1d33021 100644 --- a/Herwig/CMakeLists.txt +++ b/Herwig/CMakeLists.txt @@ -1,5 +1,10 @@ -#include_directories(${PROJECT_SOURCE_DIR}/External ${PROJECT_SOURCE_DIR}/Herwig) +add_library(FpmcHerwig SHARED herwig6500.f) +set_target_properties(FpmcHerwig PROPERTIES LINKER_LANGUAGE CXX POSITION_INDEPENDENT_CODE ON) add_library(HerwigCore OBJECT herwig6500.f) set_target_properties(HerwigCore PROPERTIES LINKER_LANGUAGE CXX POSITION_INDEPENDENT_CODE ON) +#----- installation procedure + +install(TARGETS HerwigCore DESTINATION lib) + diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt new file mode 100644 index 0000000..aad5511 --- /dev/null +++ b/Tests/CMakeLists.txt @@ -0,0 +1,32 @@ +#----- list all test executables + +file(GLOB executables_noroot RELATIVE ${PROJECT_SOURCE_DIR}/Tests *.cc) +file(GLOB executables_root RELATIVE ${PROJECT_SOURCE_DIR}/Tests *.cxx) + +foreach(exec_src ${executables_noroot}) + string(REPLACE ".cc" "" exec_bin ${exec_src}) + add_executable(${exec_bin} ${exec_src}) + set_target_properties(${exec_bin} PROPERTIES EXCLUDE_FROM_ALL true) + target_link_libraries(${exec_bin} ${LIBRARIES}) + #target_link_libraries(${exec_bin} FPMC) + #add_test(${exec_bin} ${exec_bin}) + #set_tests_properties(${exec_bin} PROPERTIES PASS_REGULAR_EXPRESSION "Passed") +endforeach() + +#----- specify the tests requiring ROOT + +find_package(ROOT QUIET) +if(${ROOT_FOUND}) + message(STATUS "ROOT found in ${ROOT_INCLUDE_DIRS}") + include_directories(${ROOT_INCLUDE_DIRS}) + link_directories(${ROOT_LIBRARY_DIR}) + + foreach(exec_src ${executables_root}) + string(REPLACE ".cxx" "" exec_bin ${exec_src}) + add_executable(${exec_bin} ${exec_src}) + set_target_properties(${exec_bin} PROPERTIES EXCLUDE_FROM_ALL true) + target_link_libraries(${exec_bin} FPMC ${ROOT_LIBRARIES}) + endforeach() +else() + message(STATUS "ROOT not found! skipping these tests!") +endif() diff --git a/Tests/dummy_hwaend.h b/Tests/dummy_hwaend.h new file mode 100644 index 0000000..9704e73 --- /dev/null +++ b/Tests/dummy_hwaend.h @@ -0,0 +1,6 @@ +#ifndef Tests_dummy_hwaend_h +#define Tests_dummy_hwaend_h + +void hwaend_() {} + +#endif diff --git a/Tests/test_solibrary.sh b/Tests/test_solibrary.sh new file mode 100755 index 0000000..2b3fb7a --- /dev/null +++ b/Tests/test_solibrary.sh @@ -0,0 +1,7 @@ +#!/bin/sh + +g++ test_standalone.cc -o test_standalone \ + -L../build/ -lFPMC \ + `gsl-config --libs` \ + `clhep-config --libs` \ + -I../Fpmc -I../Herwig diff --git a/Tests/test_standalone.cc b/Tests/test_standalone.cc new file mode 100644 index 0000000..d424bad --- /dev/null +++ b/Tests/test_standalone.cc @@ -0,0 +1,30 @@ +#include "Fpmc.h" +#include "dummy_hwaend.h" + +int main() +{ + fpmc::Fpmc gen; + + //--- parameterise the run + + gen.parameters().setSqrtS( 13.e3 ); // in GeV + gen.parameters().setProcessId( 16008 ); // QED yy->mu,mu + gen.parameters().setProcessType( fpmc::ExclusiveProcess ); + gen.parameters().setInteractionType( fpmc::QED ); + gen.parameters().setIntermediateFlux( fpmc::PhotonPhotonBudnevCoherent ); + gen.parameters().setPtRange( 10. ); // no upper limit + + //--- initialise the run + + gen.initialise(); + + //--- generate some events + + hwevnt_t event; + unsigned int i=0, num_events = 100; + + do { + if ( gen.next( event ) ) i++; + } while ( i