agile is hosted by Hepforge, IPPP Durham

FHerwig.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 #include "AGILe/FHerwig/FHerwig.hh"
00003 #include "AGILe/FHerwig/FHerwigWrapper65.hh"
00004 #include "AGILe/Utils.hh"
00005 #include "AGILe/HepMCTools.hh"
00006 
00007 
00008 namespace AGILe {
00009 
00010   namespace {
00011     const size_t AUTPDF_LENGTH = 20;
00012     const int SEED_OFFSET = 1799;
00013   }
00014 
00015 
00016   // Standard constructor
00017   FHerwig::FHerwig()
00018     : _doHadronise(true), _unitWeight(true), _iprocset(false)
00019   {
00020     _myName = "Herwig";
00021 
00023     HepMC::HEPEVT_Wrapper::set_max_number_entries(FC_HWGETHEPEVTSIZE());
00024     HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
00025 
00026     // Set up particle names map: note that they must be 8 characters long
00027     _particleNames[ELECTRON]    = "E-      ";
00028     _particleNames[POSITRON]    = "E+      ";
00029     _particleNames[PROTON]      = "P       ";
00030     _particleNames[ANTIPROTON]  = "PBAR    ";
00031     _particleNames[NEUTRON]     = "N       ";
00032     _particleNames[ANTINEUTRON] = "NBAR    ";
00033     _particleNames[PHOTON]      = "GAMMA   ";
00034     _particleNames[MUON]        = "MU-     ";
00035     _particleNames[ANTIMUON]    = "MU+     ";
00036     // TAU, ANTITAU
00037     _particleNames[NU_E]        = "NU_E    ";
00038     _particleNames[NU_EBAR]     = "NU_EBAR ";
00039     _particleNames[NU_MU]       = "NU_MU   ";
00040     _particleNames[NU_MUBAR]    = "NU_MUBAR";
00041     _particleNames[NU_TAU]      = "NU_TAU  ";
00042     _particleNames[NU_TAUBAR]   = "NU_TAUBR";
00043     _particleNames[PIPLUS]      = "PI+     ";
00044     _particleNames[PIMINUS]     = "PI-     ";
00045     // PIZERO
00046     // PHOTOELECTRON, PHOTOPOSITRON,
00047     // PHOTOMUON,     PHOTOANTIMUON,
00048     // PHOTOTAU,      PHOTOANTITAU,
00049 
00050     // Start counting events at 1.
00051     _nevt = 1;
00052   }
00053 
00054 
00055   // Set up initial state from supplied params
00056   void FHerwig::setGenSpecificInitialState(int p1, double e1, int p2, double e2) {
00057     MSG_DEBUG("Setting initial state...");
00058     // Herwig's initial state particles specification must be 8 chars long
00059     for (size_t i = 0; i < 8; ++i) {
00060       FC_HWBMCH.PART1[i] = _particleNames[ParticleName(p1)].c_str()[i];
00061       FC_HWBMCH.PART2[i] = _particleNames[ParticleName(p2)].c_str()[i];
00062     }
00063     // Set momenta / energies
00064     // NB. Beam ordering is significant (segfault if "wrong") for HERA-like initial states
00065     // Subtraction of (very rough!) beam particle masses to get momenta rather than energies
00066     if (p1 == PROTON || p1 == ANTIPROTON || p1 == NEUTRON) e1 = sqrt(e1*e1 - 1.0);
00067     if (p2 == PROTON || p2 == ANTIPROTON || p2 == NEUTRON) e2 = sqrt(e2*e2 - 1.0);
00068     // And set the common block values
00069     FC_HWPROC.PBEAM1 = e1;
00070     FC_HWPROC.PBEAM2 = e2;
00071   }
00072 
00073 
00074   // Run generator initialization
00075   void FHerwig::initialize() {
00076     Generator::initialize();
00077 
00078     // Use QCD 2->2 as the default process
00079     if (!_iprocset) FC_HWPROC.IPROC = 1500;
00080 
00081     // Print banner and initialise common blocks (requires initial state and IPROC to have been set)
00082     FC_HWIGIN();
00083 
00084     // for (size_t idhw = 1; idhw <= 458; ++idhw) {
00085     //  std::cout << idhw << "  " << FC_HWPROP.IDPDG[idhw] << std::endl;
00086     // }
00087 
00088     // Set banner printout control params
00089     FC_HWPRAM.IPRINT = 1; // Type of info to print
00090     FC_HWPRAM.PRNDEF = 1; // Enable/disable ASCII output
00091 
00092     // Set some system params which HWIGIN has initialised to silly values
00093     FC_HWPROC.MAXEV = 10000000; // Maximum number of events (irrelevant!)
00094     FC_HWEVNT.MAXER = 10000000; // Maximum number of allowed errors
00095     FC_HWEVNT.MAXPR = 0; // Number of events to print
00096 
00097     // Apply stored parameters
00098     _processParams();
00099 
00100     // Compute parameter-dependent constants
00101     FC_HWUINC();
00102 
00103     // Initialise elementary process
00104     FC_HWEINI();
00105 
00106     _initialized = true;
00107   }
00108 
00109 
00110   // Set random number seed
00111   void FHerwig::setSeed(const int value) {
00112     Generator::setSeed(value);
00113     setParam("SEED", value);
00114   }
00115 
00116 
00117   // Gather user-provided parameters, applying only IPROC and the beam params
00118   bool FHerwig::setParam(const string& name, const string& value) {
00119     Generator::setParam(name, value);
00120     if (name == "IPROC") {
00121       // 1610 = gg -> H -> WW; 1706 = qq -> ttbar; 2510 = ttH -> ttWW, ...
00122       MSG_INFO("Setting process code (IPROC) = " << asInt(value));
00123       FC_HWPROC.IPROC = asInt(value);
00124       _iprocset = true;
00125     } else if (name == "PBEAM1") {
00126       MSG_INFO("Setting beam 1 momentum (PBEAM1) = " << asDouble(value));
00127       FC_HWPROC.PBEAM1 = asDouble(value);
00128     } else if (name == "PBEAM2") {
00129       MSG_INFO("Setting beam 2 momentum (PBEAM2) = " << asDouble(value));
00130       FC_HWPROC.PBEAM2 = asDouble(value);
00131     // } else if (name == "EBEAM1") {
00132     //   MSG_INFO("Setting beam 1 energy (EBEAM1) = " << asDouble(value));
00133     //   FC_HWPROC.EBEAM1 = asDouble(value);
00134     // } else if (name == "EBEAM2") {
00135     //   MSG_INFO("Setting beam 2 energy (EBEAM2) = " << asDouble(value));
00136     //   FC_HWPROC.EBEAM2 = asDouble(value);
00137     } else if (name == "LHEF") {
00138       if(FC_HWPROC.IPROC>0) {
00139     FC_HWPROC.IPROC = -FC_HWPROC.IPROC;
00140     _iprocset=true;
00141     MSG_INFO("LHEF meta-param has been set: turn IPROC negative = "<< FC_HWPROC.IPROC);
00142       } else if (!_iprocset) {
00143     FC_HWPROC.IPROC = -1;
00144     _iprocset=true;
00145     MSG_INFO("LHEF meta-param has been set: set IPROC = "<< FC_HWPROC.IPROC);
00146       }
00147       MSG_INFO("LHEF meta-param has been set: calling OPEN *now* for LHE file = " << value);
00148       call_lhefopen(value.c_str());
00149     } else {
00150       _storedParams[name] = value;
00151     }
00152     return SUCCESS;
00153   }
00154 
00155 
00156   // Set parameters.
00157   bool FHerwig::_processParams() {
00158     typedef pair<string,string> strpair;
00159     foreach (const strpair& kv, _storedParams) {
00160       const string& name = kv.first;
00161       const string& value = kv.second;
00162 
00163       // Strings
00164       if (name.find("AUTPDF") != string::npos) {
00165         if (name == "AUTPDF(1)" || name == "AUTPDF") {
00166           MSG_INFO("Setting AUTPDF(1) = " << value);
00167           const size_t nch = (value.size() < AUTPDF_LENGTH) ? value.size() : AUTPDF_LENGTH-1;
00168           for (size_t i = 0; i < AUTPDF_LENGTH; ++i) {
00169             if (i < nch) {
00170               FC_HWPRCH.AUTPDF[0][i] = value.data()[i];
00171             } else {
00172               FC_HWPRCH.AUTPDF[0][i] = ' ';
00173             }
00174           }
00175         }
00176         if (name == "AUTPDF(2)" || name == "AUTPDF") {
00177           MSG_INFO("Setting AUTPDF(2) = " << value);
00178           const size_t nch = (value.size() < AUTPDF_LENGTH) ? value.size() : AUTPDF_LENGTH-1;
00179           for (size_t i = 0; i < AUTPDF_LENGTH; ++i) {
00180             if (i < nch) {
00181               FC_HWPRCH.AUTPDF[1][i] = value.data()[i];
00182             } else {
00183               FC_HWPRCH.AUTPDF[1][i] = ' ';
00184             }
00185           }
00186         }
00187       }
00188 
00189       // Integers (other than IPROC)
00190       else if (name == "SEED") {
00191         const int baseseed = asInt(value);
00192         MSG_INFO("Setting random seeds based on this value (SEED) = " << baseseed);
00193         FC_HWEVNT.NRN[0] = baseseed;
00194         FC_HWEVNT.NRN[1] = baseseed + SEED_OFFSET;
00195         FC_HWHARD.IBRN[0] = baseseed + 2*SEED_OFFSET;
00196         FC_HWHARD.IBRN[1] = baseseed + 3*SEED_OFFSET;
00197       } else if (name == "MAXPR") {
00198         MSG_INFO("Setting max event printouts (MAXPR) = " << asInt(value));
00199         FC_HWEVNT.MAXPR = asInt(value);
00200       } else if (name == "IPRINT") {
00201         MSG_INFO("Setting printout info code (IPRINT) = " << asInt(value));
00202         FC_HWPRAM.IPRINT = asInt(value);
00203       } else if (name == "CLDIR1") {
00204         MSG_INFO("Setting smearing of perturbative quark pT in light cluster fission (CLDIR(1)) = " << asInt(value));
00205         FC_HWPRAM.CLDIR[0] = asInt(value);
00206       } else if (name == "CLDIR2") {
00207         MSG_INFO("Setting smearing of perturbative quark pT in b cluster fission (CLDIR(2)) = " << asInt(value));
00208         FC_HWPRAM.CLDIR[1] = asInt(value);
00209       } else if (name == "IOPREM") {
00210         MSG_INFO("Setting model for treatment of remnant clusters (IOPREM) = " << asInt(value));
00211         FC_HWPRAM.IOPREM = asInt(value);
00212       } else if (name == "ISPAC") {
00213         MSG_INFO("Setting ISR forced branching IR behaviour (ISPAC) = " << asInt(value));
00214         FC_HWPRAM.ISPAC = asInt(value);
00215       } else if (name == "NFLAV") {
00216         MSG_INFO("Setting number of flavours (NFLAV) = " << asInt(value));
00217         FC_HWPRAM.NFLAV = asInt(value);
00218       } else if (name == "NSTRU") {
00219         MSG_INFO("Setting internal PDF ID... are you sure? (NSTRU) = " << asInt(value));
00220         FC_HWPRAM.NSTRU = asInt(value);
00221       } else if (name == "MODPDF") {
00222         MSG_INFO("Setting PDFLIB IDs for both beams (MODPDF(1&2)) = " << asInt(value));
00223         FC_HWPRAM.MODPDF[0] = asInt(value);
00224         FC_HWPRAM.MODPDF[1] = asInt(value);
00225       } else if (name == "MODPDF(1)") {
00226         MSG_INFO("Setting PDFLIB ID for beam 1 (MODPDF(1)) = " << asInt(value));
00227         FC_HWPRAM.MODPDF[0] = asInt(value);
00228       } else if (name == "MODPDF(2)") {
00229         MSG_INFO("Setting PDFLIB ID for beam 2 (MODPDF(2)) = " << asInt(value));
00230         FC_HWPRAM.MODPDF[1] = asInt(value);
00231       } else if (name == "IFLMAX") {
00232         MSG_INFO("Setting max quark flavour in photoproduction (IFLMAX) = " << asInt(value));
00233         FC_HWHARD.IFLMAX = asInt(value);
00234       } else if (name == "LRSUD") {
00235         MSG_INFO("Setting unit for reading Sudakov table (LRSUD) = " << asInt(value));
00236         FC_HWPRAM.LRSUD = asInt(value);
00237       } else if (name == "LWSUD") {
00238         MSG_INFO("Setting unit for writing Sudakov table (LWSUD) = " << asInt(value));
00239         FC_HWPRAM.LWSUD = asInt(value);
00240       } else if (name == "IAPHIG") {
00241         MSG_INFO("Setting gg->H top loop param (IAPHIG) = " << asInt(value));
00242         FC_HWHARD.IAPHIG = asInt(value);
00243       // } else if (name == "MODBOS") {
00244       //   MSG_INFO("Setting MODBOS = " << asInt(value));
00245       //   FC_HWBOSC.MODBOS[0] = asInt(value);
00246 
00247       // Doubles
00248       } else if (name == "ALPHEM") {
00249         MSG_INFO("Setting alpha_EM (ALPHEM) = " << asDouble(value));
00250         FC_HWPRAM.ALPHEM = asDouble(value);
00251       } else if (name == "PTMIN") {
00252         MSG_INFO("Setting minimum hard pt (PTMIN) = " << asDouble(value));
00253         FC_HWHARD.PTMIN = asDouble(value);
00254       } else if (name == "PTMAX") {
00255         MSG_INFO("Setting maximum hard pt (PTMAX) = " << asDouble(value));
00256         FC_HWHARD.PTMAX = asDouble(value);
00257       } else if (name == "PTPOW") {
00258         MSG_INFO("Setting power for reweight (PTPOW) = " << asDouble(value));
00259         FC_HWHARD.PTPOW = asDouble(value);
00260       } else if (name == "Q2WWMN") {
00261         MSG_INFO("Setting minimum Q2 in equivalent photon approximation (Q2WWMN) = " << asDouble(value));
00262         FC_HWHARD.Q2WWMN = asDouble(value);
00263       } else if (name == "Q2WWMX") {
00264         MSG_INFO("Setting maximum Q2 in equivalent photon approximation (Q2WWMX) = " << asDouble(value));
00265         FC_HWHARD.Q2WWMX = asDouble(value);
00266       } else if (name == "YWWMIN") {
00267         MSG_INFO("Setting minimum photon light-cone fraction in equivalent photon approximation (YWWMIN) = " << asDouble(value));
00268         FC_HWHARD.YWWMIN = asDouble(value);
00269       } else if (name == "YWWMAX") {
00270         MSG_INFO("Setting maximum photon light-cone fraction in equivalent photon approximation (YWWMAX) = " << asDouble(value));
00271         FC_HWHARD.YWWMAX = asDouble(value);
00272       } else if ( name =="WHMIN") {
00273         MSG_INFO("Setting min hadronic mass in photon-induced (DIS) processes (WHMIN) = " << asDouble(value));
00274         FC_HWHARD.WHMIN = asDouble(value);
00275       } else if ( name =="EMMIN") {
00276         MSG_INFO("Setting minimum DY boson mass (EMMIN) = "<< asDouble(value));
00277         FC_HWHARD.EMMIN = asDouble(value);
00278       } else if ( name =="EMMAX") {
00279         MSG_INFO("Setting maximum DY boson mass (EMMIN) = "<< asDouble(value));
00280         FC_HWHARD.EMMAX = asDouble(value);
00281       } else if ( name =="EMPOW") {
00282         MSG_INFO("Setting power for EM reweighting (EMPOW) = "<< asDouble(value));
00283         FC_HWHARD.EMPOW = asDouble(value);
00284       } else if ( name =="EMSCA") {
00285         MSG_INFO("Setting ??? for EM (EMSCA) = "<< asDouble(value));
00286         FC_HWHARD.EMSCA = asDouble(value);
00287       } else if (name == "PRECO") {
00288         MSG_INFO("Setting probability of cluster colour reassignment (PRECO) = " << asDouble(value));
00289         FC_HWUCLU.PRECO = asDouble(value);
00290       } else if (name == "PDFX0") {
00291         MSG_INFO("Setting maximum x for saturation effects (PDFX0) = " << asDouble(value));
00292         FC_HW6506.PDFX0 = asDouble(value);
00293       } else if (name == "PDFPOW") {
00294         MSG_INFO("Setting suppression of PDFs below x0 (PDFPOW) = " << asDouble(value));
00295         FC_HW6506.PDFPOW = asDouble(value);
00296       } else if (name == "BTCLM") {
00297         MSG_INFO("Setting remnant cluster mass param adjustment (BTCLM) = " << asDouble(value));
00298         FC_HWPRAM.BTCLM = asDouble(value);
00299       } else if (name == "CLMAX") {
00300         MSG_INFO("Setting cluster mass offset (CLMAX) = " << asDouble(value));
00301         FC_HWPRAM.CLMAX = asDouble(value);
00302       } else if (name == "CLPOW") {
00303         MSG_INFO("Setting exponent of cluster fission (CLPOW) = " << asDouble(value));
00304         FC_HWPRAM.CLPOW = asDouble(value);
00305       } else if (name == "CLSMR1") {
00306         MSG_INFO("Setting Gaussian quark pT smearing width for light clusters (CLSMR(1)) = " << asDouble(value));
00307         FC_HWPRAM.CLSMR[0] = asDouble(value);
00308       } else if (name == "CLSMR2") {
00309         MSG_INFO("Setting Gaussian quark pT smearing width for b clusters (CLSMR(2)) = " << asDouble(value));
00310         FC_HWPRAM.CLSMR[1] = asDouble(value);
00311       } else if (name == "PRSOF") {
00312         MSG_INFO("Setting probablility of soft scatters (PRSOF) = " << asDouble(value));
00313         FC_HWPRAM.PRSOF = asDouble(value);
00314       } else if (name == "PSPLT1") {
00315         MSG_INFO("Setting cluster mass spectrum exponent for light clusters (PSPLT(1)) = " << asDouble(value));
00316         FC_HWPRAM.PSPLT[0] = asDouble(value);
00317       } else if (name == "PSPLT2") {
00318         MSG_INFO("Setting cluster mass spectrum exponent for b clusters (PSPLT(2)) = " << asDouble(value));
00319         FC_HWPRAM.PSPLT[1] = asDouble(value);
00320       }
00321       // Following params are useful for MPI/ISR tunes:
00322       else if (name == "PTRMS") {
00323         MSG_INFO("Setting intrinsic kT (PTRMS) = " << asDouble(value));
00324         FC_HWPRAM.PTRMS = asDouble(value);
00325       } else if (name == "QCDLAM") {
00326         MSG_INFO("Setting Lambda_QCD for alpha_s running (QCDLAM) = " << asDouble(value));
00327         FC_HWPRAM.QCDLAM = asDouble(value);
00328       } else if (name == "QSPAC") {
00329         MSG_INFO("Setting Q cutoff for spacelike (ISR) shower (QSPAC) = " << asDouble(value));
00330         FC_HWPRAM.QSPAC = asDouble(value);
00331       } else if (name == "VGCUT") {
00332         MSG_INFO("Setting gluon virtuality cutoff in parton showers (VGCUT) = " << asDouble(value));
00333         FC_HWPRAM.VGCUT = asDouble(value);
00334       } else if (name == "VQCUT") {
00335         MSG_INFO("Setting quark virtuality cutoff in parton showers (VQCUT) = " << asDouble(value));
00336         FC_HWPRAM.VQCUT = asDouble(value);
00337       } else if (name == "ZMXISR") {
00338         MSG_INFO("Setting max. momentum fraction for photon ISR (ZMXISR) = " << asDouble(value));
00339         FC_HWHARD.ZMXISR = asDouble(value);
00340       }
00341       // End MPI/ISR tune params
00342       else if (name == "VPCUT") {
00343         MSG_INFO("Setting photon virtuality cutoff (VPCUT) = " << asDouble(value));
00344         FC_HWPRAM.VPCUT = asDouble(value);
00345       } else if (name == "OMEGA0") {
00346         MSG_INFO("Setting omega_0 param in Mueller-Tang pomeron for IPROC=2400 (OMEGA0) = " << asDouble(value));
00347         FC_HWHARD.OMEGA0 = asDouble(value);
00348       } else if (name == "ASFIXD") {
00349         MSG_INFO("Setting a_s param in Mueller-Tang pomeron for IPROC=2400 (ASFIXD) = " << asDouble(value));
00350         FC_HWHARD.ASFIXD = asDouble(value);
00351       } else if (name == "Q2MIN") {
00352         MSG_INFO("Setting minimum DIS Q2 (Q2MIN) = " << asDouble(value));
00353         FC_HWHARD.Q2MIN = asDouble(value);
00354       } else if (name == "Q2MAX") {
00355         MSG_INFO("Setting maximum DIS Q2 (Q2MAX) = " << asDouble(value));
00356         FC_HWHARD.Q2MAX = asDouble(value);
00357       } else if (name == "Q2POW") {
00358         MSG_INFO("Setting reweighting power in DIS Q2 (Q2POW) = " << asDouble(value));
00359         FC_HWHARD.Q2POW = asDouble(value);
00360       } else if (name == "Q2WWMN") {
00361         MSG_INFO("Setting Q2WWMN (Q2WWMN) = " << asDouble(value));
00362         FC_HWHARD.Q2WWMN = asDouble(value);
00363       } else if (name == "Q2WWMX") {
00364         MSG_INFO("Setting Q2WWMX (Q2WWMX) = " << asDouble(value));
00365         FC_HWHARD.Q2WWMX = asDouble(value);
00366       } else if (name == "QLIM") {
00367         MSG_INFO("Setting QLIM (QLIM) = " << asDouble(value));
00368         FC_HWHARD.QLIM = asDouble(value);
00369       } else if (name == "YBMIN") {
00370         MSG_INFO("Setting minimum Bjorken y=Q2/xs (YBMIN) = " << asDouble(value));
00371         FC_HWHARD.YBMIN = asDouble(value);
00372       } else if (name == "YBMAX") {
00373         MSG_INFO("Setting maximum Bjorken y=Q2/xs (YBMAX) = " << asDouble(value));
00374         FC_HWHARD.YBMAX = asDouble(value);
00375       } else if (name == "YJMIN") {
00376         MSG_INFO("Setting minimum ??? (YJMIN) = " << asDouble(value));
00377         FC_HWHARD.YJMIN = asDouble(value);
00378       } else if (name == "YJMAX") {
00379         MSG_INFO("Setting maximum ??? (YJMAX) = " << asDouble(value));
00380         FC_HWHARD.YJMAX = asDouble(value);
00381       } else if (name == "YWWMIN") {
00382         MSG_INFO("Setting minimum ??? (YWWMIN) = " << asDouble(value));
00383         FC_HWHARD.YWWMIN = asDouble(value);
00384       } else if (name == "YWWMAX") {
00385         MSG_INFO("Setting maximum ??? (YWWMAX) = " << asDouble(value));
00386         FC_HWHARD.YWWMAX = asDouble(value);
00387       } else if (name == "XX1") {
00388         MSG_INFO("Setting ??? (XX(1)) = " << asDouble(value));
00389         FC_HWHARD.XX[0] = asDouble(value);
00390       } else if (name == "XX2") {
00391         MSG_INFO("Setting ??? (XX(2)) = " << asDouble(value));
00392         FC_HWHARD.XX[1] = asDouble(value);
00393       } else if (name == "XLMIN") {
00394         MSG_INFO("Setting minimum ??? (XLMIN) = " << asDouble(value));
00395         FC_HWHARD.XLMIN = asDouble(value);
00396       } else if (name == "XXMIN") {
00397         MSG_INFO("Setting minimum ??? (XXMIN) = " << asDouble(value));
00398         FC_HWHARD.XXMIN = asDouble(value);
00399       } else if (name == "WHMIN") {
00400         MSG_INFO("Setting minimum ??? (WHMIN) = " << asDouble(value));
00401         FC_HWHARD.WHMIN = asDouble(value);
00402       } else if (name == "ZJMAX") {
00403         MSG_INFO("Setting maximum ??? (ZJMAX) = " << asDouble(value));
00404         FC_HWHARD.ZJMAX = asDouble(value);
00405       } else if (name == "TMNISR") {
00406         MSG_INFO("Setting minimum ??? (TMNISR) = " << asDouble(value));
00407         FC_HWHARD.TMNISR = asDouble(value);
00408       } else if (name == "ZMXISR") {
00409         MSG_INFO("Setting maximum ??? in ISR (ZMXISR) = " << asDouble(value));
00410         FC_HWHARD.ZMXISR = asDouble(value);
00411       } else if (name == "THMAX") {
00412         MSG_INFO("Setting maximum ??? (THMAX) = " << asDouble(value));
00413         FC_HWHARD.THMAX = asDouble(value);
00414       } else if (name == "Y4JT") {
00415         MSG_INFO("Setting ??? (Y4JT) = " << asDouble(value));
00416         FC_HWHARD.Y4JT = asDouble(value);
00417       } else if (name == "SINS") {
00418         MSG_INFO("Setting ??? (SINS) = " << asDouble(value));
00419         FC_HWHARD.SINS = asDouble(value);
00420       } else if (name == "PLTCUT") {
00421         MSG_INFO("Setting lifetime cut (PLTCUT) = " << asDouble(value));
00422         FC_HWDIST.PLTCUT = asDouble(value);
00423       } else if (name == "GAMW") {
00424         MSG_INFO("Setting W width (GAMW) = " << asDouble(value));
00425         FC_HWPRAM.GAMW = asDouble(value);
00426       } else if (name == "GAMZ") {
00427         MSG_INFO("Setting Z0 width (GAMZ) = " << asDouble(value));
00428         FC_HWPRAM.GAMZ = asDouble(value);
00429       } else if (name == "GAMZP") {
00430         MSG_INFO("Setting Z0/gamma width (GAMZP) = " << asDouble(value));
00431         FC_HWPRAM.GAMZP = asDouble(value);
00432       } else if (name == "H1MIX") {
00434         MSG_INFO("Setting H1MIX (H1MIX) = " << asDouble(value));
00435         FC_HWPRAM.H1MIX = asDouble(value);
00436       } else if (name == "PHIMIX") {
00438         MSG_INFO("Setting PHIMIX (PHIMIX) = " << asDouble(value));
00439         FC_HWPRAM.PHIMIX = asDouble(value);
00440       }
00441       // NB. The -1 C-Fortran array offset does not apply for RMASS, which is declared as (0:NMXRES).
00442       else if (name == "RMASS(1)") {
00443         MSG_INFO("Setting down mass (RMASS(1)) = " << asDouble(value));
00444         FC_HWPROP.RMASS[1] = asDouble(value);
00445       } else if (name == "RMASS(2)") {
00446         MSG_INFO("Setting up mass (RMASS(2)) = " << asDouble(value));
00447         FC_HWPROP.RMASS[2] = asDouble(value);
00448       } else if (name == "RMASS(3)") {
00449         MSG_INFO("Setting strange mass (RMASS(3)) = " << asDouble(value));
00450         FC_HWPROP.RMASS[3] = asDouble(value);
00451       } else if (name == "RMASS(4)") {
00452         MSG_INFO("Setting charm mass (RMASS(4)) = " << asDouble(value));
00453         FC_HWPROP.RMASS[4] = asDouble(value);
00454       } else if (name == "RMASS(5)") {
00455         MSG_INFO("Setting bottom mass (RMASS(5)) = " << asDouble(value));
00456         FC_HWPROP.RMASS[5] = asDouble(value);
00457       } else if (name == "RMASS(6)") {
00458         MSG_INFO("Setting top mass (RMASS(6)) = " << asDouble(value));
00459         FC_HWPROP.RMASS[6] = asDouble(value);
00460       } else if (name == "RMASS(198)") {
00461         MSG_INFO("Setting W+ mass (RMASS(198)) = " << asDouble(value));
00462         FC_HWPROP.RMASS[198] = asDouble(value);
00463       } else if (name == "RMASS(199)") {
00464         MSG_INFO("Setting W- mass (RMASS(199)) = " << asDouble(value));
00465         FC_HWPROP.RMASS[199] = asDouble(value);
00466       } else if (name == "RMASS(200)") {
00467         MSG_INFO("Setting Z0 mass (RMASS(200)) = " << asDouble(value));
00468         FC_HWPROP.RMASS[200] = asDouble(value);
00469       } else if (name == "RMASS(201)") {
00470         MSG_INFO("Setting H mass (RMASS(201)) = " << asDouble(value));
00471         FC_HWPROP.RMASS[201] = asDouble(value);
00472       } //...
00473 
00474       // Booleans (done as ints in Fortran)
00475       else if (name == "GENSOF") {
00476         MSG_INFO("Setting ... (GENSOF) = " << asInt(value));
00477         FC_HWEVNT.GENSOF = asInt(value);
00478       } else if (name == "AZSOFT") {
00479         MSG_INFO("Setting use of soft gluon azimuthal corellations on/off (AZSOFT) = " << asInt(value));
00480         FC_HWPRAM.AZSOFT = asInt(value);
00481       } else if (name == "CLRECO") {
00482         MSG_INFO("Setting inclusion of colour rearrangement on/off (CLRECO) = " << asInt(value));
00483         FC_HWUCLU.CLRECO = asInt(value);
00484       } else if (name == "AZSPIN") {
00485         MSG_INFO("Setting use of gloun spin azimuthal correlations on/off (AZSPIN) = " << asInt(value));
00486         FC_HWPRAM.AZSPIN = asInt(value);
00487       } else if (name == "ZPRIME") {
00488         MSG_INFO("Setting Z' flag (ZPRIME) = " << asInt(value));
00489         FC_HWPRAM.ZPRIME = asInt(value);
00490       } else if (name == "HARDME") {
00491         MSG_INFO("Setting hard matrix element flag (HARDME) = " << asInt(value));
00492         FC_HWPRAM.HARDME = asInt(value);
00493       } else if (name == "SOFTME") {
00494         MSG_INFO("Setting soft matrix element flag (SOFTME) = " << asInt(value));
00495         FC_HWPRAM.SOFTME = asInt(value);
00496       } else if (name == "NOSPAC") {
00497         MSG_INFO("Setting spacelike showers on/off (NOSPAC) = " << asInt(value));
00498         FC_HWPRAM.NOSPAC = asInt(value);
00499       } else if (name == "PRVTX") {
00500         MSG_INFO("Setting inclusion of vertex info in event printout on/off (PRVTX) = " << asInt(value));
00501         FC_HWPRAM.PRVTX = asInt(value);
00502       } else if (name == "PRNDEC") {
00503         MSG_INFO("Setting use of decimal in event printout (PRNDEC) = " << asInt(value));
00504         FC_HWPRAM.PRNDEC = asInt(value);
00505       } else if (name == "PRNDEF") {
00506         MSG_INFO("Setting stdout printout on/off (PRNDEF) = " << asInt(value));
00507         FC_HWPRAM.PRNDEF = asInt(value);
00508       } else if (name == "PRNTEX") {
00509         MSG_INFO("Setting LaTeX output on/off (PRNTEX) = " << asInt(value));
00510         FC_HWPRAM.PRNTEX = asInt(value);
00511       } else if (name == "PRNWEB") {
00512         MSG_INFO("Setting HTML output on/off (PRNWEB) = " << asInt(value));
00513         FC_HWPRAM.PRNWEB = asInt(value);
00514       } else if (name == "MIXING") {
00515         MSG_INFO("Setting neutral B mixing on/off (MIXING) = " << asInt(value));
00516         FC_HWDIST.MIXING = asInt(value);
00517       } else if (name == "NOWGT") {
00518         MSG_INFO("Setting unweighted generation (NOWGT) = " << asInt(value));
00519         FC_HWEVNT.NOWGT = asInt(value);
00520       } else if (name == "DoHadronisation") {
00521         MSG_INFO("Do hadronisation = " << value);
00522         _doHadronise = asBool(value);//...
00523       } else if (name == "unitWeight") {
00524         MSG_INFO("Using a weight value of 1 in unweighted events = " << value);
00525         _unitWeight = asBool(value);
00526         // Error
00527       } else {
00528         MSG_ERROR("Herwig doesn't have a parameter called " << name);
00529         return FAILURE;
00530       }
00531 
00532     }
00533     return SUCCESS;
00534   }
00535 
00536 
00537   // Run the generator for one event
00538   void FHerwig::makeEvent(HepMC::GenEvent& evt) {
00539     // Loop until event works
00540     while (true) {
00541       Generator::makeEvent(evt);
00542       FC_HWUINE();  // Initialize event
00543       FC_HWEPRO();  // Generate hard subprocess
00544       FC_HWBGEN();  // Generate parton cascade
00545       FC_HWDHOB();  // Do heavy quark decays
00546       if(_doHadronise){
00547         FC_HWCFOR();  // Do cluster formation
00548         FC_HWCDEC();  // Do cluster decays
00549         FC_HWDHAD();  // Do unstable particle decays
00550         FC_HWDHVY();  // Do heavy flavor decays
00551         FC_HWMEVT();  // Add soft underlying event
00552       }
00553       FC_HWUFNE();  // Finish event
00554       if (FC_HWEVNT.IERROR == 0) break;
00555     }
00556     // Fill event from HEPEVT
00557     fillEvent(evt);
00558     // Increment an event counter (Pythia does not count for itself).
00559     _nevt++;
00560   }
00561 
00562 
00564   void FHerwig::fillEvent(HepMC::GenEvent& evt) {
00565     HepMC::IO_HERWIG hepevt;
00566     hepevt.fill_next_event(&evt);
00567     fixHepMCUnitsFromGeVmm(evt);
00568     evt.set_event_number(_nevt);
00569 
00570     // Set beam particle status = 4
00571     if (evt.valid_beam_particles()) {
00572       evt.beam_particles().first->set_status(4);
00573       evt.beam_particles().second->set_status(4);
00574     }
00575 
00576     // Rewrite status codes to match HepMC standard (2 -> 3, 195-199 -> 2)
00577     for (HepMC::GenEvent::particle_const_iterator p = evt.particles_begin(); p != evt.particles_end(); ++p) {
00578       if ((*p)->status() == 2) (*p)->set_status(3);
00579       if ((*p)->status() >= 195 && (*p)->status() <= 199) (*p)->set_status(2);
00580     }
00581 
00582     //IO_HERWIG::fillEvent does not fill the weight
00583     evt.weights().clear();
00584     if(FC_HWEVNT.NOWGT == 0 || !_unitWeight){
00585       evt.weights().push_back(FC_HWEVNT.EVWGT);
00586     }else{
00587       evt.weights().push_back(1.0);
00588     }
00589 
00590     // Cross-section
00591     #ifdef HEPMC_HAS_CROSS_SECTION
00592     HepMC::GenCrossSection xsec;
00593     const double xsecval = getCrossSection();
00594     const double xsecerr = getCrossSection() / std::sqrt(_nevt);
00595     MSG_DEBUG("Writing cross-section = " << xsecval << " +- " << xsecerr);
00596     xsec.set_cross_section(xsecval, xsecerr);
00597     evt.set_cross_section(xsec);
00598     #endif
00599   }
00600 
00601 
00602   // Return the beam number (0 or 1) of a particle.
00603   // Return -1 if not a beam particle
00604   int FHerwig::_beamNumber(PdgCode pid) {
00605     string pName = _particleNames[pid];
00606     string beamName = "";
00607     for (int i = 0; i < 8; ++i) {
00608       beamName += FC_HWBMCH.PART1[i];
00609     }
00610     if (beamName.compare(pName) == 0) return 0;
00611 
00612     for (int i = 0; i < 8; ++i) {
00613       beamName += FC_HWBMCH.PART2[i];
00614     }
00615     if (beamName.compare(pName) == 0) return 1;
00616     return -1;
00617   }
00618 
00619 
00620   string FHerwig::getPDFSet(PdgCode pid){
00621     const int beamNum = _beamNumber(pid);
00622     if (beamNum < 0) {
00623       throw runtime_error("PDFSet unknown for PDG code " + pid);
00624     }
00625     string autString = "";
00626     for (size_t ii = 0; ii != AUTPDF_LENGTH; ++ii) {
00627       autString += FC_HWPRCH.AUTPDF[beamNum][ii];
00628     }
00629     return autString;
00630   }
00631 
00632 
00633   int FHerwig::getPDFMember(PdgCode pid){
00634     const int beamNum = _beamNumber(pid);
00635     if (beamNum < 0) {
00636       throw runtime_error("PDFMember unknown for PDG code "+ pid);
00637     }
00638     const int member = FC_HWPRAM.MODPDF[beamNum];
00639     return member;
00640   }
00641 
00642 
00643   string FHerwig::getPDFScheme(PdgCode pid) const {
00644     return "LHAPDF";
00645   }
00646 
00647 
00648   const double FHerwig::getCrossSection(){
00649     _crossSection = 1000.0 * FC_HWEVNT.WGTSUM / ((float)FC_HWEVNT.NWGTS);
00650     return _crossSection;
00651   }
00652 
00653 
00655   void FHerwig::finalize() {
00656     MSG_DEBUG("Finalising...");
00657     FC_HWEFIN();
00658   }
00659 
00660 
00661 }
00662 
00663 
00664 // Class factory
00665 extern "C" {
00666   AGILe::Generator* create() { return new AGILe::FHerwig(); }
00667   void destroy(AGILe::Generator* gen) { delete gen; }
00668 }
Generated on Tue Mar 6 10:39:38 2012 for AGILe - A Generator Interface Library (+ executable) by  doxygen 1.6.3