diff --git a/Contrib/AlpGen/AlpGenHandler.cc b/Contrib/AlpGen/AlpGenHandler.cc deleted file mode 100644 --- a/Contrib/AlpGen/AlpGenHandler.cc +++ /dev/null @@ -1,1392 +0,0 @@ -#include "AlpGenHandler.h" -#include "ThePEG/Interface/ClassDocumentation.h" -#include "ThePEG/Interface/Switch.h" -#include "ThePEG/Interface/Reference.h" -#include "ThePEG/Repository/UseRandom.h" -#include "ThePEG/Interface/Parameter.h" -#include "ThePEG/Persistency/PersistentOStream.h" -#include "ThePEG/Persistency/PersistentIStream.h" -#include "Herwig/Shower/QTilde/Base/PartnerFinder.h" -#include "Herwig/PDF/HwRemDecayer.h" -#include -#include "ThePEG/Utilities/Throw.h" -#include "Herwig/Shower/QTilde/Kinematics/KinematicsReconstructor.h" -#include "fastjet/PseudoJet.hh" -#include "fastjet/ClusterSequence.hh" -#include "gsl/gsl_rng.h" -#include "gsl/gsl_randist.h" - - -using namespace Herwig; - -bool recordEntry(PPtr i,PPtr j) { - return (i->number()number()); -} -bool pTsortFunction(PPtr i,PPtr j) { - return (i->momentum().perp2()>j->momentum().perp2()); -} -bool ETsortFunction(pair i,pair j) { - return (i.first>j.first); -} -bool isMomLessThanEpsilon(Lorentz5Momentum p,Energy epsilon) { - return (abs(p.x())> alphaS_ - >> ncy_ >> ncphi_ >> ihvy_ >> nph_ >> nh_ - >> iunit(etclusmean_,GeV) >> rclus_ >> etaclmax_ >> rclusfactor_ - >> ihrd_ >> njets_ >> drjmin_ >> highestMultiplicity_ - >> ycmax_ >> ycmin_ >> jetAlgorithm_ >> vetoIsTurnedOff_ - >> inputIsNLO_ >> highestNLOMultiplicity_ >> etclusfixed_ - >> cphcal_ >> sphcal_ >> cthcal_ >> sthcal_ >> iunit(epsetclus_,GeV); -} - -ClassDescription AlpGenHandler::initAlpGenHandler; -// Definition of the static class description member. - -void AlpGenHandler::Init() { - - static ClassDocumentation documentation - ("The AlpGenHandler class performs MEPS merging " - "using the MLM procedure."); - - static Reference interfaceShowerAlpha - ("ShowerAlpha", - "The object calculating the strong coupling constant", - &AlpGenHandler::alphaS_, false, false, true, false, false); - - static Parameter interfaceNoCellsInRapidity - ("NoCellsInRapidity", - "The number of cells spanning the rapidity interval of " - "the calorimeter", - &AlpGenHandler::ncy_, 100, 1, 10000, - false, false, Interface::limited); - - static Parameter interfaceNoCellsInPhi - ("NoCellsInPhi", - "The number of cells spanning the phi interval of " - "the calorimeter", - &AlpGenHandler::ncphi_, 60, 1, 10000, - false, false, Interface::limited); - - static Parameter interfaceihvy - ("ihvy", - "heavy flavour in WQQ,ZQQ,2Q etc (4=c, 5=b, 6=t)", - &AlpGenHandler::ihvy_, -999, -999, 7, - false, false, Interface::limited); - - static Parameter interfacenph - ("nph", - "Number of photons in the AlpGen process", - &AlpGenHandler::nph_, -999, -999, 7, - false, false, Interface::limited); - - static Parameter interfacenh - ("nh", - "Number of higgses in the AlpGen process", - &AlpGenHandler::nph_, -999, -999, 7, - false, false, Interface::limited); - - static Parameter interfaceETClus - ("ETClus", - "The ET threshold defining a jet in the merging procedure", - &AlpGenHandler::etclusmean_, GeV, 20*GeV, 0*GeV, 14000*GeV, - false, false, Interface::limited); - - static Parameter interfaceRClus - ("RClus", - "The cone size used to define a jet in the merging procedure", - &AlpGenHandler::rclus_, 0.4, 0.0, 4.0, - false, false, Interface::limited); - - static Parameter interfaceEtaClusMax - ("EtaClusMax", - "The maximum |eta| used to define a jet in the merging procedure", - &AlpGenHandler::etaclmax_, 5.0, 0.0, 15.0, - false, false, Interface::limited); - - static Parameter interfaceRClusFactor - ("RClusFactor", - "The prefactor for RClus used to define the jet-parton matching " - "distance", - &AlpGenHandler::rclusfactor_, 1.5, 0.0, 4.0, - false, false, Interface::limited); - - static Parameter interfaceihrd - ("ihrd", - "The AlpGen hard process code", - &AlpGenHandler::ihrd_, -999, 0, 10000, - false, false, Interface::limited); - - static Parameter interfacenjets - ("njets", - "The number of light jets in the AlpGen process (i.e. the " - "extra ones)", - &AlpGenHandler::njets_, -999, 0, 10000, - false, false, Interface::limited); - - static Parameter interfacedrjmin - ("drjmin", - "Mimimum parton-parton R-sep used for generation.", - &AlpGenHandler::drjmin_, 0.7, 0.0, 4.0, - false, false, Interface::limited); - - static Parameter interfacehighestMultiplicity - ("highestMultiplicity", - "If true it indicates that this is the highest multiplicity input " - "ME-level configuration to be processed.", - &AlpGenHandler::highestMultiplicity_, 0, 0, 1, - false, false, Interface::limited); - - static Parameter interfacehighestNLOMultiplicity - ("highestNLOMultiplicity", - "If true it indicates that this is the highest NLO multiplicity input " - "ME-level configuration to be processed.", - &AlpGenHandler::highestNLOMultiplicity_, 0, 0, 1, - false, false, Interface::limited); - - static Parameter interfaceETClusFixed - ("ETClusFixed", - "If false, indicates that the jet merging scale, etclus_ is allowed to vary" - "according to epsetclus_", - &AlpGenHandler::etclusfixed_, 1, 0, 1, - false, false, Interface::limited); - - static Parameter interfaceEpsilonETClus - ("EpsilonETClus", - "The ET threshold defining a jet in the merging procedure", - &AlpGenHandler::epsetclus_, GeV, 2.5*GeV, 0*GeV, 100.0*GeV, - false, false, Interface::limited); - - static Switch interfaceJetAlgorithm - ("JetAlgorithm", - "Determines the jet algorithm for finding jets in parton-jet " - "matching in the MLM procedure.", - &AlpGenHandler::jetAlgorithm_, 2, false, false); - static SwitchOption AntiKt - (interfaceJetAlgorithm, - "AntiKt", - "The anti-kt jet algorithm.", - -1); - static SwitchOption CambridgeAachen - (interfaceJetAlgorithm, - "CambridgeAachen", - "The Cambridge-Aachen jet algorithm.", - 0); - static SwitchOption Kt - (interfaceJetAlgorithm, - "Kt", - "The Kt jet algorithm.", - 1); - static SwitchOption GetJet - (interfaceJetAlgorithm, - "GetJet", - "Calorimeter-based GetJet algorithm (default).", - 2); - - static Switch interfaceVetoIsTurnedOff - ("VetoIsTurnedOff", - "Allows the vetoing mechanism to be switched off.", - &AlpGenHandler::vetoIsTurnedOff_, false, false, false); - static SwitchOption VetoingIsOn - (interfaceVetoIsTurnedOff, - "VetoingIsOn", - "The MLM merging veto mechanism is switched ON.", - false); - static SwitchOption VetoingIsOff - (interfaceVetoIsTurnedOff, - "VetoingIsOff", - "The MLM merging veto mechanism is switched OFF.", - true); - - static Switch interfaceInputIsNLO - ("InputIsNLO", - "Signals whether the input LH file is tree-level accurate " - "or contains NLO (Powheg) events.", - &AlpGenHandler::inputIsNLO_, false, false, false); - static SwitchOption InputIsNotNLO - (interfaceInputIsNLO, - "InputIsNotNLO", - "The input LH events have tree-level accuracy.", - false); - static SwitchOption InputIsNLO - (interfaceInputIsNLO, - "InputIsNLO", - "The input LH events have NLO accuracy.", - true); - -} - -void AlpGenHandler::dofinish() { - ShowerHandler::dofinish(); -} - -void AlpGenHandler::doinit() { - - //print error if HardProcID is not set in input file - if(ihrd_ == -999) { cout << "Error: AlpGenHandler:ihrd not set!" << endl; exit(1); } - ShowerHandler::doinit(); - - // Compute calorimeter edges in rapidity for GetJet algorithm. - ycmax_=etaclmax_+rclus_; - ycmin_=-ycmax_; - - // Initialise calorimeter. - calini_m(); -} - -// Throws a veto according to MLM strategy ... when we finish writing it. -bool AlpGenHandler::showerHardProcessVeto() const { - if(vetoIsTurnedOff_) return false; - // Skip veto for processes in which merging is not implemented: - if(ihrd_==7||ihrd_==8||ihrd_==13) { - ostringstream wstring; - wstring << "AlpGenHandler::showerHardProcessVeto() - warning." - << "MLM merging not implemented for AlpGen " - << "processes 4Q (ihrd=7), QQh (ihrd=8), " - << "(single) top (ihrd=13) \n"; - generator()->logWarning( Exception(wstring.str(), - Exception::warning) ); - return false; - } - - // Fill preshowerISPs_ pair and preshowerFSPs_ particle pointer vector. - getPreshowerParticles(); - - // Fill showeredISHs_, showeredISPs and showeredRems pairs, as well as - // showeredFSPs_ particle pointer vector. - getShoweredParticles(); - - // Turn on some screen output debugging: 0 = none ---> 5 = very verbose. - doSanityChecks(0); - - // Dimensions of each calorimter cell in y and phi. - dely_ = (ycmax_-ycmin_)/double(ncy_); - delphi_ = 2*M_PI/double(ncphi_); - - // Fill partonsToMatch_ with only those pre-shower partons intended to - // used in jet-parton matching and fill particlesToCluster_ using only - // those final state particles (post-shower) which are supposed to go - // in the jet clustering used to do merging. - partonsToMatch_ = preshowerFSPs_; - particlesToCluster_ = showeredFSPs_ ; // <--- TO DO: add remnants in here ??? - - // Filter out all but the 'extra' light-parton progenitors and their - // associated final state particles. - caldel_m(); - - double prob(1); - //if etclusfixed_ then set the etclus_ to the fixed chosen value - if(etclusfixed_) { - etclus_ = etclusmean_; - } else { - //else, if we wish to vary etclus_, we use the probability distribution - //choose a probability between 0 and 1 - prob = rnd(); - etclus_ = etclusran_(prob); - } - - if(jetAlgorithm_==2) { - // If using GetJet fill the calorimeter cells now from particlesToCluster_ - calsim_m(); - // Make jets from the calorimeter blobs. - getjet_m(rclus_,etclus_,etaclmax_); - } else { - // Cluster particlesToCluster_ into jets with FastJet. - getFastJets(rclus_,etclus_,etaclmax_); - } - - // If there are less jets than partons then parton-jet matching is - // bound to fail: reject the event already. Also, if the input is - // an NLO event file it will 99.5% of the time contain a number of - // light partons in the F.S. equal to that in the real emission - // process in the NLO calculation, moreover, it has already - // effectively merged njets_-1 and njets jet events. So in that - // case we do not reject events on the grounds they have jet - // multiplicity less than partonsToMatch_.size() but rather less - // jets than partonsToMatch.size()-1; such events are better - // described by the lower-by-one-unit (partonsToMatch_.size()-1) - // of multiplicity NLO event file, or the lower-by-two-units - // (partonsToMatch_.size()-2) of multiplicity LO event file. - - // If it is not jet production apply rejection criterion as above. - if(ihrd_!=9) { - if(!inputIsNLO_) { - if(pjet_.size() < partonsToMatch_.size()) return true; - } else { - if(pjet_.size() < partonsToMatch_.size()-1) return true; - } - // Otherwise, in the case of jet production allow the lowest - // contributing multiplicity process (just at NLO), namely, - // dijet production, to give rise to 1-jet and even 0-jet - // events, since these can contribute to, for example, the - // inclusive jet cross section i.e. in this case the rejection - // is only applied in the case of the next-to-lowest multiplicity - // processes (>2 parton events at LO and >3 parton events at NLO). - } else { - if(!inputIsNLO_) { - // KH - March 5th - // Removed the following line giving special treatment - // also to the LO events, to maintain consistency with - // the fortran algorithm, at least for now. So now jet - // production at LO is being treated the same as all - // other processes. - // if(partonsToMatch_.size()==2 && pjet_.size()<2) return false; - if(pjet_.size() < partonsToMatch_.size()) return true; - } else { - if(partonsToMatch_.size()<=3 && pjet_.size()<2) return false; - if(pjet_.size() < partonsToMatch_.size()-1) return true; - } - } - - // Sort partonsToMatch_ from high to low pT. - sort(partonsToMatch_.begin(),partonsToMatch_.end(),pTsortFunction); - - // Match light progenitors to jets. - vector jetToPartonMap(pjet_.size(),-999); - Energy etmin(777e100*GeV); - - // If the input is NLO events then don't do any jet-parton matching! - if(!inputIsNLO_) { - - // For each parton, starting with the hardest one ... - for(unsigned int ixx=0; ixx=0) { - jetToPartonMap[jetIndexForDRmin]=ixx; - if(ixx==0||etjet_[jetIndexForDRmin]partonsToMatch_.size() - && !inputIsNLO_) return true; - if(inputIsNLO_) { - if(!highestNLOMultiplicity_) { - if(pjet_.size()>partonsToMatch_.size()-1) return true; - } else { - if(!highestMultiplicity_&&pjet_.size()>partonsToMatch_.size()) - return true; - } - } - - // Veto events where matched jets are softer than non-matched ones, - // in the inclusive (highestMultiplicity_ = true) mode, unless we - // are dealing with NLO input events. - if(highestMultiplicity_ && !inputIsNLO_ ) { - for(unsigned int ixx=0; ixxid())==4||abs(partonsToMatch_[jxx]->id())==5)) continue; - if(partonJetDeltaR(partonsToMatch_[jxx],pjet_[ixx])etmin) return true; - } - } - - } - - // Otherwise we accept the event ... - return false; - -} - -/* Function that returns the R distance - between a particle and a jet. */ -double AlpGenHandler::partonJetDeltaR(ThePEG::tPPtr partonptr, LorentzMomentum jetmom) const { - LorentzMomentum partonmom(partonptr->momentum()); - // Calculate DY, DPhi and then DR - double DY(partonmom.eta()-jetmom.eta()); - double DPhi(partonmom.phi()-jetmom.phi()); - if(DPhi>M_PI) DPhi=2*M_PI-DPhi; - double DR(sqrt(sqr(DY)+sqr(DPhi))); - return DR; -} - - -// Initialize calorimeter for calsim_m and getjet_m. Note that -// because initialization is separte calsim_m can be called more -// than once to simulate pileup of several events. -void AlpGenHandler::calini_m() const { - - // Making sure arrays are clear before filling; - cphcal_.clear(); sphcal_.clear(); - cthcal_.clear(); sthcal_.clear(); - - // Fill array holding phi values of calorimeter cell centres. - double deltaPhi(2*M_PI/ncphi_); - for(unsigned int iphi=1; iphi<=ncphi_; iphi++) { - double phi(deltaPhi*(iphi-0.5)); // Goes phi~=0 to phi~=2*pi (iphi=0--->ncphi). - cphcal_.push_back(cos(phi)); // ==> goes from +1 ---> +1 (iphi=0--->ncphi). - sphcal_.push_back(sin(phi)); // ==> goes 0 -> 1 -> 0 -> -1 -> 0 (iphi=0--->ncphi). - } - - // Fill array holding theta values of calorimeter cell centres in Y. - double deltaY((ycmax_-ycmin_)/double(ncy_)); - for(unsigned int iy=1; iy<=ncy_; iy++) { - double Y(deltaY*(iy-0.5)+ycmin_); - double th(2*atan(exp(-Y))); // Goes bwds th~=pi to fwds th~=0 (iy=0--->ncy). - cthcal_.push_back(cos(th)); // ==> goes from -1 ---> +1 (iy=0--->ncy). - sthcal_.push_back(sin(th)); // ==> goes from 0 ---> +1 ---> 0 (iy=0--->ncy). - } - return; -} - -// Get FastJets -void AlpGenHandler::getFastJets(double rjet, Energy ejcut, double etajcut) const { - - vector particlesToCluster; - for(unsigned int ipar=0; iparmomentum().eta()); - if(y>=ycmin_&&y<=ycmax_) { - int absId(abs(particlesToCluster_[ipar]->id())); - // If it's not a lepton / top / photon it may go in the jet finder. - if(!(absId>=11&&absId<=16) && absId!=6 && absId!=22) { - // input particles into fastjet pseudojet - fastjet::PseudoJet p(particlesToCluster_[ipar]->momentum().x()/GeV, - particlesToCluster_[ipar]->momentum().y()/GeV, - particlesToCluster_[ipar]->momentum().z()/GeV, - particlesToCluster_[ipar]->momentum().e()/GeV); - p.set_user_index(ipar); - particlesToCluster.push_back(p); - } - } - } - - fastjet::RecombinationScheme recombinationScheme = fastjet::E_scheme; - fastjet::Strategy strategy = fastjet::Best; - double R(rjet); - fastjet::JetDefinition theJetDefinition; - switch (jetAlgorithm_) { - case -1: theJetDefinition=fastjet::JetDefinition(fastjet::antikt_algorithm, - R, - recombinationScheme, - strategy); break; - case 0: theJetDefinition=fastjet::JetDefinition(fastjet::cambridge_algorithm, - R, - recombinationScheme, - strategy); break; - case 1: theJetDefinition=fastjet::JetDefinition(fastjet::kt_algorithm, - R, - recombinationScheme, - strategy); break; - default: theJetDefinition=fastjet::JetDefinition(fastjet::cambridge_algorithm, - R, - recombinationScheme, - strategy); break; - } - fastjet::ClusterSequence fastjetEvent(particlesToCluster,theJetDefinition); - vector inclusiveJets = fastjetEvent.inclusive_jets(); - inclusiveJets = fastjet::sorted_by_pt(inclusiveJets); - - // Fill the array of jet momenta for the rest of the veto procedure. - pjet_.clear(); - pjet_.resize(inclusiveJets.size()); - etjet_.clear(); - etjet_.resize(inclusiveJets.size()); - for(unsigned int ffj=0; ffjetajcut) { - pjet_.erase(pjet_.begin()+fj); - etjet_.erase(etjet_.begin()+fj); - fj--; - } - - // Sort jets from high to low ET. - vector > etjet_pjet; - for(unsigned int ixx=0; ixxmomentum().eta()); - if(y>=ycmin_&&y<=ycmax_) { - int absId(abs(particlesToCluster_[ipar]->id())); - // If it's not a lepton / top / photon it goes in the calorimeter. - if(!(absId>=11&&absId<=16) && absId!=6 && absId!=22) { - double phi(atan2(particlesToCluster_[ipar]->momentum().y()/GeV, - particlesToCluster_[ipar]->momentum().x()/GeV)); - if(phi<0) phi+=2*M_PI; - unsigned int iy(int((y-ycmin_)/dely_)); - unsigned int iphi(int(phi/delphi_)); - et_[iy][iphi]+=particlesToCluster_[ipar]->momentum().e()*sthcal_[iy]; - } - } - } - return; -} - -// Find highest remaining cell > etstop and sum surrounding cells -// with -- delta(y)^2+delta(phi)^2 < Rjet^2 , ET>eccut. Keep sets -// with ET>ejcut and abs(eta)=etstop) { - - // Find the cell with the highest ET from - // those not already assigned to a jet. - etmax=0*GeV; - int iymx(0), iphimx(0); - for(unsigned int iphi=0; iphietmax&&jetIdx_[iy][iphi]<0) { - etmax = et_[iy][iphi]; - iymx = iy; - iphimx = iphi; - } - - // If the remaining cell with the highest ET has ET < etstop, stop. - if(etmax(ncy_*ncphi_)) { - cout << "AlpGenHandler::getjet_m() - Fatal error." << endl; - cout << "We found " << ipass << " calo cells with the highest ET so" - << "far\nbut the calorimeter only has " << ncy_*ncphi_ << " " - << "cells in it!" << endl; - exit(10); - } - - // Add a jet vector (may get deleted if jet fails ET / eta cuts). - etjet_.push_back(0*GeV); - pjet_.push_back(Lorentz5Momentum(0.*GeV,0.*GeV,0.*GeV,0.*GeV,0.*GeV)); - - // Loop over all calo cells in range iphimx +/- nphi1 (inclusive) - // wrapping round in azimuth if required. - for(unsigned int iphi1=0; iphi1<=2*nphi1; iphi1++) { - int iphix(iphimx-nphi1+iphi1); - if(iphix<0) iphix += ncphi_; - if(iphix>=int(ncphi_)) iphix -= ncphi_; - // Loop over all calo cells in range iymx +/- ny1 (inclusive). - for(unsigned int iy1=0; iy1<=2*ny1; iy1++) { - int iyx(iymx-ny1+iy1); - // If the cell is outside the calorimeter OR if it was already - // associated to a jet then skip to the next loop. - if(iyx>=0&&iyx=eccut) { - Energy ECell(et_[iyx][iphix]/sthcal_[iyx]); - pjet_.back()+=LorentzMomentum(ECell*sthcal_[iyx]*cphcal_[iphix], // px - ECell*sthcal_[iyx]*sphcal_[iphix], // py - ECell*cthcal_[iyx],ECell); // pz, E. - // N.B. This is the same reln as in ThePEG between phi and x,y. - etjet_.back()+=et_[iyx][iphix]; - jetIdx_[iyx][iphix] = pjet_.size()-1; // Identify cell with this jet. - } - } - } - } - - // Compute the current jet's mass. - pjet_.back().rescaleMass(); - - // Throw the jet away if it's ET is less than ejcut. - if(etjet_.back()etajcut) { - pjet_.pop_back(); - etjet_.pop_back(); - } - - } - - // Sort jets from high to low ET. - vector > etjet_pjet; - for(unsigned int ixx=0; ixxchildren()); - for (unsigned int ixx=0; ixxchildren().size()==0) - tmpList_.push_back(theChildren[ixx]); - else - getDescendents(theChildren[ixx]); - return; -} -void AlpGenHandler::caldel_m() const { - - - preshowerFSPsToDelete_.clear(); - showeredFSPsToDelete_.clear(); - - for(unsigned int ixx=0; ixxparents()[0]->id())==23|| - abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxid())==ihvy_&&ixx<2) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxparents()[0]->id())==23|| - abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxparents()[0]->id())==23|| - abs(preshowerFSPs_[ixx]->parents()[0]->id())==24|| - abs(preshowerFSPs_[ixx]->id())==22|| - abs(preshowerFSPs_[ixx]->id())==25) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxparents()[0]->id())==6|| - abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - } - if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==6) { - getDescendents(preshowerFSPs_[ixx]->parents()[0]); - for(unsigned int jxx=0; jxxid())==ihvy_&&ixx<2) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxid())==4&&ixx<1)|| - abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxid())==22) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxid())==25) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxid())==22|| - abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxid())==22|| - (abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx==(preshowerFSPs_.size()-(2+nph_+1)))|| - (abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx==(preshowerFSPs_.size()-(2+nph_+2)))|| - abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxid())==22|| - abs(preshowerFSPs_[ixx]->parents()[0]->id())==6|| - abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxid())==22|| - (abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx<2)) { - preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]); - getDescendents(preshowerFSPs_[ixx]); - for(unsigned int jxx=0; jxxparticlesToCluster_.size()) { - throw Exception() - << "AlpGenHandler::caldel_m() - ERROR!\n" - << "No. of ME level partons to be matched to jets = " - << partonsToMatch_.size() << "\n" - << "No. of showered particles to build jets from = " - << particlesToCluster_.size() << "\n" - << "There should be at least as many partons to\n" - << "cluster as there are partons to match to.\n" - << Exception::eventerror; - } - - - // Acid test. - unsigned int tmpUnsignedInt(njets_); - if(!inputIsNLO_&&partonsToMatch_.size()!=tmpUnsignedInt) { - for(unsigned int ixx=0; ixxid())>=6&& - abs(partonsToMatch_[ixx]->id())!=21) - throw Exception() - << "AlpGenHandler::caldel_m() - ERROR!\n" - << "Found a parton to match to which is not a quark or gluon!" - << *partonsToMatch_[ixx] << "\n" - << Exception::eventerror; - } - throw Exception() - << "AlpGenHandler::caldel_m() - ERROR!\n" - << "No. of ME level partons to be matched to jets = " - << partonsToMatch_.size() << "\n" - << "No. of light jets (njets) in AlpGen process = " - << njets_ << "\n" - << "These should be equal." << "\n" - << Exception::eventerror; - } - - return; -} - -// This looks for all descendents of a top up to but not including -// the W and b children. -void AlpGenHandler::getTopRadiation(PPtr theParticle) const { - ParticleVector theChildren(theParticle->children()); - for (unsigned int ixx=0; ixxchildren().size()==0) - tmpList_.push_back(theChildren[ixx]); - else if(abs(theChildren[ixx]->id())==5||abs(theChildren[ixx]->id())==24) - return; - else - getTopRadiation(theChildren[ixx]); - return; -} -void AlpGenHandler::caldel_hvq() const { - - // Fill partonsToMatch_ with only those pre-shower partons intended to - // be used in heavy-quark-jet matching and fill particlesToCluster_ using - // only those final state particles (post-shower) which are supposed - // in the heavy-quark-jet clustering used to do merging. To begin with - // these are made from the corresponding sets of particles that were - // omitted from the initial jet-parton matching run. - partonsToMatch_ = preshowerFSPsToDelete_; - particlesToCluster_.resize(showeredFSPsToDelete_.size()); - for(unsigned int ixx=0; ixxid())<4||abs(partonsToMatch_[ixx]->id())>6) { - preshowerFSPsToDelete_.push_back(partonsToMatch_[ixx]); - tmpList_.clear(); - getDescendents(partonsToMatch_[ixx]); - for(unsigned int jxx=0; jxxid())==5&& - partonsToMatch_[ixx]->parents().size()>0&& - abs(partonsToMatch_[ixx]->parents()[0]->id())==6) { - preshowerFSPsToDelete_.push_back(partonsToMatch_[ixx]); - tmpList_.clear(); - getDescendents(partonsToMatch_[ixx]); - for(unsigned int jxx=0; jxxparents().size()>0&& - (abs(partonsToMatch_[ixx]->parents()[0]->id())==23|| - abs(partonsToMatch_[ixx]->parents()[0]->id())==24|| - abs(partonsToMatch_[ixx]->parents()[0]->id())==25)) { - preshowerFSPsToDelete_.push_back(partonsToMatch_[ixx]); - tmpList_.clear(); - getDescendents(partonsToMatch_[ixx]); - for(unsigned int jxx=0; jxxsubProcess()->intermediates()); - for(unsigned int ixx=0; ixxid())==6) { - partonsToMatch_.push_back(intermediates[ixx]); - tmpList_.clear(); - getTopRadiation(partonsToMatch_.back()); - for(unsigned int jxx=0; jxxid())>=4&&abs(partonsToMatch_[ixx]->id())<=6) { - theProgenitor = partonsToMatch_[ixx]; - // Follow the heavy quark line down to where it stops branching. - while(theProgenitor->children().size()>0) { - theLastProgenitor = theProgenitor; - for(unsigned int jxx=0; jxxchildren().size(); jxx++) { - if(theProgenitor->children()[jxx]->id()==theProgenitor->id()) - theProgenitor=theProgenitor->children()[jxx]; - } - // If the progenitor had children but none of them had - // the same particle id as it, then it must have undergone - // a decay rather than a branching, i.e. it is the end of - // the evolution line, so, - if(theProgenitor==theLastProgenitor) break; - } - evolvedHeavyQuarks.push_back(theProgenitor); - } - } - // Now delete the evolved heavy quark from the particlesToCluster. - for(unsigned int ixx=0; ixxsubProcess()->incoming(); - - // LH file final-state partICLEs: - preshowerFSPs_ = lastXCombPtr()->subProcess()->outgoing(); - - return; -} - -void AlpGenHandler::getShoweredParticles() const { - // Post-shower initial-state hadrons: - showeredISHs_ = eventHandler()->currentEvent()->incoming(); - - // Post-shower initial-state partons: - for(unsigned int ixx=0; ixx<(showeredISHs_.first)->children().size(); ixx++) - if(((showeredISHs_.first)->children()[ixx]->id())<6|| - ((showeredISHs_.first)->children()[ixx]->id())==21) - showeredISPs_.first=(showeredISHs_.first)->children()[ixx]; - for(unsigned int ixx=0; ixx<(showeredISHs_.second)->children().size(); ixx++) - if(((showeredISHs_.second)->children()[ixx]->id())<6|| - ((showeredISHs_.second)->children()[ixx]->id())==21) - showeredISPs_.second=(showeredISHs_.second)->children()[ixx]; - - // Post-shower final-state partICLEs plus remnants (to be removed later): - showeredFSPs_ = eventHandler()->currentEvent()->getFinalState(); - - // Post-shower final-state remnants: - for(unsigned int ixx=0; ixxPDGName()=="Rem:p+"|| - showeredFSPs_[ixx]->PDGName()=="Rem:pbar-") { - if(showeredFSPs_[ixx]->parents()[0]->parents()[0]== - showeredISHs_.first) - showeredRems_.first=showeredFSPs_[ixx]; - else if(showeredFSPs_[ixx]->parents()[0]->parents()[0]== - showeredISHs_.second) - showeredRems_.second=showeredFSPs_[ixx]; - } - } - - // Now delete found remnants from the showeredFSPs vector for consistency. - for(unsigned int ixx=0; ixxPDGName()=="Rem:p+") - showeredFSPs_.erase(showeredFSPs_.begin()+ixx); - for(unsigned int ixx=0; ixxPDGName()=="Rem:pbar-") - showeredFSPs_.erase(showeredFSPs_.begin()+ixx); - sort(showeredFSPs_.begin(),showeredFSPs_.end(),recordEntry); - - return; -} - -void AlpGenHandler::doSanityChecks(int debugLevel) const { - - // When checking momentum conservation in the form - // p_in - p_out, any momentum component bigger / less - // than + / - epsilon will result in the p_in - p_out - // vector being flagged as "non-null", triggering a - // warning that momentum conservation is violated. - Energy epsilon(0.5*GeV); - if(debugLevel>=5) epsilon=1e-9*GeV; - - // Print out what was found for the incoming and outgoing - // partons in the lastXCombPtr regardless. - if(debugLevel>=5) { - cout << "\n\n\n\n"; - cout << "****************************************************" << endl; - cout << " The following are the hard subprocess momenta from " << "\n" - << " lastXCombPtr and should be basically identical to " << "\n" - << " the input LH file momenta." << "\n\n"; - cout << " Incoming particles:" << "\n" - << *(preshowerISPs_.first) << "\n" - << *(preshowerISPs_.second) << endl; - cout << " Outgoing particles:" << endl; - for(unsigned int ixx=0; ixx=5) { - cout << "\n\n"; - cout << "****************************************************" << endl; - cout << " The following are the particles left at the end of" << "\n" - << " the showering step." << "\n\n"; - cout << " Incoming hadrons:" << "\n" - << *(showeredISHs_.first) << "\n" - << *(showeredISHs_.second) << endl; - cout << " Incoming partons:" << "\n" - << *(showeredISPs_.first) << "\n" - << *(showeredISPs_.second) << endl; - cout << " Outgoing partons:" << endl; - for(unsigned int ixx=0; ixx=4) { - Lorentz5Momentum tmpMom; - tmpMom += showeredISPs_.first->momentum(); - tmpMom += showeredISPs_.second->momentum(); - for(unsigned int ixx=0; ixxmomentum(); - if(!isMomLessThanEpsilon(tmpMom,epsilon)) - cout << "Total parton mom.in - total parton mom.out = " - << tmpMom/GeV << endl; - tmpMom = showeredISHs_.first->momentum() - - showeredRems_.first->momentum() -showeredISPs_.first->momentum(); - if(!isMomLessThanEpsilon(tmpMom,epsilon)) - cout << "First p_hadron-p_remnant-p_incoming " << tmpMom/GeV << endl; - tmpMom = showeredISHs_.second->momentum() - - showeredRems_.second->momentum()-showeredISPs_.second->momentum(); - if(!isMomLessThanEpsilon(tmpMom,epsilon)) - cout << "Second p_hadron-p_remnant-p_incoming " << tmpMom/GeV << endl; - } - - // Check if what we found to be the remnant is consistent with - // what we identified as the parent incoming hadron i.e. p+ - // goes with Rem:p+ and pbar- goes with Rem:pbar-. - if(debugLevel>=0) { - string tmpString; - tmpString=showeredRems_.first->PDGName(); - tmpString=tmpString.substr(tmpString.find_first_of(":")+1, - string::npos); - if(showeredISHs_.first->PDGName()!=tmpString) { - cout << "AlpGenHandler::showerHardProcessVeto" << "\n" - << "Fatal error in pairing of remnant and parent hadron." << "\n" - << "Remnant = " << *(showeredRems_.first) << "\n" - << "Parent hadron = " << *(showeredISHs_.first) - << endl; - cout << showeredISHs_.first->PDGName() << endl; - cout << tmpString << endl; - } - tmpString=showeredRems_.second->PDGName(); - tmpString=tmpString.substr(tmpString.find_first_of(":")+1, - string::npos); - if(showeredISHs_.second->PDGName()!=tmpString) { - cout << "AlpGenHandler::showerHardProcessVeto" << "\n" - << "Fatal error in pairing of remnant and parent hadron." << "\n" - << "Remnant = " << *(showeredRems_.second) << "\n" - << "Parent hadron = " << *(showeredISHs_.second) - << endl; - cout << showeredISHs_.second->PDGName() << endl; - cout << tmpString << endl; - } - } - return; -} - -void AlpGenHandler::printMomVec(vector momVec) { - cout << "\n\n"; - - // Label columns. - printf("%5s %9s %9s %9s %9s %9s %9s %9s %9s %9s\n", - "jet #", - "px","py","pz","E", - "eta","phi","pt","et","mass"); - - // Print out the details for each jet - for (unsigned int ixx=0; ixx ETSTOP and sum surrounding cells with -- - * DELTA(Y)**2+DELTA(PHI)**2 < RJET**2 , ET>ECCUT. - * Keep sets with ET>EJCUT and ABS(ETA)currentEvent()->... - * filling the particle pairs showeredISHs_, showeredISPs_, - * showeredRems_ and the particle pointer vector showeredFSPs_. - */ - void getShoweredParticles() const; - - /** - * Allows printing of debug output and sanity checks like - * total momentum consrvation to be carried out. - * debugLevel = -1, 0, ...5 - * = no debugging, minimal debugging, ... verbose. - */ - void doSanityChecks(int debugLevel) const; - - /** - * Given a pointer to a particle this finds all its final state - * descendents. - */ - void getDescendents(PPtr theParticle) const; - - /** - * Accumulates all descendents of tops down to the b and W - * but not including them. - */ - void getTopRadiation(PPtr theParticle) const; - - /** - * Sorts a given vector of particles by descending pT or ETJET - */ - - ParticleVector pTsort(ParticleVector unsortedVec); - pair< vector, vector > ETsort(vector unsortedetjet, vector unsortedVec); - - /* - * A function that prints a vector of Lorentz5Momenta in a fancy way - */ - void printMomVec(vector momVec); - - - /* - * A probability function for varying etclus_ about the mean value - */ - Energy etclusran_(double petc) const; - -private: - - /** - * The static object used to initialize the description of this class. - * Indicates that this is a concrete class with persistent data. - */ - static ClassDescription initAlpGenHandler; - - /** - * The assignment operator is private and must never be called. - * In fact, it should not even be implemented. - */ - AlpGenHandler & operator=(const AlpGenHandler &) = delete; - -private: - - /** - * Initial-state incoming partons prior to showering - * (i.e. from lastXCombPtr). - */ - mutable PPair preshowerISPs_; - - /** - * Final-state outgoing partICLEs prior to showering - * (i.e. from lastXCombPtr). - */ - mutable ParticleVector preshowerFSPs_; - - /** - * Final-state outgoing partICLEs prior to showering _to_be_removed_ - * from preShowerFSPs_ prior to the light-parton-light-jet matching - * step. This same list is the starting point for determining - * partonsToMatch_ for the case of merging in heavy quark production. - */ - mutable ParticleVector preshowerFSPsToDelete_; - - /** - * Initial-state incoming hadrons after shower of hard process - * (eventHandler()->currentEvent()->incoming()). - */ - mutable PPair showeredISHs_; - - /** - * Initial-state incoming partons after shower of hard process - * (look for partonic children of showeredISHs_). - */ - mutable PPair showeredISPs_; - - /** - * Final-state outgoing partICLEs after shower of hard process - * (eventHandler()->currentEvent()->getFinalState()). - */ - mutable tPVector showeredFSPs_; - - /** - * Final-state outgoing partICLEs after shower of hard process - * _to_be_removed_ from showeredFSPs_ prior to the - * light-parton-light-jet matching step. This same list is the - * starting point for determining particlesToCluster_ for the - * case of merging in heavy quark production. - */ - mutable ParticleVector showeredFSPsToDelete_; - - /** - * ONLY the final-state partons from preshowerFSPs_ that are - * supposed to enter the jet-parton matching. - */ - mutable ParticleVector partonsToMatch_; - - /* - * The shower progenitors - */ - - mutable PPtr theProgenitor; - mutable PPtr theLastProgenitor; - - /** - * ONLY the final-state particles from showeredFSPs_ (and maybe - * also showeredRems_) that are supposed to go for jet clustering. - */ - mutable tPVector particlesToCluster_; - - /** - * Final-state remnants after shower of hard process - * (look for remnants initially in showeredFSPs_). - */ - mutable PPair showeredRems_; - - /** - * Pointer to the object calculating the strong coupling - */ - ShowerAlphaPtr alphaS_; - - /** - * Information extracted from the XComb object - */ - //@{ - /** - * The fixed factorization scale used in the MEs. - */ - Energy pdfScale_; - - /** - * Centre of mass energy - */ - Energy2 sHat_; - - /** - * Constant alphaS used to generate LH events - if not already - * using CKKW scale (ickkw = 1 in AlpGen for example). - */ - double alphaSME_; - //@} - - /* - * Number of rapidity segments of the calorimeter. - */ - unsigned int ncy_; - - /* - * Number of phi segments of the calorimeter. - */ - unsigned int ncphi_; - - /* - * Heavy flavour in WQQ,ZQQ,2Q etc (4=c, 5=b, 6=t). - */ - int ihvy_; - - /* - * Number of photons in the AlpGen process. - */ - int nph_; - - /* - * Number of higgses in the AlpGen process. - */ - int nh_; - - /* - * Jet ET cut to apply in jet clustering (in merging). - */ - mutable Energy etclus_; - - /* - * Mean Jet ET cut to apply in jet clustering (in merging). - */ - Energy etclusmean_; - - /* - * maximum deviation from mean Jet ET cut to apply in jet clustering (in merging). - */ - Energy epsetclus_; - - - - /* - * Cone size used in jet clustering (in merging). - */ - double rclus_; - - /* - * Max |eta| for jets in clustering (in merging). - */ - double etaclmax_; - - /* - * Default 1.5 factor used to decide if a jet matches a parton - * in merging: if DR(parton,jet) ncphi). - * ==> Cosine goes from +1 ---> +1 (index = 0 ---> ncphi). - */ - mutable vector cphcal_; - - /* - * Sine of phi values of calorimeter cell centres. - * Goes phi~=0 to phi~=2*pi (index = 0 ---> ncphi). - * ==> Sine goes 0 -> 1 -> 0 -> -1 -> 0 (index = 0 ---> ncphi). - */ - mutable vector sphcal_; - - /* - * Cosine of theta values of calorimeter cell centres in Y. - * Goes bwds th~=pi to fwds th~=0 (index = 0 ---> ncy). - * ==> Cosine goes from -1 ---> +1 (index = 0 ---> ncy). - */ - mutable vector cthcal_; - - /* - * Sine of theta values of calorimeter cell centres in Y. - * Goes bwds th~=pi to fwds th~=0 (index = 0 ---> ncy). - * ==> Sine goes from 0 ---> +1 ---> 0 (index = 0 ---> ncy). - */ - mutable vector sthcal_; - - /* - * Transverse energy deposit in a given calorimeter cell. - * First array index corresponds to rapidity index of cell, - * second array index corresponds to phi cell index. - */ - mutable vector > et_; - - /* - * For a given calorimeter cell this holds the index of the jet - * that the cell was clustered into. - */ - mutable vector > jetIdx_; - - /* - * Vector holding the Lorentz 5 momenta of each jet. - */ - mutable vector pjet_; - - /* - * Vector holding the list of FS particles resulting from - * the particle input to getDescendents. - */ - mutable ParticleVector tmpList_; - - /* - * Variables for the C++ translation of the calini_m(), calsim_m(), - * getjet_m(...) and caldel_m() functions - */ - mutable vector etjet_; - mutable double dely_, delphi_; - -}; - -} - -#include "ThePEG/Utilities/ClassTraits.h" - -namespace ThePEG { - -/** @cond TRAITSPECIALIZATIONS */ - -/** This template specialization informs ThePEG about the - * base classes of AlpGenHandler. */ -template <> -struct BaseClassTrait { - /** Typedef of the first base class of AlpGenHandler. */ - typedef Herwig::QTildeShowerHandler NthBase; -}; - -/** This template specialization informs ThePEG about the name of - * the AlpGenHandler class and the shared object where it is defined. */ -template <> -struct ClassTraits - : public ClassTraitsBase { - /** Return a platform-independent class name */ - static string className() { return "Herwig::AlpGenHandler"; } - /** - * The name of a file containing the dynamic library where the class - * AlpGenHandler is implemented. It may also include several, space-separated, - * libraries if the class AlpGenHandler depends on other classes (base classes - * excepted). In this case the listed libraries will be dynamically - * linked in the order they are specified. - */ - static string library() { return "AlpGenHandler.so"; } -}; - -/** @endcond */ - -} - -#endif /* HERWIG_AlpGenHandler_H */ diff --git a/Contrib/AlpGen/AlpGenToLH.cc b/Contrib/AlpGen/AlpGenToLH.cc deleted file mode 100644 --- a/Contrib/AlpGen/AlpGenToLH.cc +++ /dev/null @@ -1,1548 +0,0 @@ -//A few standard headers -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace std; - -double sqr(double x); - -inline int nInt(double x) { - int theCeiling=int(ceil(x)); - int theFloor=int(floor(x)); - if((theCeiling-theFloor!=1&&theCeiling!=theFloor) - ||theCeiling-x>1.0||x-theFloor>1.0||xtheCeiling) { - cout << "nInt:\n" - << "Fatal double to integer conversion error.\n" - << "input double = " << x << "\n" - << "integer ceiling = " << theCeiling << "\n" - << "integer floor = " << theFloor << "\n" - << "Quitting ..."; - exit(1); - } - return (theCeiling-x) < (x-theFloor) ? theCeiling : theFloor; -} - -int ndnsToLHAPDF(int ndns); - -string ndnsToLHAPDF_str(int ndns); - - -double parstrToparval(string varName, - vector * parstrPtr, - vector * parvalPtr); - -void doIndividualHardProcessAssignments(int ihrd , double * nup, - vector * idup , vector * istup, - vector * mothup1, vector * mothup2, - vector * icolup1, vector * icolup2, - vector > * pup, - vector masses , int itopprc); - -//*****************************// -//*****************************// -//*****************************// -// // -// ATTENTION!!!! // -// ATTENTION!!!! // -// ATTENTION!!!! // -// // -// Remember to uncomment all // -// of the currently commented // -// commands outputting to // -// the HW++ input file for // -// the variables not yet // -// declared and interfaced // -// in the AlpGenHandler!! // -// These should only be the // -// ones commented out by _4_ // -// consecutive forward slashes // -// '////' in the function // -// writeHWPPinFile(...) // -// // -// These are only commented // -// as the AlpGenHandler isn't // -// ready for them yet! // -// // -//*****************************// -//*****************************// -//*****************************// - -void writeHWPPinFile(string prefix, int ihrd, int unwev, - int lhapdf, string lhapdfstr, int idbmup0, int idbmup1, int idwtup, - double aqcdup, int nloop, - vector * massesPtr, - vector * parvalPtr, - vector * parstrPtr); - -string trim(string theString); - -int main(int argc, char *argv[]) { - - bool usePowhegBoxConventions(true); // Control - int debugging(3); // To validate - - cout << "\n"; - cout << "------------------------------------------------------------------\n"; - cout << " AlpGenToLH: Convert Alpgen Les Houches to Herwig Les Houches \n"; - cout << " v2.0-beta \n"; - cout << "------------------------------------------------------------------\n"; - cout << "\n"; - char* prefix; - - - if(argv[1]) { prefix = argv[1]; } else { - cout << "Use: ./AlpGenToLH [input string] [number of events (optional)]\n"; exit(1); - cout << "Note: As of version 2, the .stat and _unw.par files are " - << "required to read the generation parameters.\n"; - } - int maxevents(0); int eventcount(0); - if(argc>2) { maxevents=(atoi(argv[2])); } - - string lheFilename = string(prefix) + string(".lhe"); - string unwFilename = string(prefix) + string(".unw"); - string unwparFilename = string(prefix) + string("_unw.par"); - string statFilename = string(prefix) + string(".stat"); - - cout << "Opening input files ...\n"; - cout << "-----------------------\n"; - cout << "Unweighted events in " << unwFilename << ".\n" - << "Generation settings in " - << unwparFilename << " and " << statFilename << ".\n\n"; - - ifstream unwStream; - ifstream unwparStream; - ifstream statStream; - ofstream lheStream; - unwStream.open(unwFilename.c_str()); - unwparStream.open(unwparFilename.c_str()); - statStream.open(statFilename.c_str()); - lheStream.open(lheFilename.c_str()); - if(!unwStream) { - cerr << "error: Failed to open input file " << unwFilename << "\n"; - exit(1); - } - if(!unwparStream) { - cerr << "error: Failed to open input file " << unwparFilename << "\n"; - exit(1); - } - if(!statStream) { - cerr << "error: Failed to open input file " << statFilename << "\n"; - exit(1); - } - - // ******************************************************************** // - // Dump the AlpGen *_unw.par file into the LH header (it's not so big). // - // ******************************************************************** // - - string tmpString; - lheStream << "\n"; - lheStream << "\n"; - unwparStream.close(); - - // ***************************************** // - // Read in all relevant info from *_unw.par. // - // ***************************************** // - - int ihrd; // AlpGen convention hard process code. - double mc,mb,mt,mw,mz,mh; // C, B, Top, W, Z & Higgs mass from *_unw.par - double avgwgt,errwgt; // Average weight and its error. - int unwev; // Number of unweighted events. - double totlum; // Effective luminosity. - vector parval(200,-999.0); // AlpGen parameters. - vector parstr(200,"----"); // AlpGen parameter (variable) names. - vector alpgenParticleMasses; - - unwparStream.open(unwparFilename.c_str()); - while(unwparStream) { - getline(unwparStream,tmpString); - if(tmpString.find("hard process code") != string::npos) { - tmpString=tmpString.substr(0,tmpString.find("!")); - tmpString=trim(tmpString); - ihrd=atoi(tmpString.c_str()); - } - if(tmpString.find("mc,mb,mt,mw,mz,mh") != string::npos) { - tmpString=trim(tmpString.substr(0,tmpString.find("!"))); - mc=atof((tmpString.substr(0,tmpString.find_first_of(" "))).c_str()); - alpgenParticleMasses.push_back(mc); - tmpString=trim(tmpString.substr(tmpString.find_first_of(" "))); - mb=atof((tmpString.substr(0,tmpString.find_first_of(" "))).c_str()); - alpgenParticleMasses.push_back(mb); - tmpString=trim(tmpString.substr(tmpString.find_first_of(" "))); - mt=atof((tmpString.substr(0,tmpString.find_first_of(" "))).c_str()); - alpgenParticleMasses.push_back(mt); - tmpString=trim(tmpString.substr(tmpString.find_first_of(" "))); - mw=atof((tmpString.substr(0,tmpString.find_first_of(" "))).c_str()); - alpgenParticleMasses.push_back(mw); - tmpString=trim(tmpString.substr(tmpString.find_first_of(" "))); - mz=atof((tmpString.substr(0,tmpString.find_first_of(" "))).c_str()); - alpgenParticleMasses.push_back(mz); - tmpString=trim(tmpString.substr(tmpString.find_first_of(" "))); - mh=atof((tmpString.substr(0,tmpString.find_first_of(" "))).c_str()); - alpgenParticleMasses.push_back(mh); - } - if(tmpString.find("Crosssection +- error (pb)") != string::npos) { - tmpString=trim(tmpString.substr(0,tmpString.find("!"))); - avgwgt=atof(trim(tmpString.substr(0,tmpString.find(" "))).c_str()); - errwgt=atof(trim(tmpString.substr(tmpString.find(" "))).c_str()); - } - if(tmpString.find("unwtd events, lum (pb-1)") != string::npos) { - tmpString=trim(tmpString.substr(0,tmpString.find("!"))); - unwev=atoi(trim(tmpString.substr(0,tmpString.find(" "))).c_str()); - totlum=atof(trim(tmpString.substr(tmpString.find(" "))).c_str()); - } - } - if(maxevents > unwev) { - cout << "-------------------------------\n"; - cout << "requested " << maxevents << " > " << unwev << " (contained in file), will use all events.\n"; maxevents = 0; } - - if(debugging>=4) { - cout << "\nDebugging initial reading of *_unw.par:\n"; - cout << "ihrd = " << ihrd << "\n"; - cout << "mc,mb,mt,mw,mz,mh = " - << mc << " " << mb << " " - << mt << " " << mw << " " - << mz << " " << mh << "\n"; - cout << "Cross section +/- error = " - << avgwgt << " +/- " << errwgt << "\n"; - cout << "Number of unweighted events = " << unwev << "\n"; - cout << "Effective luminosity = " << totlum << "\n"; - } - unwparStream.close(); - - unwparStream.open(unwparFilename.c_str()); - int index; - while(unwparStream) { - getline(unwparStream,tmpString); - if(tmpString.find("!")==string::npos|| - tmpString.find("hard process code")!=string::npos|| - tmpString.find("mc,mb,mt,mw,mz,mh")!=string::npos|| - tmpString.find("Crosssection +- error (pb)")!=string::npos|| - tmpString.find("unwtd events, lum (pb-1)")!=string::npos) continue; - tmpString=trim(tmpString); - if(debugging>=4) cout << "\nDebugging reading paramters in *_unw.par:\n"; - if(debugging>=4) cout << "File says: " << tmpString << "\n"; - index = atoi((tmpString.substr(0,tmpString.find_first_of(" "))).c_str()); - tmpString=trim(tmpString.substr(tmpString.find_first_of(" "))); - parval[index]=atof((tmpString.substr(0,tmpString.find_first_of(" "))).c_str()); - tmpString=trim(tmpString.substr(tmpString.find_first_of("!")+1)); - parstr[index]=tmpString; - if(debugging>=4) cout << "We say: " - << index << " " - << parval[index] << " " - << parstr[index] << "\n\n\n\n\n"; - } - unwparStream.close(); - - // Variables defined in parval array read from *_unw.par: - - // PDG codes for the beam particles. - int idbmup[2]={0,0}; - if(parstrToparval("ih1",&parstr,&parval)==1) idbmup[0] = 2212; - else if(parstrToparval("ih1",&parstr,&parval)==-1) idbmup[0] = -2212; - else idbmup[0] = 2212; - if(parstrToparval("ih2",&parstr,&parval)==1) idbmup[1] = 2212; - else if(parstrToparval("ih2",&parstr,&parval)==-1) idbmup[1] = -2212; - else idbmup[1] = 2212; - - // Energies of the beam particles --- implementation implicitly assumes - // these are equal! - double ebmup[2]={0,0}; - ebmup[0]=parstrToparval("ebeam",&parstr,&parval); - ebmup[1]=ebmup[0]; - - // LH accord pdf info variables for block: - int pdfgup[2]; - pdfgup[0]=-1; // Simply set to 1 as in POWHEG-BOX. - pdfgup[1]=-1; - int pdfsup[2]; - // LHAPDF index: (note in POWHEG-BOX it is set to just -1). - pdfsup[0]=ndnsToLHAPDF(int(parstrToparval("ndns",&parstr,&parval))); - pdfsup[1]=pdfsup[0]; - string lhastring = ndnsToLHAPDF_str(int(parstrToparval("ndns",&parstr,&parval))); - // LH accord flag defining weight scheme: - // N.B. AlpGen alpsho.f UPINIT uses idwtup = 3 (this is likely better from the - // point of view of combining events of diff multiplicity together in real life - // i.e. in ATLAS - so we should probably use it!). - int idwtup( 3); // As in POWHEG-BOX withnegweights 0 mode: unit wgts +1 only. - // int idwtup(-4); // As in POWHEG-BOX withnegweights 1 mode: +/- |xsecup| wgts only. - // Number of processes in the file (assume all one process as - // with alpsho.f UPINIT). - int nprup(1); - // Cross section, it's error, the maximum weight in the file. - double xsecup,xerrup,xmaxup; - xsecup = avgwgt; - xerrup = errwgt; - xmaxup = xsecup; - // Process id code (to be augmented by jet multiplicity - see just below). - int lprup(ihrd*100); - - // We augment the process code (for the LH file only) by the number of - // (light) jets, just in case we end up connecting many different files - // to the shower MC in parallel (otherwise it likely won't distinguish - // between X+0,1,2,3,...,n jet processes). - int njets(int(parstrToparval("njets",&parstr,&parval))); - lprup+=njets; - - // Write out some bits of info to the screen. - cout << "No. of jets: " << njets << "\n"; - cout << "Total xsec in pb (all processes): " << scientific - << xsecup << " +/- " << xerrup << "\n\n"; - - // NOW write out block: - lheStream << "\n"; - lheStream << setw(9) << idbmup[0]; - lheStream << setw(9) << idbmup[1]; - lheStream << scientific << setprecision(5) << setw(13) << ebmup[0]; - lheStream << scientific << setprecision(5) << setw(13) << ebmup[1]; - lheStream << setw(7) << pdfgup[0]; - lheStream << setw(7) << pdfgup[1]; - if(usePowhegBoxConventions) { - lheStream << setw(7) << -1; - lheStream << setw(7) << -1; - } else { - lheStream << setw(7) << pdfsup[0]; - lheStream << setw(7) << pdfsup[1]; - } - lheStream << setw(7) << idwtup; - lheStream << setw(7) << nprup << "\n"; - - lheStream << scientific << setprecision(5) << setw(13) << xsecup; - lheStream << scientific << setprecision(5) << setw(13) << xerrup; - if(usePowhegBoxConventions) - lheStream << scientific << setprecision(5) << setw(13) << 1.0; - else // else put in more info (xmaxup=xsecup). - lheStream << scientific << setprecision(5) << setw(13) << xmaxup; - lheStream << setw(7) << lprup << "\n"; - - lheStream << "\n"; - - // ************************************************* // - // All done with the init section of the LH file! // - // ************************************************* // - - // ***************************************************************** // - // To write the HW++ .in file we have everything we could possibly // - // want except maybe the QCD coupling and no. of loops for running. // - // These are the only numbers we get / use from *.stat now. // - // ***************************************************************** // - - double aqedup(-999),aqedStat(-999); - double aqcdup(-999),aqcdStat(-999); - int nloop(-999); - - // Fish around for the QCD and QED alphas in .stat. - while(statStream) { - getline(statStream,tmpString); - if(tmpString.find("as(MZ)") != string::npos) { - aqcdup=atof(trim(tmpString.substr(tmpString.find_last_of("=")+1, - tmpString.length())).c_str()); - tmpString=trim(tmpString.substr(tmpString.find_first_of("nloop="))); - tmpString=trim(tmpString.substr(tmpString.find_first_of(" "))); - tmpString=trim(tmpString.substr(0,tmpString.find_first_of("]"))); - nloop=atoi(tmpString.c_str()); - } - if(tmpString.find("aem(mZ)") != string::npos) { - tmpString=tmpString.substr(tmpString.find("aem(mZ)=")); - tmpString=tmpString.substr(tmpString.find_first_of(" ")); - tmpString=trim(tmpString); - aqedup=atof(tmpString.c_str()); - aqedup=1/aqedup; - } - } - statStream.close(); - aqedStat=aqedup; - aqcdStat=aqcdup; - - // Write out a couple more bits of info to the screen. - cout << "aqcdup [as(MZ)] from stat file: " << aqcdup << "\n"; - cout << "nloop for as from stat file : " << nloop << "\n"; - cout << "aqedup [inverse] from stat file: " << 1/aqedup << "\n"; - cout << "\n"; - - writeHWPPinFile(prefix,ihrd,unwev,pdfsup[0],lhastring, idbmup[0], idbmup[1], idwtup, - aqcdup,nloop, - &alpgenParticleMasses, - &parval,&parstr); - - // *********************************************************** // - // All done writing the HW++ .in file! // - // *********************************************************** // - - // *********************************************************** // - // Start reading AlpGen events and writing them as LH events: // - // *********************************************************** // - - int nupMax(20); - - // First line of an event contains: - double nup(0),idprup(0),xwgtup(0),scalup(0); - - // Subsequent lines for particle content contain: - vector idup,istup; - vector mothup1,mothup2; // WARNING: not implemented ... YET! - vector icolup1,icolup2; - vector vtimup,spinup; // WARNING: not implemented (but safe). - vector > pup; - pup.resize(5); - for (int ixx = 0; ixx <= 4; ++ixx) - pup[ixx].resize(nupMax); - - // Initialise pup matrix (not really necessary): - for(unsigned int jup=0; jup < nupMax; jup++) - for(unsigned int ixx=0; ixx<5; ixx++) pup[ixx][jup]=0; - - // Control reading AlpGen file: - int iup(0); - int counter(0); - bool readInWholeEvent(false),beginNewEvent(true); - double stringdoub; - - // Needed as input to doIndividualHardProcessAssignments only: - int itopprc(nInt(parstrToparval("itopprc",&parstr,&parval))); - - while(unwStream && (eventcount < maxevents || maxevents==0)) { // So long as we haven't hit the EOF do ... - - if(beginNewEvent) { - // Rest / set control variables: - beginNewEvent=false; - readInWholeEvent=false; - counter=0; - iup=0; - // Reset variables for first line of LH event: - nup=0; idprup=0; xwgtup=0; scalup=0; aqedup=0; aqcdup=0; - // Reset variables for all individual particles in LH event: - idup.clear() ; istup.clear(); - mothup1.clear(); mothup2.clear(); - icolup1.clear(); icolup2.clear(); - vtimup.clear() ; spinup.clear(); - for(unsigned int jup=0; jup < nupMax; jup++) - for(unsigned int ixx=0; ixx<5; ixx++) pup[ixx][jup]=0; - } - - // Read in next thing starting from last position in - // the file (int/real) as a double - unwStream >> stringdoub; - - // counter counts the number of numbers read-in for a given - // event. On starting to read a new event counter will be 0. - counter++; - - switch (counter) { - case 1: if(int(stringdoub)!=0) { -// if(int(stringdoub)%100==0) -// cout << "Processed " << fixed << setprecision(0) -// << stringdoub/unwev*100 -// << " % of events ..." << "\r" << flush; - } - break; - case 2: idprup = stringdoub; break; - case 3: nup = stringdoub; break; - case 4: xwgtup = stringdoub; break; - case 5: scalup = stringdoub; break; - // N.B. There are no aqedup / aqcdup variables from AlpGen. - - // Initial state particle +z direction: - case 6: idup.push_back(stringdoub) ; - istup.push_back(-1); break; - case 7: icolup1.push_back(stringdoub); - mothup1.push_back(0); break; // ATTENTION: not given by AlpGen - case 8: icolup2.push_back(stringdoub); - mothup2.push_back(0); break; // ATTENTION: not given by AlpGen - case 9: pup[2][iup]=stringdoub ; - pup[3][iup]=fabs(stringdoub) ; iup++ ; break; - - // Initial state particle -z direction: - case 10: idup.push_back(stringdoub) ; - istup.push_back(-1); break; - case 11: icolup1.push_back(stringdoub); - mothup1.push_back(0); break; // ATTENTION: not given by AlpGen - case 12: icolup2.push_back(stringdoub); - mothup2.push_back(0); break; // ATTENTION: not given by AlpGen - case 13: pup[2][iup]=stringdoub; - pup[3][iup]=fabs(stringdoub) ; iup++ ; break; - } - - if(debugging<5) idprup=lprup; - - if(counter<14) continue; - - // Final state particles: - if(counter==0+7*iup) idup.push_back(stringdoub); - // istup gets assigned later on to just -1/+1 (I.S. / F.S.). - if(counter==1+7*iup) { icolup1.push_back(stringdoub); - mothup1.push_back(1.); } // ATTENTION: not given by AlpGen - if(counter==2+7*iup) { icolup2.push_back(stringdoub); - mothup2.push_back(2.); } // ATTENTION: not given by AlpGen - if(counter==3+7*iup) pup[0][iup] = stringdoub; - if(counter==4+7*iup) pup[1][iup] = stringdoub; - if(counter==5+7*iup) pup[2][iup] = stringdoub; - if(counter==6+7*iup) pup[4][iup] = stringdoub; - if(counter==6+7*iup) istup.push_back(1.); - if(counter==6+7*iup) iup+=1; - if(counter==7*nup-1) readInWholeEvent = true; - - for(int jup = 0; jup < nup; jup++) { - pup[3][jup] = sqrt( sqr(pup[4][jup]) + sqr(pup[0][jup]) - + sqr(pup[1][jup]) + sqr(pup[2][jup]) ); - vtimup.push_back(0.); // ATTENTION: not implemented - so taking - spinup.push_back(9.); // POWHEG-BOX default values (should be v.safe). - } - - // ***************************************************************** // - // Now consider assignments specific to individual hard processes: // - // ***************************************************************** // - - if(readInWholeEvent) - doIndividualHardProcessAssignments(ihrd, &nup, &idup, &istup, - &mothup1 , &mothup2, - &icolup1 , &icolup2, - &pup , - alpgenParticleMasses, itopprc); - - if(readInWholeEvent) { - lheStream << "\n"; - if(debugging>=5) { - lheStream << nup << "\t" << idprup << "\t" << xwgtup << "\t" - << scalup << "\t" << "0.007297352" << "\t" << "0.118\n"; - } else { - // Bit about signing of xwgtup here is redundant as - // AlpGen only gives +ve weight events ... - double signOfXwgtup = xwgtup >= 0 ? 1 : -1; - xwgtup = idwtup==3 ? 1 : xsecup*signOfXwgtup; - // N.B. There are no aqedup / aqcdup variables from AlpGen - // events only the stat file has related information. - aqedup=aqedStat; - aqcdup=aqcdStat; - if(usePowhegBoxConventions) aqedup=-1; - lheStream << setw(7) << int(nup); - lheStream << setw(7) << int(idprup); - lheStream << scientific << setprecision(5) << setw(13) << xwgtup; - lheStream << scientific << setprecision(5) << setw(13) << scalup; - lheStream << scientific << setprecision(5) << setw(13) << aqedup; - lheStream << scientific << setprecision(5) << setw(13) << aqcdup; - lheStream << "\n"; - } - for(int jup = 0; jup < nup; jup++) { - if(debugging>=5) { - lheStream << idup[jup] << "\t" << istup[jup] << "\t" - << mothup1[jup] << "\t" << mothup2[jup] << "\t" - << icolup1[jup] << "\t" << icolup2[jup] << "\t"; - lheStream << pup[0][jup] << "\t" << pup[1][jup] << "\t" - << pup[2][jup] << "\t" << pup[3][jup] << "\t" - << pup[4][jup] << "\t"; - lheStream << "0" << "\t" << "1\n"; - } else { - if(icolup1[jup]!=0) icolup1[jup]+=500; - if(icolup2[jup]!=0) icolup2[jup]+=500; - lheStream << setw(8) << int(idup[jup]); - lheStream << setw(6) << int(istup[jup]) - << setw(6) << int(mothup1[jup]) - << setw(6) << int(mothup2[jup]) - << setw(6) << int(icolup1[jup]) - << setw(6) << int(icolup2[jup]); - lheStream << scientific << setprecision(9) << setw(17) << pup[0][jup]; - lheStream << scientific << setprecision(9) << setw(17) << pup[1][jup]; - lheStream << scientific << setprecision(9) << setw(17) << pup[2][jup]; - lheStream << scientific << setprecision(9) << setw(17) << pup[3][jup]; - lheStream << scientific << setprecision(9) << setw(17) << pup[4][jup]; - lheStream << scientific << setprecision(5) << setw(13) << vtimup[jup]; - lheStream << scientific << setprecision(3) << setw(11) << vtimup[jup]; - lheStream << "\n"; - } - } - lheStream << "\n"; - eventcount++; - cout << "Processed " << eventcount << " events ..." << "\r" << flush; - beginNewEvent=true; - } - } - - cout << "\n\n"; - - if( maxevents!=0 ) { - cout << "All done (" << maxevents << " out of " << unwev << " events).\n"; - } else { - cout << "All done (" << eventcount << " events).\n"; - } - - cout << "\n\n" - << "Wrote a LH event file " << lheFilename - << " and a HW++ MLM merging input file " - << prefix+string(".in") << ".\n\n"; - return 0; - -} - - -inline double sqr(double x) { return x*x; } - -int ndnsToLHAPDF(int ndns) { - // The information in this function is based on - // subroutine PRNTSF from alplib/alppdf.f, LHAPDF's - // PDFsets.index and, finally, the .stat output that - // results when the relevant ndns value is entered - // in the input file. - string Set("no PDF set found"); - double Lambda_4(0),Lambda_5_2loop(0); - string Scheme("no PDF scheme"); - int LHAPDFindex(-999); - string tmpString(""); - - if(ndns==1) { - Set = "CTEQ4M" ; Lambda_4 = 0.298 ; Lambda_5_2loop = 0.202 ; Scheme = "MS" ; - LHAPDFindex = 19150; - } else if(ndns==2) { - Set = "CTEQ4L" ; Lambda_4 = 0.298 ; Lambda_5_2loop = 0.202 ; Scheme = "MS" ; - LHAPDFindex = 19170; - } else if(ndns==3) { - Set = "CTEQ4HJ" ; Lambda_4 = 0.298 ; Lambda_5_2loop = 0.202 ; Scheme = "MS" ; - LHAPDFindex = -99999; - } else if(ndns==4) { - Set = "CTEQ5M" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 19050; - } else if(ndns==5) { - Set = "CTEQ5L" ; Lambda_4 = 0.192 ; Lambda_5_2loop = 0.144 ; Scheme = "MS" ; - LHAPDFindex = 19070; - } else if(ndns==6) { - Set = "CTEQ5HJ" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = -99999; - } else if(ndns==7) { - Set = "CTEQ6M" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex= 10050; - } else if(ndns==8) { - Set = "CTEQ6L" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 10041; - } else if(ndns==9) { - Set = "CTEQ6L1" ; Lambda_4 = 0.215 ; Lambda_5_2loop = 0.167 ; Scheme = "MS" ; - LHAPDFindex = 10042; - } else if(ndns>=10&&ndns<=50) { - Set = "CTEQ6xx" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 10150+(ndns-10); - } else if(ndns==101) { - Set = "MRST99" ; Lambda_4 = 0.321 ; Lambda_5_2loop = 0.220 ; Scheme = "MS" ; - LHAPDFindex = -99999; - } else if(ndns==102) { - Set = "MRST01" ; Lambda_4 = 0.342 ; Lambda_5_2loop = 0.239 ; Scheme = "MS" ; - LHAPDFindex = -99999; - } else if(ndns==103) { - Set = "MRST01" ; Lambda_4 = 0.310 ; Lambda_5_2loop = 0.214 ; Scheme = "MS" ; - LHAPDFindex = -99999; - } else if(ndns==104) { - Set = "MRST01" ; Lambda_4 = 0.378 ; Lambda_5_2loop = 0.267 ; Scheme = "MS" ; - LHAPDFindex = -99999; - } else if(ndns==105) { - Set = "MRST01J" ; Lambda_4 = 0.378 ; Lambda_5_2loop = 0.267 ; Scheme = "MS" ; - LHAPDFindex = -99999; - } else if(ndns==106) { - Set = "MRST02LO" ; Lambda_4 = 0.215 ; Lambda_5_2loop = 0.167 ; Scheme = "MS" ; - LHAPDFindex = -99999; - } else if(ndns==201) { - Set = "MSTW2008lo" ; Lambda_4 = 0.322 ; Lambda_5_2loop = 0.255 ; Scheme = "MS" ; - LHAPDFindex = 21000; - } else if(ndns==202) { - Set = "MSTW2008nlo" ; Lambda_4 = 0.365 ; Lambda_5_2loop = 0.255 ; Scheme = "MS" ; - LHAPDFindex = 21100; - } else if(ndns>=203&&ndns<=242) { - Set = "MSTW2008lo68cl"; Lambda_4 = 0.322 ; Lambda_5_2loop = 0.255 ; Scheme = "MS" ; - LHAPDFindex = 21000+(ndns-202); - } else if(ndns==243) { - Set = "MRST LO*" ; Lambda_4 = 0.365 ; Lambda_5_2loop = 0.255 ; Scheme = "MS" ; - LHAPDFindex = 20650; - } else if(ndns==244) { - Set = "MRST LO**" ; Lambda_4 = 0.280 ; Lambda_5_2loop = 0.190 ; Scheme = "MS" ; - LHAPDFindex = 20651; - } else if(ndns==301 ) { - Set = "CTQ6.6" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 10550; - } else if(ndns>=302&&ndns<=345) { - Set = "CTQ66" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 10550+(ndns-301); - } else if(ndns==346) { - Set = "CT09MC1" ; Lambda_4 = 0.215 ; Lambda_5_2loop = 0.167 ; Scheme = "MS" ; - LHAPDFindex = 10771; - } else if(ndns==347) { - Set = "CT09MC2" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 10772; - } - - /*cout << "-------------------------------\n"; - cout << "ndnsToLHAPDF found: \n"; - cout << "PDF set = " << Set << "\n"; - cout << "ndns index = " << ndns << "\n"; - cout << "LHAPDF index = " << LHAPDFindex << "\n"; - cout << "-------------------------------\n\n";*/ - return LHAPDFindex; - -} - - -string ndnsToLHAPDF_str(int ndns) { - // The information in this function is based on - // subroutine PRNTSF from alplib/alppdf.f, LHAPDF's - // PDFsets.index and, finally, the .stat output that - // results when the relevant ndns value is entered - // in the input file. - string Set("no PDF set found"); - double Lambda_4(0),Lambda_5_2loop(0); - string Scheme("no PDF scheme"); - int LHAPDFindex(-999); - string tmpString(""); - string lhastring(""); - - if(ndns==1) { - Set = "CTEQ4M" ; Lambda_4 = 0.298 ; Lambda_5_2loop = 0.202 ; Scheme = "MS" ; - LHAPDFindex = 19150; - lhastring = "cteq6"; - } else if(ndns==2) { - Set = "CTEQ4L" ; Lambda_4 = 0.298 ; Lambda_5_2loop = 0.202 ; Scheme = "MS" ; - LHAPDFindex = 19170; - lhastring = "cteq6l1"; - } else if(ndns==3) { - Set = "CTEQ4HJ" ; Lambda_4 = 0.298 ; Lambda_5_2loop = 0.202 ; Scheme = "MS" ; - LHAPDFindex = -99999; - lhastring = "cteq6l1"; - } else if(ndns==4) { - Set = "CTEQ5M" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 19050; - lhastring = "cteq6l1"; - } else if(ndns==5) { - Set = "CTEQ5L" ; Lambda_4 = 0.192 ; Lambda_5_2loop = 0.144 ; Scheme = "MS" ; - LHAPDFindex = 19070; - lhastring = "cteq6l1"; - - } else if(ndns==6) { - Set = "CTEQ5HJ" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = -99999; - lhastring = "cteq6l1"; - } else if(ndns==7) { - Set = "CTEQ6M" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex= 10050; - lhastring = "cteq6"; - } else if(ndns==8) { - Set = "CTEQ6L" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 10041; - lhastring = "cteq6l1"; - } else if(ndns==9) { - Set = "CTEQ6L1" ; Lambda_4 = 0.215 ; Lambda_5_2loop = 0.167 ; Scheme = "MS" ; - LHAPDFindex = 10042; - lhastring = "cteq6l1"; - } else if(ndns>=10&&ndns<=50) { - Set = "CTEQ6xx" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 10150+(ndns-10); - lhastring = "cteq6l1"; - } else if(ndns==101) { - Set = "MRST99" ; Lambda_4 = 0.321 ; Lambda_5_2loop = 0.220 ; Scheme = "MS" ; - LHAPDFindex = -99999; - lhastring = "mrst99"; - } else if(ndns==102) { - Set = "MRST01" ; Lambda_4 = 0.342 ; Lambda_5_2loop = 0.239 ; Scheme = "MS" ; - LHAPDFindex = -99999; - lhastring = "mrst01"; - } else if(ndns==103) { - Set = "MRST01" ; Lambda_4 = 0.310 ; Lambda_5_2loop = 0.214 ; Scheme = "MS" ; - LHAPDFindex = -99999; - } else if(ndns==104) { - Set = "MRST01" ; Lambda_4 = 0.378 ; Lambda_5_2loop = 0.267 ; Scheme = "MS" ; - LHAPDFindex = -99999; - lhastring = "mrst01"; - } else if(ndns==105) { - Set = "MRST01J" ; Lambda_4 = 0.378 ; Lambda_5_2loop = 0.267 ; Scheme = "MS" ; - LHAPDFindex = -99999; - lhastring = "mrst01j"; - } else if(ndns==106) { - Set = "MRST02LO" ; Lambda_4 = 0.215 ; Lambda_5_2loop = 0.167 ; Scheme = "MS" ; - LHAPDFindex = -99999; - lhastring = "mrst02lo"; - } else if(ndns==201) { - Set = "MSTW2008lo" ; Lambda_4 = 0.322 ; Lambda_5_2loop = 0.255 ; Scheme = "MS" ; - LHAPDFindex = 21000; - } else if(ndns==202) { - Set = "MSTW2008nlo" ; Lambda_4 = 0.365 ; Lambda_5_2loop = 0.255 ; Scheme = "MS" ; - LHAPDFindex = 21100; - lhastring = "mstw2008nlo"; - } else if(ndns>=203&&ndns<=242) { - Set = "MSTW2008lo68cl"; Lambda_4 = 0.322 ; Lambda_5_2loop = 0.255 ; Scheme = "MS" ; - LHAPDFindex = 21000+(ndns-202); - lhastring = "mstw2008lo68cl"; - } else if(ndns==243) { - Set = "MRST LO*" ; Lambda_4 = 0.365 ; Lambda_5_2loop = 0.255 ; Scheme = "MS" ; - LHAPDFindex = 20650; - lhastring = "MRST2007lomod"; - } else if(ndns==244) { - Set = "MRST LO**" ; Lambda_4 = 0.280 ; Lambda_5_2loop = 0.190 ; Scheme = "MS" ; - LHAPDFindex = 20651; - lhastring = "MRSTMCal"; - - } else if(ndns==301 ) { - Set = "CTQ6.6" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 10550; - lhastring = "cteq66"; - } else if(ndns>=302&&ndns<=345) { - Set = "CTQ66" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 10550+(ndns-301); - lhastring = "cteq66"; - } else if(ndns==346) { - Set = "CT09MC1" ; Lambda_4 = 0.215 ; Lambda_5_2loop = 0.167 ; Scheme = "MS" ; - LHAPDFindex = 10771; - lhastring = "CT09MC1"; - } else if(ndns==347) { - Set = "CT09MC2" ; Lambda_4 = 0.326 ; Lambda_5_2loop = 0.226 ; Scheme = "MS" ; - LHAPDFindex = 10772; - lhastring = "CT09MC2"; - } - - cout << "-------------------------------\n"; - cout << "ndnsToLHAPDF found: \n"; - cout << "PDF set = " << Set << "\n"; - cout << "ndns index = " << ndns << "\n"; - cout << "LHAPDF index = " << LHAPDFindex << "\n"; - cout << "WARNING: YOU MAY HAVE TO ENTER THE PDF NAME MANUALLY IN THE INPUT FILES!" << endl; - cout << "-------------------------------\n\n"; - return lhastring; - -} - - -double parstrToparval(string varName, - vector * parstrPtr, - vector * parvalPtr) { - for(unsigned int index=0; indexsize(); index++) - if(varName==parstrPtr->at(index)) - return parvalPtr->at(index); - - return -999.0; -} - -string trim(string theString) { - int endStr = theString.find_last_not_of(" "); - int beginStr = theString.find_first_not_of(" "); - if(beginStr==0&&endStr==theString.length()-1) return theString; // No lead / trail spaces. - theString = theString.substr(beginStr,endStr-beginStr+1); - return theString; -} - -void writeHWPPinFile(string prefix, int ihrd, int unwev, - int lhapdf, string lhapdfstr, int idbmup0, int idbmup1, int idwtup, - double aqcdup, int nloop, - vector * massesPtr, - vector * parvalPtr, - vector * parstrPtr) { - ofstream hwpp; - hwpp.open(string(prefix+".in").c_str()); - - hwpp << "#############################################################\n"; - hwpp << "# Create an event generator taking the default EventGenerator #\n"; - hwpp << "# as the starting point ... #\n"; - hwpp << "#############################################################\n"; - hwpp << "cd /Herwig/Generators\n"; - hwpp << "# Copy the default EventGenerator with its settings to a new \n"; - hwpp << "# which will be the basis of the one we use for showering: \n"; - hwpp << "cp EventGenerator theGenerator\n"; - hwpp << "\n"; - hwpp << "#############################################################\n"; - hwpp << "# Create a LH event handler (set up & assigned below) ... #\n"; - hwpp << "#############################################################\n"; - hwpp << "cd /Herwig/EventHandlers\n"; - hwpp << "library LesHouches.so\n"; - hwpp << "create ThePEG::LesHouchesEventHandler theLesHouchesHandler\n"; - hwpp << "\n"; - hwpp << "#############################################################\n"; - hwpp << "# Create a LH reader (set up & assigned below) ... #\n"; - hwpp << "#############################################################\n"; - hwpp << "cd /Herwig/EventHandlers\n"; - hwpp << "library LesHouches.so\n"; - hwpp << "create ThePEG::LesHouchesFileReader theLHReader\n"; - hwpp << "\n"; - hwpp << "#############################################################\n"; - hwpp << "# Create an AlpGenHandler (set up & assigned below) ... #\n"; - hwpp << "#############################################################\n"; - hwpp << "cd /Herwig/Shower\n"; - hwpp << "library AlpGenHandler.so\n"; - hwpp << "create Herwig::AlpGenHandler AlpGenHandler\n"; - hwpp << "set /Herwig/Shower/AlpGenHandler:ShowerModel /Herwig/Shower/ShowerModel\n"; - hwpp << "set /Herwig/Shower/AlpGenHandler:SplittingGenerator /Herwig/Shower/SplittingGenerator\n"; - hwpp << "\n"; - hwpp << "#############################################################\n"; - hwpp << "# Create an LHAPDF (set up & assigned below) ... #\n"; - hwpp << "#############################################################\n"; - hwpp << "cd /Herwig/Partons\n"; - hwpp << "create ThePEG::LHAPDF thePDFset ThePEGLHAPDF.so\n"; - hwpp << "\n"; - hwpp << "############################################################\n"; - hwpp << "# Create a cuts object ... #\n"; - hwpp << "############################################################\n"; - hwpp << "cd /Herwig/EventHandlers\n"; - hwpp << "create ThePEG::Cuts /Herwig/Cuts/NoCuts\n"; - hwpp << "\n"; - hwpp << "#############################################################\n"; - hwpp << "# Setup the LH event handler ... #\n"; - hwpp << "#############################################################\n"; - hwpp << "cd /Herwig/EventHandlers\n"; - hwpp << "insert theLesHouchesHandler:LesHouchesReaders 0 theLHReader\n"; - if(idwtup==3) { - hwpp << "set theLesHouchesHandler:WeightOption UnitWeight\n"; - } else if(idwtup==-3) { - hwpp << "set theLesHouchesHandler:WeightOption NegUnitWeight\n"; - } else if(idwtup==4) { - hwpp << "set theLesHouchesHandler:WeightOption VarWeight\n"; - } else { - hwpp << "set theLesHouchesHandler:WeightOption VarNegWeight\n"; - } - // hwpp << "set theLesHouchesHandler:PartonExtractor " - // << "/Herwig/Partons/QCDExtractor\n"; - hwpp << "set theLesHouchesHandler:PartonExtractor /Herwig/Partons/PPExtractor\n"; - - hwpp << "set theLesHouchesHandler:CascadeHandler " - << "/Herwig/Shower/AlpGenHandler\n"; - hwpp << "set theLesHouchesHandler:HadronizationHandler " - << "/Herwig/Hadronization/ClusterHadHandler\n"; - hwpp << "set theLesHouchesHandler:DecayHandler " - << "/Herwig/Decays/DecayHandler\n"; - hwpp << "\n"; - /* hwpp << "#############################################################\n"; - hwpp << "# Set up the Evolver to veto hard emissions > scalup ... #\n"; - hwpp << "#############################################################\n"; - hwpp << "cd /Herwig/Shower\n"; - hwpp << "# MaxTry 100 sets the maximum number of times to try \n"; - hwpp << "# showering a given shower tree to 100. \n"; - hwpp << "# HardVetoMode Yes to veto emissions with pT greater than pT_max.\n"; - hwpp << "# HardVetoScaleSource Read means pT_max comes from hepeup.SCALUP.\n"; - hwpp << "# This is what you need to set _along_with_ HardVetoMode Yes in \n"; - hwpp << "# the case of Powheg external events _AND_ mc@nlo (we know this \n"; - hwpp << "# from looking at the *MCinput file that mc@nlo generates). \n"; - hwpp << "# MeCorrMode No turns off ME corrs. IntrinsicPtGaussian 2.2*GeV \n"; - hwpp << "# is the RMS of intrinsic pT of Gaussian distribution: \n"; - hwpp << "# 2*(1-Beta)*exp(-sqr(intrinsicpT/RMS))/sqr(RMS) \n"; - // hwpp << "set Evolver:MaxTry 100\n"; - // hwpp << "set Evolver:HardVetoMode Yes\n"; - // hwpp << "set Evolver:HardVetoScaleSource Read\n"; - // hwpp << "set Evolver:HardVetoReadOption PrimaryCollision\n"; - // hwpp << "set Evolver:MECorrMode No\n"; - hwpp << "# Intrinsic pT etc should be set as part of a tune i.e. it \n"; - hwpp << "# should either be left alone (default) or set by reading in \n"; - hwpp << "# one of the tunes before theGenerator is created by copying \n"; - hwpp << "# EventGenerator (second line of this file). The following \n"; - hwpp << "# settings were extracted from LHC-UE-EE-3-CTEQ6L1.in - In \n"; - hwpp << "# light of some bad experience with MPI we prefer to set \n"; - hwpp << "# these here manually rather than read try to read that .in. \n"; - if(idbmup0 == 2212 && idbmup1 == -2212) { - hwpp << "#set /Herwig/UnderlyingEvent/KtCut:MinKT 2.26 \n"; - hwpp << "#set /Herwig/UnderlyingEvent/UECuts:MHatMin 4.52 \n"; - hwpp << "#set /Herwig/Shower/ShowerHandler:IntrinsicPtGaussian 1.9*GeV \n"; - } else { - hwpp << "#set /Herwig/UnderlyingEvent/KtCut:MinKT 2.752 \n"; - hwpp << "#set /Herwig/UnderlyingEvent/UECuts:MHatMin 5.504 \n"; - hwpp << "#set /Herwig/Shower/Evolver:IntrinsicPtGaussian 2.34*GeV \n"; - } - hwpp << "# Colour reconnection (re)settings \n"; - hwpp << "#set /Herwig/Hadronization/ColourReconnector:ColourReconnection Yes \n"; - hwpp << "#set /Herwig/Hadronization/ColourReconnector:ReconnectionProbability 0.61\n"; - hwpp << "# Colour Disrupt settings \n"; - hwpp << "#set /Herwig/Partons/RemnantDecayer:colourDisrupt 0.75 \n"; - hwpp << "# Inverse hadron radius \n"; - hwpp << "#set /Herwig/UnderlyingEvent/MPIHandler:InvRadius 1.35 \n"; - hwpp << "#set /Herwig/UnderlyingEvent/MPIHandler:softInt Yes \n"; - hwpp << "#set /Herwig/UnderlyingEvent/MPIHandler:twoComp Yes \n"; - hwpp << "#set /Herwig/UnderlyingEvent/MPIHandler:DLmode 2 \n"; - hwpp << "\n"; - hwpp << "#############################################################\n"; - hwpp << "# Set up kinematics reconstructor (relevant only to mc@nlo) #\n"; - hwpp << "#############################################################\n"; - hwpp << "# Options for QTildeReconstructor - not needed for Powheg but\n"; - hwpp << "# critical for mc@nlo. If using Powheg you may either leave \n"; - hwpp << "# the next two settings as they are or comment them out as \n"; - hwpp << "# you wish, for mc@nlo though you must leave them as they \n"; - hwpp << "# are! ReconstructionOption General was the old default, it \n"; - hwpp << "# ignores the colour structure for all processes - mc@nlo \n"; - hwpp << "# will give you garbage unless you set this! \n"; - hwpp << "# ReconstructionOption Colour is the new default - use the \n"; - hwpp << "# colour structure of the process to determine the \n"; - hwpp << "# reconstruction procedure. InitialInitialBoostOption \n"; - hwpp << "# determines how the boost from the system before ISR to that\n"; - hwpp << "# after ISR is applied. mc@nlo requires the old kinematics \n"; - hwpp << "# reconstruction method: \n"; - hwpp << "# InitialInitialBoostOption LongTransBoost - first apply a \n"; - hwpp << "# longitudinal and then a transverse boost. Whereas the \n"; - hwpp << "# default method now is InitialInitialBoostOption OneBoost - \n"; - hwpp << "# apply one boost from old CMS to new CMS. Both options \n"; - hwpp << "# should work with Powheg but probably it's best to use the \n"; - hwpp << "# defaults in that case by simply commenting this setting. \n"; - hwpp << "#set KinematicsReconstructor:ReconstructionOption General \n"; - hwpp << "#set KinematicsReconstructor:InitialInitialBoostOption LongTransBoost\n"; - hwpp << "\n";*/ - hwpp << "#############################################################\n"; - hwpp << "# Set up the AlpGenHandler ... #\n"; - hwpp << "#############################################################\n"; - hwpp << "cd /Herwig/Shower\n"; - hwpp << "set AlpGenHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler\n"; - hwpp << "set AlpGenHandler:RemDecayer /Herwig/Partons/RemnantDecayer\n"; - hwpp << "set /Herwig/Shower/ShowerHandler:MaxPtIsMuF Yes\n"; - hwpp << "set /Herwig/Shower/ShowerHandler:RestrictPhasespace Yes\n"; - hwpp << "set /Herwig/Shower/ShowerHandler:MaxTry 100\n"; - hwpp << "set /Herwig/Shower/PartnerFinder:PartnerMethod Random\n"; - hwpp << "set /Herwig/Shower/PartnerFinder:ScaleChoice Partner\n"; - hwpp << "set /Herwig/Shower/ShowerHandler:Interactions QCD\n"; - hwpp << "set /Herwig/Shower/ShowerHandler:SpinCorrelations No\n"; - hwpp << "set /Herwig/Shower/ShowerHandler:SoftCorrelations No\n"; - hwpp << "set /Herwig/Shower/ShowerHandler:ReconstructionOption CutOff\n"; - hwpp << "set /Herwig/Shower/ShowerHandler:Interactions QCD\n"; - - // hwpp << "set AlpGenHandler:Evolver Evolver\n"; - // hwpp << "set AlphaQCD:AlphaMZ " << aqcdup << "\n"; - // hwpp << "set AlphaQCD:NumberOfLoops " << nloop << "\n"; - hwpp << "set AlpGenHandler:ShowerAlpha AlphaQCD\n"; - hwpp << "# Calorimeter granularity settings used by GetJet algorithm\n"; - hwpp << "set AlpGenHandler:NoCellsInRapidity 100\n"; - hwpp << "set AlpGenHandler:NoCellsInPhi 60\n"; - // AlpGen hard process code. - hwpp << "# AlpGen hard process code.\n"; - hwpp << "set AlpGenHandler:ihrd " << ihrd << "\n"; - // Number of (light) jets. - int njets(int(parstrToparval("njets",parstrPtr,parvalPtr))); - hwpp << "# No. of light jets in AlpGen process (the \"extra\" ones).\n"; - hwpp << "set AlpGenHandler:njets " << njets << "\n"; - // Mimimum jet pT use for generation. - double ptjmin(parstrToparval("ptjmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:ptjmin " << ptjmin << "*GeV\n"; - // Mimimum parton-parton R-sep used for generation. - double drjmin(parstrToparval("drjmin",parstrPtr,parvalPtr)); - hwpp << "# Mimimum parton-parton R-sep used for generation.\n"; - hwpp << "set AlpGenHandler:drjmin " << drjmin << "\n"; - // Max |eta| for partons in generation. - double etajmax(parstrToparval("etajmax",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:etajmax " << etajmax << "\n"; - // Also probably want these variables fed to AlpGenHandler too --- - // they get set in the alpsho.f AHspar routine (note the list below - // does not include some variables from AHspar because they are already - // included in the above eg the PDFs are already handled so I removed - // ndns and also ptjmin drjmin are written out for the AlpGenHandler above). - int ickkw(int(parstrToparval("ickkw",parstrPtr,parvalPtr))); - //// hwpp << "set AlpGenHandler:ickkw " << ickkw << "\n"; - int ihvy(int(parstrToparval("ihvy",parstrPtr,parvalPtr))); - hwpp << "# heavy flavour in WQQ,ZQQ,2Q etc (4=c, 5=b, 6=t):\n"; - hwpp << "set AlpGenHandler:ihvy " << ihvy << "\n"; - int ihvy2(int(parstrToparval("ihvy2",parstrPtr,parvalPtr))); - //// hwpp << "set AlpGenHandler:ihvy2 " << ihvy2 << "\n"; - int itopprc(nInt(parstrToparval("itopprc",parstrPtr,parvalPtr))); - //// hwpp << "set AlpGenHandler:itopprc " << itopprc << "\n"; - int nw(int(parstrToparval("nw",parstrPtr,parvalPtr))); - if(ihrd==13) { - nw=1; // N.B. nw is reassigned in this way - if(itopprc>=3) nw=2; // by UPEVNT (after UPINIT). - } - //// hwpp << "set AlpGenHandler:nw " << nw << "\n"; - int nz(int(parstrToparval("nz",parstrPtr,parvalPtr))); - //// hwpp << "set AlpGenHandler:nz " << nz << "\n"; - int nh(int(parstrToparval("nh",parstrPtr,parvalPtr))); - hwpp << "# Number of Higgses in the AlpGen process:\n"; - hwpp << "set AlpGenHandler:nh " << nh << "\n"; - int nph(int(parstrToparval("nph",parstrPtr,parvalPtr))); - hwpp << "# Number of photons in the AlpGen process:\n"; - hwpp << "set AlpGenHandler:nph " << nph << "\n"; - double ptbmin(parstrToparval("ptbmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:ptbmin " << ptbmin << "\n"; - double ptcmin(parstrToparval("ptcmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:ptcmin " << ptcmin << "\n"; - double ptlmin(parstrToparval("ptlmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:ptlmin " << ptlmin << "\n"; - double metmin(parstrToparval("metmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:metmin " << metmin << "\n"; - double ptphmin(parstrToparval("ptphmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:ptphmin " << ptphmin << "\n"; - double etabmax(parstrToparval("etabmax",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:etabmax " << etabmax << "\n"; - double etacmax(parstrToparval("etacmax",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:etacmax " << etacmax << "\n"; - double etalmax(parstrToparval("etalmax",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:etalmax " << etalmax << "\n"; - double etaphmax(parstrToparval("etaphmax",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:etaphmax " << etaphmax << "\n"; - double drbmin(parstrToparval("drbmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:drbmin " << drbmin << "\n"; - double drcmin(parstrToparval("drcmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:drcmin " << drcmin << "\n"; - double drlmin(parstrToparval("drlmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:drlmin " << drlmin << "\n"; - double drphjmin(parstrToparval("drphjmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:drphjmin " << drphjmin << "\n"; - double drphlmin(parstrToparval("drphlmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:drphlmin " << drphlmin << "\n"; - double drphmin(parstrToparval("drphmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:drphmin " << drphmin << "\n"; - double mllmin(parstrToparval("mllmin",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:mllmin " << mllmin << "\n"; - double mllmax(parstrToparval("mllmax",parstrPtr,parvalPtr)); - //// hwpp << "set AlpGenHandler:mllmax " << mllmax << "\n"; - hwpp << "\n"; - hwpp << "#############################################################\n"; - hwpp << "# Set up the LH reader ... #\n"; - hwpp << "#############################################################\n"; - hwpp << "cd /Herwig/EventHandlers\n"; - hwpp << "set theLHReader:WeightWarnings false\n"; - hwpp << "# Input event file name:\n"; - hwpp << "set theLHReader:FileName " << prefix << ".lhe\n"; - hwpp << "set theLHReader:MomentumTreatment RescaleEnergy\n"; - hwpp << "# set theLHReader:IgnoreIDPRUP Yes\n"; - hwpp << "set theLHReader:Cuts /Herwig/Cuts/NoCuts\n"; - hwpp << "#set theLHReader:OverSampling ForbidOverSampling\n"; - hwpp << "\n"; - hwpp << "#############################################################\n"; - hwpp << "# Set up the LHAPDF ... #\n"; - hwpp << "#############################################################\n"; - hwpp << "cd /Herwig/Partons\n"; - hwpp << "# Don't try and find PDF index out from the LH file ...\n"; - hwpp << "set /Herwig/EventHandlers/theLHReader:InitPDFs false\n"; - hwpp << "# Instead set them explicitly here:\n"; - // hwpp << "set thePDFset:PDFNumber " << lhapdf << "\n"; - hwpp << "set thePDFset:PDFName " << lhapdfstr << "\n"; - - hwpp << "set thePDFset:RemnantHandler HadronRemnants\n"; - hwpp << "set /Herwig/EventHandlers/theLHReader:PDFA thePDFset\n"; - hwpp << "set /Herwig/EventHandlers/theLHReader:PDFB thePDFset\n"; - hwpp << "set /Herwig/Particles/p+:PDF thePDFset\n"; - hwpp << "set /Herwig/Particles/pbar-:PDF thePDFset\n"; - hwpp << "# The PDF for beam particles A/B - overrides particle's own " - << "PDF above\n"; - hwpp << "set /Herwig/Shower/AlpGenHandler:PDFA thePDFset\n"; - hwpp << "set /Herwig/Shower/AlpGenHandler:PDFB thePDFset\n"; - hwpp << "set /Herwig/Shower/ShowerHandler:PDFA thePDFset\n"; - hwpp << "set /Herwig/Shower/ShowerHandler:PDFB thePDFset\n"; - hwpp << "\n"; - hwpp << "####################################################\n"; - hwpp << "# Set up the generator ... #\n"; - hwpp << "####################################################\n"; - hwpp << "cd /Herwig/Generators\n"; - hwpp << "set theGenerator:EventHandler " - << "/Herwig/EventHandlers/theLesHouchesHandler\n"; - hwpp << "set theGenerator:NumberOfEvents " << unwev << "\n"; - hwpp << "set theGenerator:RandomNumberGenerator:Seed 31122001\n"; - hwpp << "set theGenerator:PrintEvent 10\n"; - hwpp << "set theGenerator:MaxErrors 10000\n"; - hwpp << "\n"; - hwpp << "###################################################################\n"; - hwpp << "# ReDefine particle data like it is in the AlpGen parameter file. #\n"; - hwpp << "###################################################################\n"; - hwpp << "\ncd /Herwig/Particles/ \n"; - // 'if' statements needed here to protect against mcat(0)>1.0) { - hwpp << "set c:NominalMass " << massesPtr->at(0) << "*GeV\n"; - hwpp << "set cbar:NominalMass " << massesPtr->at(0) << "*GeV\n"; - } - // Ditto. - if(massesPtr->at(0)>4.0) { - hwpp << "set b:NominalMass " << massesPtr->at(1) << "*GeV\n"; - hwpp << "set bbar:NominalMass " << massesPtr->at(1) << "*GeV\n"; - } - hwpp << "set t:NominalMass " << massesPtr->at(2) << "*GeV\n"; - hwpp << "set tbar:NominalMass " << massesPtr->at(2) << "*GeV\n"; - hwpp << "set W+:NominalMass " << massesPtr->at(3) << "*GeV\n"; - hwpp << "set W-:NominalMass " << massesPtr->at(3) << "*GeV\n"; - hwpp << "set Z0:NominalMass " << massesPtr->at(4) << "*GeV\n"; - hwpp << "set h0:NominalMass " << massesPtr->at(5) << "*GeV\n"; - hwpp << "\n\n\n\n\n"; - hwpp << "######################################################### \n"; - hwpp << "######################################################### \n"; - hwpp << "## # \n"; - hwpp << "## --- USER SERVICEABLE PART BELOW HERE ONLY ! --- # \n"; - hwpp << "## # \n"; - hwpp << "######################################################### \n"; - hwpp << "######################################################### \n"; - hwpp << "\n\n\n\n\n"; - hwpp << "######################################################### \n"; - hwpp << "# Option to off shower / hadronization / decays / MPI. # \n"; - hwpp << "######################################################### \n"; - hwpp << "cd /Herwig/EventHandlers \n"; - hwpp << "# set theLesHouchesHandler:CascadeHandler NULL \n"; - hwpp << "#set theLesHouchesHandler:HadronizationHandler NULL \n"; - hwpp << "#set theLesHouchesHandler:DecayHandler NULL \n"; - hwpp << "# The handler for multiple parton interactions \n"; - hwpp << "#set /Herwig/Shower/AlpGenHandler:MPIHandler NULL \n"; - hwpp << "\n\n"; - hwpp << "######################################################### \n"; - hwpp << "# Recommended key MLM merging parameters below - change # \n"; - hwpp << "# for systematic error studies and / or at your peril. # \n"; - hwpp << "######################################################### \n"; - hwpp << "cd /Herwig/Shower\n"; - hwpp << "# Is this the highest multiplicity ME in merging? \n"; - hwpp << "# 0 = no, 1 = yes . \n"; - hwpp << "set AlpGenHandler:highestMultiplicity 0 \n"; - hwpp << "# Jet ET cut to apply in jet clustering in merging.\n"; - double etclus(max(ptjmin+5,1.2*ptjmin)); - hwpp << "set AlpGenHandler:ETClus " << etclus << "*GeV\n"; - hwpp << "# Cone size used in clustering in merging.\n"; - double rclus(drjmin); - hwpp << "set AlpGenHandler:RClus " << rclus << "\n"; - hwpp << "# Max |eta| for jets in clustering in merging.\n"; - double etaclmax(etajmax); - hwpp << "set AlpGenHandler:EtaClusMax " << etaclmax << "\n"; - hwpp << "# Default 1.5 factor used to decide if a jet matches a parton\n"; - hwpp << "# in merging: if DR(parton,jet) * idup , vector * istup, - vector * mothup1, vector * mothup2, - vector * icolup1, vector * icolup2, - vector > * pup, - vector masses , int itopprc) { - - int iwch(0); - - // W/Z/gamma b bbar + jets ( wcjet*, ihrd=10 / wphjet*, ihrd=14 / wphqq*, - // ihrd=15 ), or W/Z + jets ( wjet*, ihrd=3 / zjet*, ihrd=4 ): - // In this case we add to the list of particles the single intermediate - // boson (at the end of the list) appropriately, and assign the relevant - // parent-daughter and colour flow indices. - if (ihrd<=4||ihrd==10||ihrd==14||ihrd==15) { - iwch=0; // <--- used to determine type: W/Z/gamma - for(int iup=int(*nup)-2; iup 11 + (-12) => -(1)+0 = -1 => W- - // positron+nu -> -11+ 12 => -(-1)+0 = -1 => W+ - // u dbar -> 2 -1 => 0 -(-1) = 1 => W+ - // c dbar -> 4 -1 => W+ - // etc. - } - // Now we start adding the intermediate boson entry: - int iup=nInt(*nup); - if(iwch>0) (*idup).push_back( 24); - else if(iwch<0) (*idup).push_back(-24); - else (*idup).push_back( 23); - (*istup).push_back(2); - (*mothup1).push_back(1); - (*mothup2).push_back(2); - (*pup).push_back(vector(5)); - double tmp = (*pup)[3][iup-2]+(*pup)[3][iup-1]; // Vector boson energy. - (*pup)[3][iup] = tmp; - tmp = sqr(tmp); - for(unsigned int ixx=0; ixx<=2; ixx++) { - (*pup)[ixx][iup] = (*pup)[ixx][iup-2]+(*pup)[ixx][iup-1]; - tmp = tmp-sqr((*pup)[ixx][iup]); // Vector boson mass^2 when loop ends. - } - (*pup)[4][iup] = sqrt(tmp); // Set vector boson mass. - (*icolup1).push_back(0); // Set 1st colour line for vector boson. - (*icolup2).push_back(0); // Set 2nd colour line for vector boson. - (*nup) = (*nup)+1; // Increment number of particles to be read in the event. - } - // nW + mZ + kH + jets ( vbjet* / ihrd=5 ): - else if(ihrd==5) { - unsigned int ivstart(0),ivend(0); - // Identify the range of the Z and W bosons in the event (AlpGen conventions). - // Note the Z's and W's are the only things that decay - Higgs's and photons - // do not. - vector bosonIndex; - for(unsigned int ixx=0; ixx<(*nup); ixx++) - if(abs((*idup)[ixx])==24||(*idup)[ixx]==23) { - (*istup)[ixx]=2; // Set the W/Z boson status to intermediate. - bosonIndex.push_back(ixx+1); // W/Z boson index in LH event record (1->nup). - } - unsigned int bosonCounter(nInt(*nup)-2*bosonIndex.size()); - for(unsigned int ixx=0; ixx(5)); - } - for(unsigned int iw=1; iw<=2; iw++) { - int iwdec(nInt(*nup)-5+2*iw); // First of decay products for this W (iw). - int iwup(nInt(*nup)+iw); // Where the reco. W will go in evt.record (LH index). - int ibup(iwup+2); // Where the reco. b will go in evt.record (under the Ws). - int iwch(0); - for(unsigned int iup=iwdec; iup<=iwdec+1; iup++) { - (*mothup1)[iup-1]=iwup; - (*mothup2)[iup-1]=0; - iwch=iwch-int((*idup)[iup-1])%2; // iwch = charge of W boson. - // electron+nubar -> 11 + (-12) = -1 => W- - // d + ubar -> 1 + (-2) = -1 => W- - // positron+nu -> -11+ 12 = 1 => W+ - // u + dbar -> 2 + (-1) = 1 => W+ - } - // Make space for the b and W: - // Fill in b and W LH record entries: - if(iwch>0) { - (*idup)[iwup-1]=24; - (*idup)[ibup-1]=5; - (*mothup1)[iwup-1]=it; - (*mothup2)[iwup-1]=0; - (*mothup1)[ibup-1]=it; - (*mothup2)[ibup-1]=0; - } else if (iwch<0) { - (*idup)[iwup-1]=-24; - (*idup)[ibup-1]=-5; - (*mothup1)[iwup-1]=itb; - (*mothup2)[iwup-1]=0; - (*mothup1)[ibup-1]=itb; - (*mothup2)[ibup-1]=0; - } - (*istup)[iwup-1]=2; // The W is an intermediate and the - (*istup)[ibup-1]=1; // b is a final-state particle. - // Now reconstruct W boson momentum - double tmp=(*pup)[3][iwdec-1]+(*pup)[3][iwdec]; - (*pup)[3][iwup-1]=tmp; // W energy. - tmp=sqr(tmp); - for(unsigned int ixx=0; ixx<=2; ixx++) { - (*pup)[ixx][iwup-1] = (*pup)[ixx][iwdec-1] // W's 3-mom. - + (*pup)[ixx][iwdec]; - tmp=tmp-sqr((*pup)[ixx][iwup-1]); // Equals m^2 at the end of the loop. - } - (*pup)[4][iwup-1]=sqrt(tmp); // W mass. - // Reconstruct b momentum - int itmp(nInt((*mothup1)[iwup-1])); - tmp=(*pup)[3][itmp-1]-(*pup)[3][iwup-1]; - (*pup)[3][ibup-1]=tmp; // b energy. - tmp=sqr(tmp); - for(unsigned int ixx=0; ixx<=2; ixx++) { - (*pup)[ixx][ibup-1] = (*pup)[ixx][nInt((*mothup1)[iwup-1])-1] // b's 3mom. - - (*pup)[ixx][iwup-1]; - tmp=tmp-sqr((*pup)[ixx][ibup-1]); // Equals m^2 at the end of the loop. - } - (*pup)[4][ibup-1]=sqrt(tmp); // b mass. - (*icolup1)[iwup-1]=0; // W has no colour - (*icolup2)[iwup-1]=0; // lines. - (*icolup1)[ibup-1]=(*icolup1)[nInt((*mothup1)[iwup-1])-1]; // b shares top - (*icolup2)[ibup-1]=(*icolup2)[nInt((*mothup1)[iwup-1])-1]; // colour line. - } - (*nup)+=4; - } - // H t tbar + jets ( QQh*, ihrd=8 ): - else if (ihrd==8&&abs((*idup)[3])==6) { - // Redefine the tops as intermediates in the event record. - (*istup)[3]=2; - (*istup)[4]=2; - unsigned int it(4),itb(5); // Index of t & tbar in evt.record (LH index). - if((*idup)[3]!=6) swap(it,itb); - // Reconstruct intermediate W's from the decay products. - for(unsigned int ixx=0; ixx<4;ixx++) { - (*idup).push_back(0); - (*istup).push_back(0); - (*mothup1).push_back(0); - (*mothup2).push_back(0); - (*icolup1).push_back(0); - (*icolup2).push_back(0); - (*pup).push_back(vector(5)); - } - for(unsigned int iw=1; iw<=2; iw++) { - int iwdec(nInt(*nup)-5+2*iw); // First of decay products for this W (iw). - int iwup(nInt(*nup)+iw); // Where the reco. W will go in evt.record (LH index). - int ibup(iwup+2); // Where the reco. b will go in evt.record (under the Ws). - int iwch(0); - for(unsigned int iup=iwdec; iup<=iwdec+1; iup++) { - (*mothup1)[iup-1]=iwup; - (*mothup2)[iup-1]=0; - iwch=iwch-int((*idup)[iup-1])%2; // iwch = charge of W boson. - // electron+nubar -> 11 + (-12) = -1 => W- - // d + ubar -> 1 + (-2) = -1 => W- - // positron+nu -> -11+ 12 = 1 => W+ - // u + dbar -> 2 + (-1) = 1 => W+ - } - if(iwch>0) { - (*idup)[iwup-1]=24; - (*idup)[ibup-1]=5; - (*mothup1)[iwup-1]=it; - (*mothup2)[iwup-1]=0; - (*mothup1)[ibup-1]=it; - (*mothup2)[ibup-1]=0; - } else if (iwch<0) { - (*idup)[iwup-1]=-24; - (*idup)[ibup-1]=-5; - (*mothup1)[iwup-1]=itb; - (*mothup2)[iwup-1]=0; - (*mothup1)[ibup-1]=itb; - (*mothup2)[ibup-1]=0; - } - (*istup)[iwup-1]=2; // The W is an intermediate and the - (*istup)[ibup-1]=1; // b is a final-state particle. - // Now reconstruct W boson momentum - double tmp=(*pup)[3][iwdec-1]+(*pup)[3][iwdec]; - (*pup)[3][iwup-1]=tmp; // W energy. - tmp=sqr(tmp); - for(unsigned int ixx=0; ixx<=2; ixx++) { - (*pup)[ixx][iwup-1] = (*pup)[ixx][iwdec-1] // W's 3-mom. - + (*pup)[ixx][iwdec]; - tmp=tmp-sqr((*pup)[ixx][iwup-1]); // Equals m^2 at the end of the loop. - } - (*pup)[4][iwup-1]=sqrt(tmp); // W mass. - // Reconstruct b momentum - tmp=(*pup)[3][nInt((*mothup1)[iwup-1])-1]-(*pup)[3][iwup-1]; - (*pup)[3][ibup-1]=tmp; // b energy. - tmp=sqr(tmp); - for(unsigned int ixx=0; ixx<=2; ixx++) { - (*pup)[ixx][ibup-1] = (*pup)[ixx][nInt((*mothup1)[iwup-1])-1] // b's 3mom. - - (*pup)[ixx][iwup-1]; - tmp=tmp-sqr((*pup)[ixx][ibup-1]); // Equals m^2 at the end of the loop. - } - (*pup)[4][ibup-1]=sqrt(tmp); // b mass. - (*icolup1)[iwup-1]=0; // W has no colour - (*icolup2)[iwup-1]=0; // lines. - (*icolup1)[ibup-1]=(*icolup1)[nInt((*mothup1)[iwup-1])-1]; // b shares top - (*icolup2)[ibup-1]=(*icolup2)[nInt((*mothup1)[iwup-1])-1]; // colour line. - } - (*nup)+=4; - } - // Single top production ( top*, ihrd=13): - else if (ihrd==13) { - int nw=1; - if(itopprc>=3) nw=2; - // Assign a mass to the incoming bottom quark, if there is one, - // rescaling the energy to accommodate the mass. - for(unsigned int ixx=1; ixx<=2; ixx++) { - if(abs((*idup)[ixx-1])==5) { - (*pup)[4][ixx-1]=masses[1]; - (*pup)[3][ixx-1]=sqrt(sqr((*pup)[2][ixx-1])+sqr((*pup)[4][ixx-1])); - } - } - // Set the top status to that of an intermediate. - (*istup)[2]=2; - // Get the index of the t / tbar in the evt. record (LH convention: 1->nup). - unsigned int it=0; - unsigned itb=0; - if((*idup)[2]==6) it=3; - else if((*idup)[2]==-6) itb=3; - else { - cout << "doIndividualHardProcessAssignments:\n" - << "wrong assumption about top position.\n" - << "Quitting ..."; - exit(1); - } - // Reconstruct intermediate W's from the decay products. - // iwdec is the index of the first W decay product - unsigned int iwdec(nInt(*nup)-1); // LH conventions: 1->nup. - if(nw==2) iwdec=nInt(*nup)-3; - // The W and b will go at the end of the record. - unsigned int iwup(nInt(*nup)+1); - unsigned int ibup(nInt(*nup)+2); - int iwch(0); - for(unsigned int iup=iwdec; iup<=iwdec+1; iup++) { - (*mothup1)[iup-1]=iwup; - (*mothup2)[iup-1]=0; - iwch=iwch-int((*idup)[iup-1])%2; // iwch = charge of W boson. - // electron+nubar -> 11 + (-12) = -1 => W- - // d + ubar -> 1 + (-2) = -1 => W- - // positron+nu -> -11+ 12 = 1 => W+ - // u + dbar -> 2 + (-1) = 1 => W+ - } - for(unsigned int ixx=0; ixx<2;ixx++) { - (*idup).push_back(0); - (*istup).push_back(0); - (*mothup1).push_back(0); - (*mothup2).push_back(0); - (*icolup1).push_back(0); - (*icolup2).push_back(0); - (*pup).push_back(vector(5)); - } - if(iwch>0) { - (*idup)[iwup-1]=24; - (*idup)[ibup-1]=5; - (*mothup1)[iwup-1]=it; - (*mothup2)[iwup-1]=0; - (*mothup1)[ibup-1]=it; - (*mothup2)[ibup-1]=0; - } else if (iwch<0) { - (*idup)[iwup-1]=-24; - (*idup)[ibup-1]=-5; - (*mothup1)[iwup-1]=itb; - (*mothup2)[iwup-1]=0; - (*mothup1)[ibup-1]=itb; - (*mothup2)[ibup-1]=0; - } - (*istup)[iwup-1]=2; - (*istup)[ibup-1]=1; - // Now reconstruct W boson momentum - double tmp=(*pup)[3][iwdec-1]+(*pup)[3][iwdec]; - (*pup)[3][iwup-1]=tmp; // W energy. - tmp=sqr(tmp); - for(unsigned int ixx=0; ixx<=2; ixx++) { - (*pup)[ixx][iwup-1] = (*pup)[ixx][iwdec-1] // W's 3-mom. - + (*pup)[ixx][iwdec]; - tmp=tmp-sqr((*pup)[ixx][iwup-1]); // Equals m^2 at the end of the loop. - } - (*pup)[4][iwup-1]=sqrt(tmp); // W mass. - // Reconstruct b momentum - tmp=(*pup)[3][nInt((*mothup1)[iwup-1])-1]-(*pup)[3][iwup-1]; - (*pup)[3][ibup-1]=tmp; // b energy. - tmp=sqr(tmp); - for(unsigned int ixx=0; ixx<=2; ixx++) { - (*pup)[ixx][ibup-1] = (*pup)[ixx][nInt((*mothup1)[iwup-1])-1] // b's 3mom. - - (*pup)[ixx][iwup-1]; - tmp=tmp-sqr((*pup)[ixx][ibup-1]); // Equals m^2 at the end of the loop. - } - (*pup)[4][ibup-1]=sqrt(tmp); // b mass. - (*icolup1)[iwup-1]=0; - (*icolup2)[iwup-1]=0; - (*icolup1)[ibup-1]=(*icolup1)[nInt((*mothup1)[iwup-1])-1]; - (*icolup2)[ibup-1]=(*icolup2)[nInt((*mothup1)[iwup-1])-1]; - (*nup)=(*nup)+2; - if(nw==2) { - // W DECAY - // iwdec is the index of the first W decay product. - iwdec=nInt(*nup)-3; - // iwup is the index of the W in the event record (LH conventions: 1->nup). - iwup=nInt(*nup)-6; - iwch=0; - int iwch(0); - for(unsigned int iup=iwdec; iup<=iwdec+1; iup++) { - (*mothup1)[iup-1]=iwup; - (*mothup2)[iup-1]=0; - iwch=iwch-int((*idup)[iup-1])%2; // iwch = charge of W boson. - // electron+nubar -> 11 + (-12) = -1 => W- - // d + ubar -> 1 + (-2) = -1 => W- - // positron+nu -> -11+ 12 = 1 => W+ - // u + dbar -> 2 + (-1) = 1 => W+ - } - (*istup)[iwup-1]=2; - (*icolup1)[iwup-1]=0; - (*icolup2)[iwup-1]=0; - } - } - - return; -} diff --git a/Contrib/AlpGen/BasicLesHouchesFileReader.cc b/Contrib/AlpGen/BasicLesHouchesFileReader.cc deleted file mode 100644 --- a/Contrib/AlpGen/BasicLesHouchesFileReader.cc +++ /dev/null @@ -1,466 +0,0 @@ -// -*- C++ -*- -// -// BasicLesHouchesFileReader.cc is a part of Herwig - A multi-purpose -// Monte Carlo event generator. -// Copyright (C) 2002-2017 The Herwig Collaboration -// -// Herwig is licenced under version 3 of the GPL, see COPYING for details. -// Please respect the MCnet academic guidelines, see GUIDELINES for details. -// -// -// This is the implementation of the non-inlined, non-templated member -// functions of the BasicLesHouchesFileReader class. -// -#include "BasicLesHouchesFileReader.h" -#include "ThePEG/Interface/ClassDocumentation.h" -#include "ThePEG/Interface/Reference.h" -#include "ThePEG/Interface/Switch.h" -#include "ThePEG/Interface/Parameter.h" -#include "ThePEG/Utilities/Throw.h" -#include "ThePEG/PDT/DecayMode.h" -#include "ThePEG/Persistency/PersistentOStream.h" -#include "ThePEG/Persistency/PersistentIStream.h" -#include "ThePEG/PDF/PartonExtractor.h" -#include "ThePEG/PDF/NoPDF.h" -#include "ThePEG/Cuts/Cuts.h" -#include "ThePEG/EventRecord/TmpTransform.h" -#include "ThePEG/Utilities/UtilityBase.h" - -using namespace Herwig; - -BasicLesHouchesFileReader:: -BasicLesHouchesFileReader(const BasicLesHouchesFileReader & x) - : LesHouchesReader(x), neve(x.neve), ieve(0), - LHFVersion(x.LHFVersion), outsideBlock(x.outsideBlock), - headerBlock(x.headerBlock), initComments(x.initComments), - initAttributes(x.initAttributes), eventComments(x.eventComments), - eventAttributes(x.eventAttributes), - theFileName(x.theFileName),overSampling_(x.overSampling_) {} - -BasicLesHouchesFileReader::~BasicLesHouchesFileReader() {} - -IBPtr BasicLesHouchesFileReader::clone() const { - return new_ptr(*this); -} - -IBPtr BasicLesHouchesFileReader::fullclone() const { - return new_ptr(*this); -} - -bool BasicLesHouchesFileReader::preInitialize() const { - return true; -} - -void BasicLesHouchesFileReader::doinit() { - LesHouchesReader::doinit(); -} - -void BasicLesHouchesFileReader::initialize(LesHouchesEventHandler & eh) { - LesHouchesReader::initialize(eh); - if ( LHFVersion.empty() ) - Throw() - << "The file associated with '" << name() << "' does not contain a " - << "proper formatted Les Houches event file. The events may not be " - << "properly sampled." << Exception::warning; -} - - -long BasicLesHouchesFileReader::scan() { - - open(); - - // Shall we write the events to a cache file for fast reading? If so - // we write to a temporary file if the caches events should be - // randomized. - if ( cacheFileName().length() ) openWriteCacheFile(); - - // Keep track of the number of events scanned. - long neve = 0; - long cuteve = 0; - bool negw = false; - - // If the open() has not already gotten information about subprocesses - // and cross sections we have to scan through the events. - if ( !heprup.NPRUP || cacheFile() || abs(heprup.IDWTUP) != 1 ) { // why scan if IDWTUP != 1? - - HoldFlag<> isScanning(scanning); - - double oldsum = 0.0; - vector lprup; - vector newmax; - vector oldeve; - vector neweve; - for ( int i = 0; ( maxScan() < 0 || i < maxScan() ) && readEvent(); ++i ) { - if ( !checkPartonBin() ) Throw() - << "Found event in LesHouchesReader '" << name() - << "' which cannot be handeled by the assigned PartonExtractor '" - << partonExtractor()->name() << "'." << Exception::runerror; - vector::iterator idit = - find(lprup.begin(), lprup.end(), hepeup.IDPRUP); - int id = lprup.size(); - if ( idit == lprup.end() ) { - lprup.push_back(hepeup.IDPRUP); - newmax.push_back(0.0); - neweve.push_back(0); - oldeve.push_back(0); - } else { - id = idit - lprup.begin(); - } - ++neve; - ++oldeve[id]; - oldsum += hepeup.XWGTUP; - if ( cacheFile() ) { - if ( eventWeight() == 0.0 ) { - ++cuteve; - continue; - } - cacheEvent(); - } - ++neweve[id]; - newmax[id] = max(newmax[id], abs(eventWeight())); - if ( eventWeight() < 0.0 ) negw = true; - } - - xSecWeights.resize(oldeve.size(), 1.0); - for ( int i = 0, N = oldeve.size(); i < N; ++i ) - if ( oldeve[i] ) xSecWeights[i] = double(neweve[i])/double(oldeve[i]); - - if ( maxScan() < 0 || neve > NEvents() ) NEvents(neve - cuteve); - - if ( lprup.size() == heprup.LPRUP.size() ) { - for ( int id = 0, N = lprup.size(); id < N; ++id ) { - vector::iterator idit = - find(heprup.LPRUP.begin(), heprup.LPRUP.end(), hepeup.IDPRUP); - if ( idit == heprup.LPRUP.end() ) { - Throw() - << "When scanning events, the LesHouschesReader '" << name() - << "' found undeclared processes." << Exception::warning; - heprup.NPRUP = 0; - break; - } - int idh = idit - heprup.LPRUP.begin(); - heprup.XMAXUP[idh] = newmax[id]; - } - } - if ( heprup.NPRUP == 0 ) { - // No heprup block was supplied or something went wrong. - heprup.NPRUP = lprup.size(); - heprup.LPRUP.resize(lprup.size()); - heprup.XMAXUP.resize(lprup.size()); - for ( int id = 0, N = lprup.size(); id < N; ++id ) { - heprup.LPRUP[id] = lprup[id]; - heprup.XMAXUP[id] = newmax[id]; - } - } else if ( abs(heprup.IDWTUP) != 1 ) { - // Try to fix things if abs(heprup.IDWTUP) != 1. - double sumxsec = 0.0; - for ( int id = 0; id < heprup.NPRUP; ++id ) sumxsec += heprup.XSECUP[id]; - weightScale = picobarn*neve*sumxsec/oldsum; - } - } - - if ( cacheFile() ) closeCacheFile(); - - if ( negw ) heprup.IDWTUP = min(-abs(heprup.IDWTUP), -1); - - return neve; - -} - -void BasicLesHouchesFileReader::open() { - if ( filename().empty() ) - throw LesHouchesFileError() - << "No Les Houches file name. " - << "Use 'set " << name() << ":FileName'." - << Exception::runerror; - cfile.open(filename()); - if ( !cfile ) - throw LesHouchesFileError() - << "The BasicLesHouchesFileReader '" << name() << "' could not open the " - << "event file called '" << theFileName << "'." - << Exception::runerror; - - cfile.readline(); - if ( !cfile.find(" attributes = - StringUtils::xmlAttributes("LesHouchesEvents", cfile.getline()); - LHFVersion = attributes["version"]; - if ( LHFVersion.empty() ) return; - - bool readingHeader = false; - bool readingInit = false; - headerBlock = ""; - - // Loop over all lines until we hit the tag. - while ( cfile.readline() && !cfile.find("") ) { - if ( cfile.find("> heprup.IDBMUP.first >> heprup.IDBMUP.second - >> heprup.EBMUP.first >> heprup.EBMUP.second - >> heprup.PDFGUP.first >> heprup.PDFGUP.second - >> heprup.PDFSUP.first >> heprup.PDFSUP.second - >> heprup.IDWTUP >> heprup.NPRUP ) ) { - heprup.NPRUP = -42; - LHFVersion = ""; - return; - } - heprup.resize(); - - for ( int i = 0; i < heprup.NPRUP; ++i ) { - cfile.readline(); - if ( !( cfile >> heprup.XSECUP[i] >> heprup.XERRUP[i] - >> heprup.XMAXUP[i] >> heprup.LPRUP[i] ) ) { - heprup.NPRUP = -42; - LHFVersion = ""; - return; - } - } - } - else if ( cfile.find("momentum().plus()/ - beams().first->momentum().plus(); - - if ( reweightPDF && - inPDF.first && outPDF.first && inPDF.first != outPDF.first ) { - if ( hepeup.XPDWUP.first <= 0.0 ) - hepeup.XPDWUP.first = - inPDF.first->xfx(inData.first, incoming().first->dataPtr(), - sqr(hepeup.SCALUP*GeV), x1); - double xf = outPDF.first->xfx(inData.first, incoming().first->dataPtr(), - sqr(hepeup.SCALUP*GeV), x1); - lastweight *= xf/hepeup.XPDWUP.first; - hepeup.XPDWUP.first = xf; - } - - double x2 = incoming().second->momentum().minus()/ - beams().second->momentum().minus(); - - if ( reweightPDF && - inPDF.second && outPDF.second && inPDF.second != outPDF.second ) { - if ( hepeup.XPDWUP.second <= 0.0 ) - hepeup.XPDWUP.second = - inPDF.second->xfx(inData.second, incoming().second->dataPtr(), - sqr(hepeup.SCALUP*GeV), x2); - double xf = - outPDF.second->xfx(inData.second, incoming().second->dataPtr(), - sqr(hepeup.SCALUP*GeV), x2); - lastweight *= xf/hepeup.XPDWUP.second; - hepeup.XPDWUP.second = xf; - } - - if ( cutEarly() ) { - if ( !cuts().initSubProcess((incoming().first->momentum() + - incoming().second->momentum()).m2(), - 0.5*log(x1/x2)) ) lastweight = 0.0; - tSubProPtr sub = getSubProcess(); - TmpTransform tmp(sub, Utilities::getBoostToCM(sub->incoming())); - if ( !cuts().passCuts(*sub) ) lastweight = 0.0; - } - - return true; -} - -double BasicLesHouchesFileReader::getEvent() { - if ( cacheFile() ) { - if (overSampling_) { - if ( !uncacheEvent() ) reopen(); - } else { - if ( !uncacheEvent() || stats.attempts()==NEvents() ) - throw LesHouchesReopenWarning() - << "More events requested than available in LesHouchesReader " - << name() << Exception::runerror; - } - } else { - if (overSampling_) { - if ( !readEvent() ) reopen(); - } else { - if ( !readEvent() || stats.attempts()==NEvents() ) - throw LesHouchesReopenWarning() - << "More events requested than available in LesHouchesReader " - << name() << Exception::runerror; - } - } - ++position; - - double max = maxWeights[hepeup.IDPRUP]*maxFactor; - return max != 0.0? eventWeight()/max: 0.0; - -} - -void BasicLesHouchesFileReader::skip(long n) { - HoldFlag<> skipflag(skipping); - if(overSampling_) while ( n-- ) getEvent(); -} - -bool BasicLesHouchesFileReader::doReadEvent() { - if ( !cfile ) return false; - if ( LHFVersion.empty() ) return false; - if ( heprup.NPRUP < 0 ) return false; - eventComments = ""; - outsideBlock = ""; - hepeup.NUP = 0; - hepeup.XPDWUP.first = hepeup.XPDWUP.second = 0.0; - - // Keep reading lines until we hit the next event or the end of - // the event block. Save any inbetween lines. Exit if we didn't - // find an event. - while ( cfile.readline() && !cfile.find("> hepeup.NUP >> hepeup.IDPRUP >> hepeup.XWGTUP - >> hepeup.SCALUP >> hepeup.AQEDUP >> hepeup.AQCDUP ) ) - return false; - hepeup.resize(); - // Read all particle lines. - for ( int i = 0; i < hepeup.NUP; ++i ) { - if ( !cfile.readline() ) return false; - if ( !( cfile >> hepeup.IDUP[i] >> hepeup.ISTUP[i] - >> hepeup.MOTHUP[i].first >> hepeup.MOTHUP[i].second - >> hepeup.ICOLUP[i].first >> hepeup.ICOLUP[i].second - >> hepeup.PUP[i][0] >> hepeup.PUP[i][1] >> hepeup.PUP[i][2] - >> hepeup.PUP[i][3] >> hepeup.PUP[i][4] - >> hepeup.VTIMUP[i] >> hepeup.SPINUP[i] ) ) - return false; - if(std::isnan(hepeup.PUP[i][0])||std::isnan(hepeup.PUP[i][1])|| - std::isnan(hepeup.PUP[i][2])||std::isnan(hepeup.PUP[i][3])|| - std::isnan(hepeup.PUP[i][4])) - throw Exception() - << "nan's as momenta in Les Houches file " - << Exception::eventerror; - } - - // Now read any additional comments. - while ( cfile.readline() && !cfile.find("") ) - eventComments += cfile.getline() + "\n"; - - if ( !cfile ) return false; - return true; - -} - -void BasicLesHouchesFileReader::close() { - cfile.close(); -} - -void BasicLesHouchesFileReader::persistentOutput(PersistentOStream & os) const { - os << neve << LHFVersion << outsideBlock << headerBlock << initComments - << initAttributes << eventComments << eventAttributes << theFileName - << overSampling_; -} - -void BasicLesHouchesFileReader::persistentInput(PersistentIStream & is, int) { - is >> neve >> LHFVersion >> outsideBlock >> headerBlock >> initComments - >> initAttributes >> eventComments >> eventAttributes >> theFileName - >> overSampling_; - ieve = 0; -} - -ClassDescription -BasicLesHouchesFileReader::initBasicLesHouchesFileReader; -// Definition of the static class description member. - -void BasicLesHouchesFileReader::Init() { - - static ClassDocumentation documentation - ("Herwig::BasicLesHouchesFileReader is an base class to be used for objects " - "which reads event files from matrix element generators. This class is " - "able to read plain event files conforming to the Les Houches Event File " - "accord."); - - static Parameter interfaceFileName - ("FileName", - "The name of a file containing events conforming to the Les Houches " - "protocol to be read into ThePEG. A file name ending in " - ".gz will be read from a pipe which uses " - "zcat. If a file name ends in | the " - "preceeding string is interpreted as a command, the output of which " - "will be read through a pipe.", - &BasicLesHouchesFileReader::theFileName, "", false, false); - - static Switch interfaceOverSampling - ("OverSampling", - "Allow / Forbid reading of LH events more than once by the " - "LH reader, allowing / protecting against statistical problems.", - &BasicLesHouchesFileReader::overSampling_, true, false, false); - static SwitchOption AllowOverSampling - (interfaceOverSampling, - "AllowOverSampling", - "The reader will read events in the file more than once if more " - "events are needed to generate the requested number than that in " - "the LH file.", - true); - static SwitchOption ForbidOverSampling - (interfaceOverSampling, - "ForbidOverSampling", - "The reader will NOT read events in the file more than once if more " - "events are needed to generate the requested number than that in " - "the LH file - instead it will stop when all have been read.", - false); - - interfaceFileName.fileType(); - interfaceFileName.rank(11); - -} - diff --git a/Contrib/AlpGen/BasicLesHouchesFileReader.fh b/Contrib/AlpGen/BasicLesHouchesFileReader.fh deleted file mode 100644 --- a/Contrib/AlpGen/BasicLesHouchesFileReader.fh +++ /dev/null @@ -1,22 +0,0 @@ -// -*- C++ -*- -// -// This is the forward declaration of the BasicLesHouchesFileReader class. -// -#ifndef HERWIG_BasicLesHouchesFileReader_FH -#define HERWIG_BasicLesHouchesFileReader_FH - -#include "ThePEG/Config/Pointers.h" - -namespace Herwig { - -class BasicLesHouchesFileReader; - -} - -namespace ThePEG { - -ThePEG_DECLARE_POINTERS(Herwig::BasicLesHouchesFileReader,BasicLesHouchesFileReaderPtr); - -} - -#endif diff --git a/Contrib/AlpGen/BasicLesHouchesFileReader.h b/Contrib/AlpGen/BasicLesHouchesFileReader.h deleted file mode 100644 --- a/Contrib/AlpGen/BasicLesHouchesFileReader.h +++ /dev/null @@ -1,333 +0,0 @@ -// -*- C++ -*- -// -// BasicLesHouchesFileReader.h is a part of Herwig - A multi-purpose -// Monte Carlo event generator. -// Copyright (C) 2002-2017 The Herwig Collaboration -// -// Herwig is licenced under version 3 of the GPL, see COPYING for details. -// Please respect the MCnet academic guidelines, see GUIDELINES for details. -// -// -// This is the implementation of the non-inlined, non-templated member -// functions of the BasicLesHouchesFileReader class. -// -#ifndef HERWIG_BasicLesHouchesFileReader_H -#define HERWIG_BasicLesHouchesFileReader_H -// This is the declaration of the BasicLesHouchesFileReader class. - -#include "ThePEG/LesHouches/LesHouchesReader.h" -#include "BasicLesHouchesFileReader.fh" -#include "ThePEG/PDT/Decayer.h" -#include "ThePEG/Utilities/CFileLineReader.h" -#include - -namespace Herwig { - -using namespace ThePEG; - -/** - * BasicLesHouchesFileReader derives from the LesHouchesReader base class - * to be used for objects which read event files from matrix element - * generators. It extends LesHouchesReader by defining a file handle to be - * read from, which is opened and closed by the open() and close() - * functions. Note that the file handle is a standard C filehandle and - * not a C++ stream. This is because there is no standard way in C++ - * to connect a pipe to a stream for reading eg. gzipped files. This - * class is able to read plain event files conforming to the Les - * Houches Event File accord. - * - * @see \ref BasicLesHouchesFileReaderInterfaces "The interfaces" - * defined for BasicLesHouchesFileReader. - * @see LesHouchesReader - * @see BasicLesHouchesFileReader - */ -class BasicLesHouchesFileReader: public LesHouchesReader { - -public: - - /** @name Standard constructors and destructors. */ - //@{ - /** - * Default constructor. - */ - BasicLesHouchesFileReader() : neve(0), ieve(0) {} - - /** - * Copy-constructor. Note that a file which is opened in the object - * copied from will have to be reopened in this. - */ - BasicLesHouchesFileReader(const BasicLesHouchesFileReader &); - - /** - * Destructor. - */ - virtual ~BasicLesHouchesFileReader(); - //@} - -public: - - /** @name Virtual functions specified by the LesHouchesReader base class. */ - //@{ - /** - * Initialize. This function is called by the LesHouchesEventHandler - * to which this object is assigned. - */ - virtual void initialize(LesHouchesEventHandler & eh); - - /** - * Calls readEvent() or uncacheEvent() to read information into the - * LesHouches common block variables. This function is called by the - * LesHouchesEventHandler if this reader has been selectod to - * produce an event. - * - * @return the weight asociated with this event. If negative weights - * are allowed it should be between -1 and 1, otherwise between 0 - * and 1. If outside these limits the previously estimated maximum - * is violated. Note that the estimated maximum then should be - * updated from the outside. - */ - virtual double getEvent(); - - /** - * Calls doReadEvent() and performs pre-defined reweightings. A - * sub-class overrides this function it must make sure that the - * corresponding reweightings are done. - */ - virtual bool readEvent(); - - /** - * Skip \a n events. Used by LesHouchesEventHandler to make sure - * that a file is scanned an even number of times in case the events - * are not ramdomly distributed in the file. - */ - virtual void skip(long n); - - /** - * Scan the file or stream to obtain information about cross section - * weights and particles etc. This function should fill the - * variables corresponding to the /HEPRUP/ common block. The - * function returns the number of events scanned. - */ - virtual long scan(); - - /** - * Open a file with events. Derived classes should overwrite it and - * first calling it before reading in the run information into the - * corresponding protected variables. - */ - virtual void open(); - - /** - * Close the file from which events have been read. - */ - virtual void close(); - - /** - * Read the next event from the file or stream into the - * corresponding protected variables. Return false if there is no - * more events or if this was not a LHF event file. - */ - virtual bool doReadEvent(); - //@} - - /** - * Return the name of the file from where to read events. - */ - string filename() const { return theFileName; } - -public: - - /** @name Functions used by the persistent I/O system. */ - //@{ - /** - * Function used to write out object persistently. - * @param os the persistent output stream written to. - */ - void persistentOutput(PersistentOStream & os) const; - - /** - * Function used to read in object persistently. - * @param is the persistent input stream read from. - * @param version the version number of the object when written. - */ - void persistentInput(PersistentIStream & is, int version); - //@} - - /** - * Standard Init function used to initialize the interfaces. - */ - static void Init(); - -protected: - - /** @name Clone Methods. */ - //@{ - /** - * Make a simple clone of this object. - * @return a pointer to the new object. - */ - virtual IBPtr clone() const; - - /** Make a clone of this object, possibly modifying the cloned object - * to make it sane. - * @return a pointer to the new object. - */ - virtual IBPtr fullclone() const; - //@} - - /** @name Standard (and non-standard) Interfaced functions. */ - //@{ - /** - * Initialize this object after the setup phase before saving an - * EventGenerator to disk. - * @throws InitException if object could not be initialized properly. - */ - virtual void doinit(); - - /** - * Return true if this object needs to be initialized before all - * other objects because it needs to extract PDFs from the event file. - */ - virtual bool preInitialize() const; - //@ - -protected: - - /** - * The wrapper around the C FILE stream from which to read - */ - CFileLineReader cfile; - -protected: - - /** - * The number of events in this file. - */ - long neve; - - /** - * The current event number. - */ - long ieve; - - /** - * If the file is a standard Les Houches formatted file (LHF) this - * is its version number. If empty, this is not a Les Houches - * formatted file - */ - string LHFVersion; - - /** - * If LHF. All lines (since the last open() or readEvent()) outside - * the header, init and event tags. - */ - string outsideBlock; - - /** - * If LHF. All lines from the header block. - */ - string headerBlock; - - /** - * If LHF. Additional comments found in the init block. - */ - string initComments; - - /** - * If LHF. Map of attributes (name-value pairs) found in the init - * tag. - */ - map initAttributes; - - /** - * If LHF. Additional comments found with the last read event. - */ - string eventComments; - - /** - * If LHF. Map of attributes (name-value pairs) found in the last - * event tag. - */ - map eventAttributes; - -private: - - /** - * The name of the file from where to read events. - */ - string theFileName; - - /** - * Determines whether events in the LH file are or are not read - * more than once in order to generate the requested number of - * events. - */ - bool overSampling_; - -private: - - /** - * Describe an abstract base class with persistent data. - */ - static ClassDescription initBasicLesHouchesFileReader; - - /** - * Private and non-existent assignment operator. - */ - BasicLesHouchesFileReader & operator=(const BasicLesHouchesFileReader &) = delete; - -public: - - /** @cond EXCEPTIONCLASSES */ - /** Exception class used by BasicLesHouchesFileReader if reading the file - * fails. */ - class LesHouchesFileError: public Exception {}; - /** @endcond */ - -}; - -} - - -#include "ThePEG/Utilities/ClassTraits.h" - -namespace ThePEG { - -/** @cond TRAITSPECIALIZATIONS */ - -/** - * This template specialization informs ThePEG about the - * base class of BasicLesHouchesFileReader. - */ -template <> -struct BaseClassTrait { - /** Typedef of the base class of BasicLesHouchesFileReader. */ - typedef LesHouchesReader NthBase; -}; - -/** - * This template specialization informs ThePEG about the name of the - * BasicLesHouchesFileReader class and the shared object where it is - * defined. - */ -template <> -struct ClassTraits - : public ClassTraitsBase { - /** - * Return the class name. - */ - static string className() { return "Herwig::BasicLesHouchesFileReader"; } - /** - * Return the name of the shared library to be loaded to get access - * to the BasicLesHouchesFileReader class and every other class it uses - * (except the base class). - */ - static string library() { return "BasicLesHouchesFileReader.so"; } - -}; - -/** @endcond */ - -} - -#endif /* HERWIG_BasicLesHouchesFileReader_H */ diff --git a/Contrib/AlpGen/Makefile.in b/Contrib/AlpGen/Makefile.in deleted file mode 100644 --- a/Contrib/AlpGen/Makefile.in +++ /dev/null @@ -1,54 +0,0 @@ -# -*- Makefile -*- (for emacs) -# -# This Makefile is intended for compiling Herwig++ plugins -# -# This Makefile received very little testing, -# any bug reports are very welcome! -# -# Location of include files -THEPEGINCLUDE = -HERWIGINCLUDE = -GSLINCLUDE = -FASTJETINCLUDE = -# Messy hack, not guaranteed to work: -FASTJETLIB = `echo $(FASTJETINCLUDE) | sed "s/-I//" | sed "s/ //"`/../lib/ -LDFLAGS = -L$(FASTJETLIB) -lfastjet -SHARED_FLAG = -INCLUDE = $(THEPEGINCLUDE) $(HERWIGINCLUDE) $(GSLINCLUDE) $(FASTJETINCLUDE) -# C++ flags -CXX = -CXXFLAGS = - -ALLCCFILES=$(shell echo *.cc) - -default : AlpGenHandler.so AlpGenToLH.exe - -%.o : %.cc %.h - $(CXX) -fPIC $(CPPFLAGS) $(INCLUDE) $(CXXFLAGS) -c -shared $< -o $@ - -BasicLesHouchesFileReader.so : BasicLesHouchesFileReader.o - $(CXX) -fPIC $(CPPFLAGS) $(INCLUDE) $(CXXFLAGS) $(SHARED_FLAG) $(LDFLAGS) \ - -Wl,-undefined -Wl,dynamic_lookup \ - BasicLesHouchesFileReader.o -o BasicLesHouchesFileReader.so - -AlpGenHandler.so : AlpGenHandler.o - $(CXX) -fPIC $(CPPFLAGS) $(INCLUDE) $(CXXFLAGS) $(SHARED_FLAG) $(LDFLAGS) \ - -Wl,-undefined -Wl,dynamic_lookup \ - AlpGenHandler.o -o AlpGenHandler.so - -AlpGenToLH.exe : AlpGenToLH.o - $(CXX) AlpGenToLH.o -o AlpGenToLH.exe - -MLMHistogram.so : Validation/MLMHistogram.o - $(CXX) -fPIC $(CPPFLAGS) $(INCLUDE) $(CXXFLAGS) $(SHARED_FLAG) $(LDFLAGS) \ - -Wl,-undefined -Wl,dynamic_lookup \ - Validation/MLMHistogram.o -o MLMHistogram.so - -MLMAnalysis.so : Validation/MLMAnalysis.o - $(CXX) -fPIC $(CPPFLAGS) $(INCLUDE) $(CXXFLAGS) $(SHARED_FLAG) $(LDFLAGS) \ - -Wl,-undefined -Wl,dynamic_lookup \ - MLMHistogram.so Validation/MLMAnalysis.o -o MLMAnalysis.so - -clean: - rm -f $(ALLCCFILES:.cc=.o) \ - BasicLesHouchesFileReader.so AlpGenHandler.so AlpGenToLH.exe MLMAnalysis.so MLMHistogram.so \ No newline at end of file