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
1.6.3