diff --git a/Config/RopeUserHooks.h b/Config/RopeUserHooks.h --- a/Config/RopeUserHooks.h +++ b/Config/RopeUserHooks.h @@ -1,371 +1,371 @@ #ifndef THEP8I_RopeUserHooks_H #define THEP8I_RopeUserHooks_H #include "Pythia8/Pythia.h" #include "ThePEG/EventRecord/Particle.h" #include "TheP8I/Hadronization/ParameterHandler.h" #include "TheP8I/Hadronization/Ropewalk.h" #include "ThePEG/Utilities/Throw.h" namespace TheP8I { class RopeUserHooks : public Pythia8::UserHooks { // Convenient typedefs typedef map PytPars; typedef map DipMass; typedef map FlavourEnd; public: RopeUserHooks() : dPtrSave(NULL), _ph(NULL), _h(-1.0), _m0(0*GeV), _pTcut(100*GeV), _r0(1*femtometer), _window(false) { } ~RopeUserHooks() { } virtual bool canChangeFragPar() { return true; } virtual bool doChangeFragPar(Pythia8::StringFlav* flavPtr, Pythia8::StringZ* zPtr, Pythia8::StringPT * pTPtr, int endFlavour, double m2HadIn, vector iParton) { // We're not using it here (void) iParton; // Get new parameters m2Had = m2HadIn; PytPars newPar = fetchParameters(endFlavour); if(newPar.find("null") != newPar.end()) { Throw() << "Problem fetching parameters in RopeUserHook. " << "Ropes switched off in this string." << Exception::warning; _h = 1.0; newPar = fetchParameters(endFlavour); _h = -1.0; } for(PytPars::iterator itr = newPar.begin(); itr!=newPar.end();++itr) settingsPtr->parm(itr->first,itr->second); // Re-initialize all three - flavPtr->init(*settingsPtr,rndmPtr); - zPtr->init(*settingsPtr,*particleDataPtr,rndmPtr); - pTPtr->init(*settingsPtr,*particleDataPtr,rndmPtr); + flavPtr->init(*settingsPtr,particleDataPtr,rndmPtr,infoPtr); + zPtr->init(*settingsPtr,*particleDataPtr,rndmPtr,infoPtr); + pTPtr->init(*settingsPtr,particleDataPtr,rndmPtr, infoPtr); //settingsPtr->listChanged(); return true; } /* virtual Pythia8::Vec4 shoveTheHadron(Pythia8::Vec4 pHadIn){ if(!dPtrSave) return pHadIn; Lorentz5Momentum pH(pHadIn.px()*GeV, pHadIn.py()*GeV, pHadIn.pz()*GeV, pHadIn.e()*GeV); LorentzMomentum pOut = dPtrSave->shoveOverlappers(dipFrac,m2Had, _r0, pH); return Pythia8::Vec4(pOut.x()/GeV, pOut.y()/GeV, pOut.z()/GeV, pOut.t()/GeV); } */ void setParameterHandler(ParameterHandler * ph){ _ph = ph; } void setEnhancement(double h){ // Will overrule others if set. _h = h; } void setWindow(bool window, Energy pTcut){ // Can set a window where rope should not apply _window = window; _pTcut = pTcut; } bool setDipoles(const pair > * dm, Energy m0, Length r0, bool throwHadr = false, double alpha = 1.0){ // If we set the dipoles, we should also use them. _h = -1.0; _m0 = m0; _r0 = r0; _throwHadr = throwHadr; _alpha = alpha; vector dip = dm->second; // Closed gluon loop -- we will return false and let StringFragmentation figure out what to do if(dip.front()->pa->id() == 21 && dip.front()->pc->id() == 21 && dip.back()->pa->id() == 21 && dip.back()->pc->id() == 21) return false; // Get flavour of quark ends -- get the one that is not a gluon. // We always assign left to the front of the // dipole vector and right to the back. DipMass leftC; DipMass rightC; leftC.clear(); rightC.clear(); if(dip.size()==1){ _leftP = dip[0]->pc; _rightP = dip[0]->pa; } else{ _leftP = (dip.front()->pa->id() != 21 ? dip.front()->pa : dip.front()->pc); _rightP = (dip.back()->pc->id() != 21 ? dip.back()->pc : dip.back()->pa); } if (_leftP->id() == _rightP->id() ) Throw() << "Flavours in strings corrupted. " << "This is a serious error - please contact the authors." << Exception::abortnow; // A dipole gets all of quark mass, half of gluon // The end quarks energy m0 // Start by adding half of the first quark // from both left and right LorentzMomentum momL(0*GeV,0*GeV,0*GeV,0*GeV); LorentzMomentum momR(0*GeV,0*GeV,0*GeV,0*GeV); for(size_t i = 0, N = dip.size(); i < N; ++i ){ Ropewalk::Dipole * dPtrL = dip[i]; Ropewalk::Dipole * dPtrR = dip[N - i - 1]; if(!dPtrL || dPtrL->broken || !dPtrR || dPtrR->broken){ Throw() << "Broken dipole in RopeUserHooks. " << "This is a serious error - please contact the authors." << Exception::abortnow; return false; } momL += 0.5*(dPtrL->pc->momentum()); momL += 0.5*(dPtrL->pa->momentum()); momR += 0.5*(dPtrR->pc->momentum()); momR += 0.5*(dPtrR->pa->momentum()); if(i == 0){ momL += 0.5 * _leftP->momentum(); momR += 0.5 * _rightP->momentum(); } if(i == N - 1){ momR += 0.5 * _leftP->momentum(); momL += 0.5 * _rightP->momentum(); } leftC.insert(DipMass::value_type(momL.m2(),dPtrL)); rightC.insert(DipMass::value_type(momR.m2(),dPtrR)); } flavourDipoles.clear(); flavourDipoles.insert( FlavourEnd::value_type(_leftP,leftC) ); flavourDipoles.insert( FlavourEnd::value_type(_rightP,rightC) ); return true; } bool setDipoles(const pair > * dm, Energy m0, Length r0, bool throwHadr = false, double alpha = 1.0){ // If we set the dipoles, we should also use them. _h = -1.0; _m0 = m0; _r0 = r0; _throwHadr = throwHadr; _alpha = alpha; vector dip = dm->second; // Closed gluon loop -- we will return false and let StringFragmentation figure out what to do if(dip.front()->pa->id() == 21 && dip.front()->pc->id() == 21 && dip.back()->pa->id() == 21 && dip.back()->pc->id() == 21) return false; // Get flavour of quark ends -- get the one that is not a gluon. // We always assign left to the front of the // dipole vector and right to the back. DipMass leftC; DipMass rightC; leftC.clear(); rightC.clear(); if(dip.size()==1){ _leftP = dip[0]->pc; _rightP = dip[0]->pa; } else{ _leftP = (dip.front()->pa->id() != 21 ? dip.front()->pa : dip.front()->pc); _rightP = (dip.back()->pc->id() != 21 ? dip.back()->pc : dip.back()->pa); } if (_leftP->id() == _rightP->id() ) Throw() << dip.size() << " " << dip.front()->pa->id() << " " << dip.front()->pc->id() << " " << dip.front()->pc->id() << " " << dip.front()->pa->id() << "\n" << "Flavours in strings corrupted (setDipoles). " << "This is a serious error - please contact the authors." << Exception::abortnow; // A dipole gets all of quark mass, half of gluon // The end quarks energy m0 // Start by adding half of the first quark // from both left and right LorentzMomentum momL(0*GeV,0*GeV,0*GeV,0*GeV); LorentzMomentum momR(0*GeV,0*GeV,0*GeV,0*GeV); for(size_t i = 0, N = dip.size(); i < N; ++i ){ Ropewalk::Dipole * dPtrL = dip[i]; Ropewalk::Dipole * dPtrR = dip[N - i - 1]; if(!dPtrL || dPtrL->broken || !dPtrR || dPtrR->broken){ Throw() << "Broken dipole in RopeUserHooks. " << "This is a serious error - please contact the authors." << Exception::abortnow; return false; } momL += 0.5*(dPtrL->pc->momentum()); momL += 0.5*(dPtrL->pa->momentum()); momR += 0.5*(dPtrR->pc->momentum()); momR += 0.5*(dPtrR->pa->momentum()); if(i == 0){ momL += 0.5 * _leftP->momentum(); momR += 0.5 * _rightP->momentum(); } if(i == N - 1){ momR += 0.5 * _leftP->momentum(); momL += 0.5 * _rightP->momentum(); } leftC.insert(DipMass::value_type(momL.m2(),dPtrL)); rightC.insert(DipMass::value_type(momR.m2(),dPtrR)); } flavourDipoles.clear(); flavourDipoles.insert( FlavourEnd::value_type(_leftP,leftC) ); flavourDipoles.insert( FlavourEnd::value_type(_rightP,rightC) ); return true; } public: /** * Exception class. */ struct RopeException: public Exception {}; private: PytPars fetchParameters(int endFlavour){ dPtrSave = NULL; // Did we remember to set the callback pointer? // Do not just construct a new one, as there will be no sensible initial values, resulting // in very subtle errors. if(!(_ph)){ Throw() << "Missing parameter handler in RopeUserHooks. " << "This is a serious error - please contact the authors." << Exception::abortnow; PytPars p; p.insert(make_pair("null",0.)); return p; } // If we have manually set enhancement, we should just use that, and avoid complicated behavior. if(_h > 0) return _ph->GetEffectiveParameters(_h); // Test the string ends... if( (endFlavour == _leftP->id() && endFlavour == _rightP->id()) || (endFlavour != _leftP->id() && endFlavour != _rightP->id()) ){ Throw() << "Flavours in strings corrupted (fetchParameters). " << "This is a serious error - please contact the authors." << Exception::abortnow; PytPars p; p.insert(make_pair("null",0.)); return p; } // Find the right end tcPPtr thisEnd = (_leftP->id() == endFlavour ? _leftP : _rightP); FlavourEnd::iterator feItr = flavourDipoles.find( thisEnd); if( feItr == flavourDipoles.end()){ Throw() << "Could not find dipole end. " << "Ropes switched off in this string." << Exception::warning; PytPars p; p.insert(make_pair("null",0.)); return p; } DipMass dms = feItr->second; // Find the first dipole where the total invariant mass sq. from the string // (left or right) exceeds m2Had DipMass::iterator dmsItr = dms.lower_bound(m2Had*GeV2); if( dmsItr == dms.end()){ Throw() << "Could not get correct dipole. mpythia: " << m2Had << " " << (--dmsItr)->first/GeV2 << ". " << "Ropes switched off in this string." << Exception::warning; PytPars p; p.insert(make_pair("null",0.)); return p; } Ropewalk::Dipole * dPtr = dmsItr->second; if(!dPtr){ Throw() << "Missing dipole in RopeUserHooks. " << "This is a serious error - please contact the authors." << Exception::abortnow; PytPars p; p.insert(make_pair("null",0.)); return p; } dPtrSave = dPtr; // Find out how long in we are on the dipole double mBig = dmsItr->first/GeV2; if(m2Had == 0) dipFrac = 0; else if( dmsItr == dms.begin() ) dipFrac = sqrt(m2Had/mBig); else{ double mSmall = double((--dmsItr)->first/GeV2); dipFrac = (sqrt(m2Had) - sqrt(mSmall)) / (sqrt(mBig) - sqrt(mSmall)); } // We need dipFrag to be fraction from the pc parton in the dipole. // Check if the starting end is pa or pc. If it is already a pc, do noting, // otherwise take fraction from other side if( thisEnd == dms.lower_bound(0*GeV2)->second->pa) dipFrac = 1 - dipFrac; // use the window here if(_window && ((dPtr->pc->momentum().perp() > _pTcut && dipFrac < 0.5) || (dPtr->pa->momentum().perp() > _pTcut && dipFrac >=0.5)) ){ return _ph->GetEffectiveParameters(1.0); } dPtr->reinit(dipFrac,_r0,_m0); // cout << "p = " << dPtr->p << " q = " << dPtr->q << " kappa " << dPtr->kappaEnhancement() << endl; if (_throwHadr){ (*dPtr).hadr = true; return _ph->GetEffectiveParameters(dPtr->firstKappa(_alpha)); } return _ph->GetEffectiveParameters(dPtr->kappaEnhancement()); } Ropewalk::Dipole* dPtrSave; tcPPtr _leftP, _rightP; FlavourEnd flavourDipoles; ParameterHandler * _ph; double _h, _alpha, dipFrac, m2Had; Energy _m0, _pTcut; Length _r0; bool _throwHadr, _window; }; } #endif diff --git a/Hadronization/Makefile.am b/Hadronization/Makefile.am --- a/Hadronization/Makefile.am +++ b/Hadronization/Makefile.am @@ -1,21 +1,19 @@ if PYTHIA8WORKS mySOURCES = StringFragmentation.cc TheP8EventShapes.cc StringPipe.cc OverlapPythiaHandler.cc ParameterHandler.cc RandomAverageHandler.cc RandomHandler.cc Ropewalk.cc HistKeeper.cc Histogram.cc Bin.cc DOCFILES = StringFragmentation.h TheP8EventShapes.h StringPipe.h OverlapPythiaHandler.h ParameterHandler.h RandomAverageHandler.h RandomHandler.h Ropewalk.h HistKeeper.h Histogram.h Bin.h AUTOGENERATED = StringFragmentation-var.h StringFragmentation-init.h \ StringFragmentation-input.h StringFragmentation-output.h \ StringFragmentation-interfaces.h INCLUDEFILES = $(DOCFILES) $(AUTOGENERATED) -AM_CPPFLAGS += -I/usr/local/new/include/ - noinst_LTLIBRARIES = libTheP8IString.la libTheP8IString_la_SOURCES = $(mySOURCES) $(INCLUDEFILES) libTheP8IString_la_CPPFLAGS = $(AM_CPPFLAGS) $(PYTHIA8_CPPFLAGS) endif include $(top_srcdir)/Config/Makefile.aminclude diff --git a/Hadronization/OverlapPythiaHandler.cc b/Hadronization/OverlapPythiaHandler.cc --- a/Hadronization/OverlapPythiaHandler.cc +++ b/Hadronization/OverlapPythiaHandler.cc @@ -1,216 +1,216 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the OverlapPythiaHandler class. // #include "OverlapPythiaHandler.h" using namespace TheP8I; OverlapPythiaHandler::OverlapPythiaHandler(){ } OverlapPythiaHandler::~OverlapPythiaHandler(){ for ( vector::iterator itr = _overlapPythias.begin(); itr != _overlapPythias.end(); ++itr ) itr->nullify(); _overlapPythias.clear(); } OverlapPythiaHandler::OverlapPythiaHandler(HadronizationHandler * sf, vectorarguments) : _sf(sf){ // Set the arguments for(vector::iterator itr = arguments.begin(); itr!=arguments.end();++itr){ if(itr->find("StringPT:sigma")!=string::npos){ sigma = PythiaParameter(*itr); } else if(itr->find("StringZ:aLund")!=string::npos){ a = PythiaParameter(*itr); } else if(itr->find("StringZ:bLund")!=string::npos){ b = PythiaParameter(*itr); } else if(itr->find("StringFlav:probStoUD")!=string::npos){ rho = PythiaParameter(*itr); } else if(itr->find("StringFlav:probSQtoQQ")!=string::npos){ x = PythiaParameter(*itr); } else if(itr->find("StringFlav:probQQ1toQQ0")!=string::npos){ y = PythiaParameter(*itr); } else if(itr->find("StringFlav:probQQtoQ")!=string::npos){ xi = PythiaParameter(*itr); } else if(itr->find("OverlapStrings:fragMass")!=string::npos){ double m = PythiaParameter(*itr); m2 = m*m; } else if(itr->find("OverlapStrings:baryonSuppression")!=string::npos){ _bsparameter = PythiaParameter(*itr); } else{ _pythiaSettings.push_back(*itr); } } // Sanity check of parameters if (rho < 0 || rho > 1 || x < 0 || x > 1 || y < 0 || y > 1 || xi < 0 || xi > 1 || sigma < 0 || sigma > 1 || a < 0.0 || a > 2.0 || b < 0.2 || b > 2.0){ cout << "Did you set up sensible initial Pythia 8 values? Remember:" << endl; cout << "0 < a < 2; 0.2 < b < 2; 0 < sigma < 1; 0 < xi < 1; 0 < y < 1; 0 < x < 1; 0 < rho < 1" << endl; } } void OverlapPythiaHandler::CreateSinglePythia(double h){ if(!CalculateEffectiveParameters(h)) cout << "Unexpected error setting up parameters for overlap string mode!" << endl; /* cout << "Settings and parameters for overlap string hadronization" << endl; cout << "bs = " << _bsparameter << " m_cut = " << sqrt(m2) << endl; cout << " h = " << h << " rho_eff = " << rho_eff << " xi_eff = " << xi_eff << " a_eff = " << a_eff << " b_eff = " << b_eff << " x_eff = " << x_eff << " y_eff = " << y_eff << endl; */ vector newPythiaSettings = _pythiaSettings; newPythiaSettings.push_back("StringPT:sigma = " + convert(sigma_eff)); newPythiaSettings.push_back("StringZ:aLund = " + convert(a_eff)); newPythiaSettings.push_back("StringZ:bLund = " + convert(b_eff)); newPythiaSettings.push_back("StringFlav:probStoUD = " + convert(rho_eff)); newPythiaSettings.push_back("StringFlav:probSQtoQQ = " + convert(x_eff)); newPythiaSettings.push_back("StringFlav:probQQ1toQQ0 = " + convert(y_eff)); newPythiaSettings.push_back("StringFlav:probQQtoQ = " + convert(xi_eff)); PythiaPtr tmp; _overlapPythias.push_back(tmp); _overlapPythias.back().pPtr = new Pythia8Interface(); _overlapPythias.back().pPtr->init(*_sf,newPythiaSettings); _overlapPythias.back().h = h; ComparePythias comparePythias; sort(_overlapPythias.begin(),_overlapPythias.end(),comparePythias); } Pythia8Interface * OverlapPythiaHandler::GetPythiaPtr(double h, bool recycle){ // Broken calculations could give h = nan. This recursive function would thus // loop infinitely, as nothing == nan. In that case, throw a Veto, and discard the event // after making sufficient noise. if(std::isnan(h)){ cout << "OverlapPythiaHandler::GetPythiaPtr: h is nan. Fix preceeding calculation before running again." << endl; throw Veto(); } // Option to delete all unused Pythia ptrs every ef. 10000 events. // Usefull for Heavy Ion, which contains rare circumstances if(recycle){ vector newVec; for(vector::iterator itr = _overlapPythias.begin(); itr!=_overlapPythias.end();++itr){ if ( itr->h !=h ){ itr->nullify( ); } else{ newVec.push_back(*itr); } } _overlapPythias = newVec; } // Check if we already have for(vector::iterator itr = _overlapPythias.begin(); itr!=_overlapPythias.end();++itr){ if(itr->h==h){ itr->used(); return itr->getpPtr(); } } //cout << "Creating new Pythia with h = " << h << endl; CreateSinglePythia(h); return GetPythiaPtr(h); } vector OverlapPythiaHandler::GetPythiaParameters(double h){ vector ret; ret.clear(); if(!CalculateEffectiveParameters(h)) cout << "Something went wrong calculating effective Pythia parameters!" << endl; ret.push_back(h); ret.push_back(a_eff); ret.push_back(b_eff); ret.push_back(rho_eff); ret.push_back(xi_eff); ret.push_back(x_eff); ret.push_back(y_eff); return ret; } bool OverlapPythiaHandler::CalculateEffectiveParameters(double h){ if (h <= 0) return false; double hinv = 1.0 / h; rho_eff = pow(rho, hinv); x_eff = pow(x, hinv); y_eff = pow(3.0 * y, hinv) / 3.0; sigma_eff = sigma * sqrt(h); double C1 = _bsparameter * xi * (2 + rho) / (3 + 2*x*rho + 9*y + 6*x*rho*y + x*x*rho*rho + 3*y*x*x*rho*rho); xi_eff = (pow(C1, hinv)/_bsparameter) * (3 + 2*x_eff*rho_eff + 9*y_eff + 6*x_eff*rho_eff*y_eff + x_eff*x_eff*rho_eff*rho_eff + 3*y_eff*x_eff*x_eff*rho_eff*rho_eff) / (2 + rho_eff); if (xi_eff > 1.0) xi_eff = 1.0; C1 = b / (2 + rho); b_eff = C1 * (2 + rho_eff); if (b_eff < 0.2) b_eff = 0.2; if (b_eff > 2.0) b_eff = 2.0; double N = IFragmentationF(a, b); double N_eff = IFragmentationF(a, b_eff); int s = sign(N - N_eff); double da = 0.1; a_eff = a - s * da; do { N_eff = IFragmentationF(a_eff, b_eff); if (sign(N - N_eff) != s) { s = sign(N - N_eff); da /= 10.0; } a_eff -= s * da; if (a_eff < 0.0) {a_eff = 0.0; break;} if (a_eff > 2.0) {a_eff = 2.0; break;} } while (da > 0.00001); return true; } double OverlapPythiaHandler::IFragmentationF(double a, double b){ const int N = 100000; const double dz = 1.0 / N; //increment of z variable double Integral = 0.0; for (double z = dz; z < 1; z += dz) { //numerical integral calculation Integral += pow(1 - z, a) * exp(- b * m2 / z) / z; } return Integral * dz; } double OverlapPythiaHandler::PythiaParameter(string par){ size_t found = par.find("="); string number = par.substr(found+2,par.size()); - string::iterator end_pos = remove(par.begin(), par.end(), ' '); - return atof(number.c_str()); + remove(par.begin(), par.end(), ' '); + return atof(number.c_str()); } int OverlapPythiaHandler::sign(double num){ return (num < 0) ? -1 : 1; } string OverlapPythiaHandler::convert(double d) { ostringstream os; os << d; return os.str(); } diff --git a/Hadronization/Ropewalk.cc b/Hadronization/Ropewalk.cc --- a/Hadronization/Ropewalk.cc +++ b/Hadronization/Ropewalk.cc @@ -1,1262 +1,1266 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the Ropewalk class. // #include "Ropewalk.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/EventRecord/Step.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/CurrentGenerator.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/Utilities/Selector.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/Utilities/UtilityBase.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/Handlers/StepHandler.h" using namespace TheP8I; Ropewalk::Ropewalk(const vector & singlets, Length width, Energy scale, double jDP, bool throwaway, bool shoving, double rcIn, StringTension gaIn, double geIn, double ycIn, double dyIn, Length tsIn, Length dtIn, bool lorentzIn, bool cylinderIn, Energy pTcutIn, double mstringyIn, bool longSoftIn, HistKeeper* hk, bool deb) : R0(width), m0(scale), junctionDiquarkProb(jDP), rCutOff(rcIn*width), gAmplitude(gaIn), gExponent(geIn), yCutOff(ycIn), deltay(dyIn), tshove(tsIn), deltat(dtIn), lorentz(lorentzIn), cylinderShove(cylinderIn), stringpTCut(pTcutIn), minStringy(mstringyIn), longSoft(longSoftIn), hkPtr(hk), debug(deb), strl0(0.0), strl(0.0), avh(0.0), avp(0.0), avq(0.0) { for ( int is = 0, Ns = singlets.size(); is < Ns; ++is ){ if(!longSoft) stringdipoles[cloneToFinal(singlets[is])]; else{ tcPVector pVec = singlets[is].partons(); double maxrap = 0; double minrap = 0; Energy pTmax = 0*GeV; for(int id = 0, Nd = int(pVec.size()); id < Nd; ++id){ double y = pVec[id]->rapidity(); Energy pT = pVec[id]->momentum().perp(); if(id == 0){ maxrap = y; minrap = y; pTmax = pT; } else{ if(pTmax < pT) pTmax = pT; if(maxrap < y) maxrap = y; if(minrap > y) minrap = y; } } if(hkPtr){ hkPtr->hist("stringlength")->fill(maxrap-minrap,hkPtr->currentWeight()); hkPtr->hist("stringptMax")->fill(pTmax/GeV,hkPtr->currentWeight()); if(maxrap - minrap > minStringy) hkPtr->hist("stringpTMaxCut")->fill(pTmax/GeV,hkPtr->currentWeight()); } if(maxrap-minrap > minStringy && pTmax < stringpTCut) stringdipoles[cloneToFinal(singlets[is])]; else spectatorDipoles[cloneToFinal(singlets[is])]; } } getDipoles(); calculateOverlaps(false); if(throwaway) doBreakups(); lambdaBefore = lambdaSum(); if ( debug ) CurrentGenerator::log() << singlets.size() << "s " << dipoles.size() << "d "; if(hkPtr){ for(DipoleMap::iterator itr = stringdipoles.begin(); itr!=stringdipoles.end();++itr){ vector dips = itr->second; for(int id = 0, Nd = dips.size(); id < Nd; ++id){ Dipole* d = dips[id]; LorentzMomentum mpa = d->pa->momentum(); LorentzMomentum mpc = d->pc->momentum(); LorentzMomentum mpaRot = (d->R)*mpa; LorentzMomentum mpcRot = (d->R)*mpc; if(d->pa->id() == 21){ hk->hist("normalGluonPt")->fill(mpa.perp()/GeV, hk->currentWeight()); hk->hist("normalGluonRestFramePt")->fill( mpaRot.perp()/GeV,hk->currentWeight()); } if(abs(d->pa->id()) < 10 ){ hk->hist("normalQuarkPt")->fill(mpa.perp()/GeV, hk->currentWeight()); hk->hist("normalQuarkRestFramePt")->fill( mpaRot.perp()/GeV,hk->currentWeight()); } if(abs(d->pc->id()) < 10 ){ hk->hist("normalQuarkPt")->fill(mpc.perp()/GeV, hk->currentWeight()); hk->hist("normalQuarkRestFramePt")->fill( mpcRot.perp()/GeV,hk->currentWeight()); } if(abs(d->pc->id()) > 100 ){ hk->hist("normalDiQuarkPt")->fill(mpc.perp()/GeV, hk->currentWeight()); hk->hist("normalDiQuarkRestFramePt")->fill( mpcRot.perp()/GeV,hk->currentWeight()); } if(abs(d->pa->id()) > 100 ){ hk->hist("normalDiQuarkPt")->fill(mpa.perp()/GeV, hk->currentWeight()); hk->hist("normalDiQuarkRestFramePt")->fill( mpaRot.perp()/GeV,hk->currentWeight()); } } } } if(shoving){ shoveTheDipoles(); // update color singlets DipoleMap newStringDipoles; for(DipoleMap::iterator itr = stringdipoles.begin(); itr != stringdipoles.end(); ++itr){ delete itr->first; // The old dipoles vector oldDipoles = itr->second; tcParticleSet updatedParticles; for(int id = 0, Nd = oldDipoles.size(); id < Nd; ++id){ Dipole* d = oldDipoles[id]; const ParticleVector children = d->pa->children(); if(children.empty()) updatedParticles.insert(d->pa); for(int ic = 0, Nc = children.size(); ic < Nc; ++ic) updatedParticles.insert(children[ic]); if(d->pc->id() != 21){ const ParticleVector children = d->pc->children(); if(children.empty()) updatedParticles.insert(d->pc); for(int ic = 0, Nc = children.size(); ic < Nc; ++ic) updatedParticles.insert(children[ic]); } } tcPPtr p = *updatedParticles.begin(); if(p->colourLine()) newStringDipoles[new ColourSinglet(p->colourLine(), updatedParticles)]; else if(p->antiColourLine()) newStringDipoles[new ColourSinglet(p->antiColourLine(), updatedParticles)]; else Throw() << "Updating colour singlets after shoving failed." << "This is a serious error, please contact the authors." << Exception::abortnow; } stringdipoles = newStringDipoles; getDipoles(); calculateOverlaps(false); if(hkPtr){ for(DipoleMap::iterator itr = stringdipoles.begin(); itr!=stringdipoles.end();++itr){ vector dips = itr->second; for(int id = 0, Nd = dips.size(); id < Nd; ++id){ Dipole* d = dips[id]; LorentzMomentum mpa = d->pa->momentum(); LorentzMomentum mpc = d->pc->momentum(); LorentzMomentum mpaRot = (d->R)*mpa; LorentzMomentum mpcRot = (d->R)*mpc; if(d->pa->id() == 21){ hk->hist("excitedEndGluonPt")->fill(mpa.perp()/GeV, hk->currentWeight()); hk->hist("excitedEndGluonRestFramePt")->fill( mpaRot.perp()/GeV,hk->currentWeight()); } if(abs(d->pa->id()) < 10 ){ hk->hist("excitedQuarkPt")->fill(mpa.perp()/GeV, hk->currentWeight()); hk->hist("excitedQuarkRestFramePt")->fill( mpaRot.perp()/GeV,hk->currentWeight()); } if(abs(d->pc->id()) < 10 ){ hk->hist("excitedQuarkPt")->fill(mpc.perp()/GeV, hk->currentWeight()); hk->hist("excitedQuarkRestFramePt")->fill( mpcRot.perp()/GeV,hk->currentWeight()); } if(abs(d->pc->id()) > 100 ){ hk->hist("excitedDiQuarkPt")->fill(mpc.perp()/GeV, hk->currentWeight()); hk->hist("excitedDiQuarkRestFramePt")->fill( mpcRot.perp()/GeV,hk->currentWeight()); } if(abs(d->pa->id()) > 100 ){ hk->hist("excitedDiQuarkPt")->fill(mpa.perp()/GeV, hk->currentWeight()); hk->hist("excitedDiQuarkRestFramePt")->fill( mpaRot.perp()/GeV,hk->currentWeight()); } } } } } } Ropewalk::~Ropewalk() { for ( DipoleMap::iterator it = stringdipoles.begin(); it != stringdipoles.end(); ++it ) delete it->first; } ColourSinglet * Ropewalk::cloneToFinal(const ColourSinglet & cs) { tcParticleSet final; for ( int i = 0, N = cs.partons().size(); i < N; ++i ) final.insert(cs.partons()[i]->final()); tcPPtr p = *final.begin(); if(hkPtr){ hkPtr->hist("nPartons")->fill(final.size(), hkPtr->currentWeight()); } if ( p->colourLine() ) return new ColourSinglet(p->colourLine(), final); if ( p->antiColourLine() ) return new ColourSinglet(p->antiColourLine(), final); Throw() << "Cloning ColourSinglets failed in Ropewalk. " << "This is a serious error - please contact the authors." << Exception::abortnow; return 0; } LorentzPoint Ropewalk:: propagate(const LorentzPoint & b, const LorentzMomentum & p) const { return b + p/(p.e()*m0/hbarc); // return b; } void Ropewalk::getDipoles() { dipoles.clear(); // First extract all dipoles from all strings (pieces). for ( DipoleMap::iterator it = stringdipoles.begin(); it != stringdipoles.end(); ++it ) { ColourSinglet & string = *it->first; for ( int isp = 1, Nsp = string.nPieces(); isp <= Nsp; ++isp ) { const ColourSinglet::StringPiece & sp = string.piece(isp); if ( sp.size() < 2 ) continue; bool forward = sp[0]->colourLine() && sp[0]->colourLine() == sp[1]->antiColourLine(); for ( int i = 0, N = sp.size(); i < (N - 1); ++i ) { if ( forward ) dipoles.push_back(Dipole(sp[i], sp[i + 1], string, m0)); else dipoles.push_back(Dipole(sp[i + 1], sp[i], string, m0)); } if ( !forward && sp.front()->colourLine() && sp.front()->colourLine() == sp.back()->antiColourLine() ) dipoles.push_back(Dipole(sp.front(), sp.back(), string, m0)); else if ( forward && sp.front()->antiColourLine() && sp.front()->antiColourLine() == sp.back()->colourLine() ) dipoles.push_back(Dipole(sp.back(), sp.front(), string, m0)); } } for ( int i = 0, N = dipoles.size(); i < N; ++i ) { stringdipoles[dipoles[i].string].push_back(&dipoles[i]); strl0 += dipoles[i].yspan(1.0*GeV); } } double Ropewalk::limitrap(const LorentzMomentum & p) const { if ( p.z() == ZERO ) return 0.0; Energy mt = sqrt(max(p.perp2() + p.m2(), sqr(m0))); double rap = log((p.t() + abs(p.z()))/mt); return p.z() > ZERO? rap: -rap; } Ropewalk::OverlappingDipole:: OverlappingDipole(const Dipole & d, const LorentzRotation R, const Ropewalk * rw) : dipole(&d), dir(1) { // Get the boost to other dipole's rest frame and calculate // rapidities. bc = R*rw->propagate(d.pc->vertex(), d.pc->momentum()); ba = R*rw->propagate(d.pa->vertex(), d.pa->momentum()); yc = rw->limitrap(R*d.pc->momentum()); ya = rw->limitrap(R*d.pa->momentum()); if ( yc < ya ) dir = -1; } void Ropewalk::calculateOverlaps(bool doPropagate) { // Then calculate all overlaps between pairs of dipoles. for ( int i1 = 0, Nd = dipoles.size(); i1 < Nd; ++i1 ) { Dipole & d1 = dipoles[i1]; d1.clear(); if ( d1.s() <= sqr(m0) ) continue; // Get the boost to dipoles rest frame and calculate rapidities. LorentzMomentum pc1 = d1.pc->momentum(); LorentzMomentum pa1 = d1.pa->momentum(); //LorentzRotation R = Utilities::boostToCM(make_pair(&pc1, &pa1)); if(doPropagate){ d1.bc = d1.R*propagate(d1.pc->vertex(), d1.pc->momentum()); d1.ba = d1.R*propagate(d1.pa->vertex(), d1.pa->momentum()); } else{ d1.bc = d1.pc->vertex(); d1.ba = d1.pa->vertex(); } double yc1 = limitrap(pc1); double ya1 = limitrap(pa1); double n = 0.0; double m = 0.0; if ( yc1 <= ya1 ) continue; // don't care about too small strings. for ( int i2 = 0; i2 < Nd; ++i2 ) { if ( i1 == i2 ) continue; // Boost to the rest frame of dipole 1. if ( dipoles[i2].s() <= sqr(m0) ) continue; OverlappingDipole od(dipoles[i2], d1.R, this); // Ignore if not overlapping in rapidity. if ( min(od.yc, od.ya) > yc1 || max(od.yc, od.ya) < ya1 || od.yc == od.ya ) continue; // Calculate rapidity overlap. double yomax = min(max(od.yc, od.ya), yc1); double yomin = max(min(od.yc, od.ya), ya1); // Sample a random point in the rapidity overlap region and get // positions in space and check overlap. double yfrac = (yomax - yomin)/(yc1 - ya1); double y = UseRandom::rnd(yomin, yomax); if ( !od.overlap(y, d1.ba + (d1.bc - d1.ba)*(y - ya1)/(yc1 - ya1), R0) ) yfrac = 0.0; // Sum overlaps. if ( od.dir > 0 ) { n += yfrac; } else { m += yfrac; } od.yfrac = yfrac*od.dir; d1.overlaps.push_back(od); } d1.n = int(n + UseRandom::rnd()); d1.m = int(m + UseRandom::rnd()); d1.initWeightMultiplet(); avp += d1.p*d1.yspan(1.0*GeV); avq += d1.q*d1.yspan(1.0*GeV); } } double Ropewalk::Dipole::reinit(double ry, Length R0, Energy m0) { // First calculate the rapidity in this restframe. double y = yspan(m0)*(ry - 0.5); // Then get the position in space of the string at this rapidity. LorentzPoint b = ba + (bc - ba)*ry; p = 1; q = 0; // Now go through all overlapping dipoles. for ( int i = 0, N = overlaps.size(); i < N; ++i ) { if ( overlaps[i].dipole->broken ) continue; if( overlaps[i].dipole->hadr ) continue; if ( overlaps[i].overlap(y, b, R0) ) { if ( overlaps[i].dir > 0 ) ++p; else ++q; } } return kappaEnhancement(); } void Ropewalk::Dipole::initMultiplet() { typedef pair Plet; int ns = 0; int ms = 0; while ( ns < n || ms < m ) { Plet diff; double octetveto = double(ns + ms + 1 - p - q)/double(ns + ms + 1); if ( double(n - ns) > UseRandom::rnd()*double(n + m - ns - ms) ) { Selector sel; sel.insert(multiplicity(p + 1, q), Plet(1, 0)); sel.insert(multiplicity(p, q - 1)*octetveto, Plet(0, -1)); sel.insert(multiplicity(p - 1, q + 1), Plet(-1, 1)); diff = sel[UseRandom::rnd()]; ++ns; } else { Selector sel; sel.insert(multiplicity(p, q + 1), Plet(0, 1)); sel.insert(multiplicity(p - 1, q)*octetveto, Plet(-1, 0)); sel.insert(multiplicity(p + 1, q - 1), Plet(1, -1)); diff = sel[UseRandom::rnd()]; ++ms; } if ( diff.first == -diff.second ) nb++; p += diff.first; q += diff.second; } /* *** ATTENTION *** maybe better if Christian implements this */ } void Ropewalk::Dipole::initNewMultiplet() { typedef pair Plet; int ns = 0; int ms = 0; for ( int i = 0, N = overlaps.size(); i < N; ++i ) { if ( abs(overlaps[i].yfrac) < UseRandom::rnd() ) continue; Plet diff(0,0); double octetveto = double(ns + ms + 1 - p - q)/double(ns + ms + 1); if ( overlaps[i].dir > 0 ) { Selector sel; sel.insert(multiplicity(p + 1, q), Plet(1, 0)); sel.insert(multiplicity(p, q - 1)*octetveto, Plet(0, -1)); sel.insert(multiplicity(p - 1, q + 1), Plet(-1, 1)); diff = sel[UseRandom::rnd()]; ++ns; } else { Selector sel; sel.insert(multiplicity(p, q + 1), Plet(0, 1)); sel.insert(multiplicity(p - 1, q)*octetveto, Plet(-1, 0)); sel.insert(multiplicity(p + 1, q - 1), Plet(1, -1)); diff = sel[UseRandom::rnd()]; ++ms; } if ( diff.first == -diff.second ) nb++; p += diff.first; q += diff.second; } n = ns; m = ms; /* *** ATTENTION *** maybe better if Christian implements this */ } void Ropewalk::Dipole::initWeightMultiplet() { typedef pair Plet; int ns = 0; int ms = 0; double mab = sqrt(s())/GeV; for ( int i = 0, N = overlaps.size(); i < N; ++i ) { if ( abs(overlaps[i].yfrac) < UseRandom::rnd() ) continue; const Dipole * od = overlaps[i].dipole; double mcd = sqrt(od->s())/GeV; double w8 = 1.0/sqr(mab*mcd); //The weight for 3 + 3bar -> 8 double w6 = w8; // The weight for 3 + 3 -> 6 Plet diff(0,0); double octetveto = double(ns + ms + 1 - p - q)/double(ns + ms + 1); if ( overlaps[i].dir > 0 ) { - double mac = (pc->momentum() + od->pc->momentum()).m()/GeV; - double mbd = (pa->momentum() + od->pa->momentum()).m()/GeV; - double w1 = octetveto/sqr(mac*mbd); // The weight for 3 + 3bar -> 1 - double w3 = 1.0/(mab*mcd*mac*mbd); //The weight for 3 + 3 -> 3bar - Selector sel; - sel.insert(multiplicity(p + 1, q) + - multiplicity(p, q - 1)*w8/(w8 + w1) + - multiplicity(p - 1, q + 1)*w6/(w6 + w3), Plet(1, 0)); - sel.insert(multiplicity(p, q - 1)*w1/(w1 + w8), Plet(0, -1)); - sel.insert(multiplicity(p - 1, q + 1)*w3/(w3 + w6), Plet(-1, 1)); + double mac = max((pc->momentum() + od->pc->momentum()).m()/GeV, 1.0e-10); + double mbd = max((pa->momentum() + od->pa->momentum()).m()/GeV, 1.0e-10); + Selector sel; + if ( mac <= ZERO || mbd <= ZERO ) { + sel.insert(multiplicity(p + 1, q), Plet(1, 0)); + } else { + double w1 = octetveto/sqr(mac*mbd); // The weight for 3 + 3bar -> 1 + double w3 = 1.0/(mab*mcd*mac*mbd); //The weight for 3 + 3 -> 3bar + sel.insert(multiplicity(p + 1, q) + + multiplicity(p, q - 1)*w8/(w8 + w1) + + multiplicity(p - 1, q + 1)*w6/(w6 + w3), Plet(1, 0)); + sel.insert(multiplicity(p, q - 1)*w1/(w1 + w8), Plet(0, -1)); + sel.insert(multiplicity(p - 1, q + 1)*w3/(w3 + w6), Plet(-1, 1)); + } diff = sel[UseRandom::rnd()]; ++ns; } else { - double mac = (pa->momentum() + od->pc->momentum()).m()/GeV; - double mbd = (pc->momentum() + od->pa->momentum()).m()/GeV; + double mac = max((pa->momentum() + od->pc->momentum()).m()/GeV, 1.0e-10); + double mbd = max((pc->momentum() + od->pa->momentum()).m()/GeV, 1.0e-10); if ( mac*mbd <= 0.0 ) { ++q; ++ms; continue; } double w1 = octetveto/sqr(mac*mbd); // The weight for 3 + 3bar -> 1 double w3 = 1.0/(mab*mcd*mac*mbd); //The weight for 3 + 3 -> 3bar Selector sel; sel.insert(multiplicity(p, q + 1) + multiplicity(p - 1, q)*w8/(w8 + w1) + multiplicity(p + 1, q - 1)*w6/(w6 + w3), Plet(0, 1)); sel.insert(multiplicity(p - 1, q)*w1/(w1 + w8), Plet(-1, 0)); sel.insert(multiplicity(p + 1, q - 1)*w3/(w3 + w6), Plet(1, -1)); diff = sel[UseRandom::rnd()]; ++ms; } if ( diff.first == -diff.second ) nb++; p += diff.first; q += diff.second; } n = ns; m = ms; /* *** ATTENTION *** maybe better if Christian implements this */ } double Ropewalk::getNb(){ double nba = 0; for ( int i = 0, N = dipoles.size(); i < N; ++i ){ nba += double(dipoles[i].nb); //cout << dipoles[i].n << " " << dipoles[i].m << " " << dipoles[i].p << " " << dipoles[i].q << endl; //cout << double(dipoles[i].n + dipoles[i].m + 1 - dipoles[i].p - dipoles[i].q) << endl; //cout << " --- " << endl; //nba += double(dipoles[i].n + dipoles[i].m + 1 - dipoles[i].p - dipoles[i].q); } // cout << nba << endl; return nba; } double Ropewalk::getkW(){ double kw = 0; for ( int i = 0, N = dipoles.size(); i < N; ++i ) kw += dipoles[i].kappaEnhancement()*log(dipoles[i].s()/sqr(m0)); return kw; } pair Ropewalk::getMN(){ pair ret = make_pair(0,0); for (int i = 0, N = dipoles.size(); i < N; ++i){ ret.first += dipoles[i].m; ret.second += dipoles[i].n; } return ret; } double Ropewalk::lambdaSum(double cutoff){ double ret = 0; for (int i = 0, N = dipoles.size(); i < N; ++i) ret += dipoles[i].s() > cutoff*GeV2 ? log(dipoles[i].s()/sqr(m0)) : 0; return ret; } double Ropewalk::Dipole::breakupProbability() const { if ( !breakable() || n + m <= 0.0 || n + m + 1 == p + q ) return 0.0; double sum = 0.0; for (int j = 0, N = overlaps.size(); j < N; ++j ) if ( overlaps[j].dipole->breakable() ) sum += abs(overlaps[j].yfrac); if ( sum <= 0.0 ) return 1.0; return double(n + m + 1 - p - q)/(sum + 1.0); } Step & Ropewalk::step() const { return *(CurrentGenerator()->currentStepHandler()->newStep()); } Step & Ropewalk::Dipole::step() const { return *(CurrentGenerator()->currentStepHandler()->newStep()); } // Helper class to describe a pair of excitations struct Exc { // The constructor Exc(double yIn, Energy gmassIn, PPtr p1In, PPtr p2In, Ropewalk::Dipole* dip1In, Ropewalk::Dipole* dip2In) : y(yIn), gluonMass(gmassIn), pp1(p1In), pp2(p2In), dip1(dip1In), dip2(dip2In) { // Set a sensible starting vertex pp1->setVertex(dip1->pointAtRap(y)); pp2->setVertex(dip2->pointAtRap(y)); // Set a pointer to the excitation in the dipole dip1->addExcitation(y, pp1); dip2->addExcitation(y, pp2); } // Give the excitation a kick in the x and y direction, void shove(Energy dpx, Energy dpy, HistKeeper* hkPtr){ // The old momenta LorentzMomentum p2 = pp2->momentum(); LorentzMomentum p1 = pp1->momentum(); // The new momenta, after the shove Energy mt2new = sqrt(sqr(p2.x() - dpx) + sqr(p2.y() - dpy) + sqr(gluonMass)); Energy e2new = mt2new*cosh(y); Energy p2znew = mt2new*sinh(y); LorentzMomentum p2new(p2.x() - dpx, p2.y() - dpy, p2znew, e2new); Energy mt1new = sqrt(sqr(p1.x() + dpx) + sqr(p1.y() + dpy) + sqr(gluonMass)); Energy e1new = mt1new*cosh(y); Energy p1znew = mt1new*sinh(y); LorentzMomentum p1new(p1.x() + dpx, p1.y() + dpy, p1znew, e1new); // The differences between the two LorentzMomentum deltap1 = p1new - p1; LorentzMomentum deltap2 = p2new - p2; // We now want to check if we can add these // two excitations to the dipoles if( dip2->recoil(deltap2) ) { if ( dip1->recoil(deltap1) ) { pp1->setMomentum(p1new); pp2->setMomentum(p2new); if(hkPtr){ LorentzMomentum ex2Rot = (dip2->R)*deltap2; LorentzMomentum ex1Rot = (dip1->R)*deltap1; hkPtr->hist("shoveInRestFrame")->fill(ex2Rot.perp()/GeV,hkPtr->currentWeight()); hkPtr->hist("shoveInRestFrame")->fill(ex1Rot.perp()/GeV,hkPtr->currentWeight()); } } else { dip2->recoil(-deltap2); if(hkPtr){ LorentzMomentum ex2Rot = (dip2->R)*deltap2; LorentzMomentum ex1Rot = (dip1->R)*deltap1; hkPtr->hist("rejectedShoveInRestFrame")->fill(ex2Rot.perp()/GeV,hkPtr->currentWeight()); hkPtr->hist("rejectedShoveInRestFrame")->fill(ex1Rot.perp()/GeV,hkPtr->currentWeight()); } } } else if(hkPtr){ LorentzMomentum ex2Rot = (dip2->R)*deltap2; LorentzMomentum ex1Rot = (dip1->R)*deltap1; hkPtr->hist("rejectedShoveInRestFrame")->fill(ex2Rot.perp()/GeV,hkPtr->currentWeight()); hkPtr->hist("rejectedShoveInRestFrame")->fill(ex1Rot.perp()/GeV,hkPtr->currentWeight()); } } LorentzPoint direction(){ // This gives the push direction return (pp2->vertex()-pp1->vertex()); } double y; Energy gluonMass; PPtr pp1, pp2; Ropewalk::Dipole* dip1; Ropewalk::Dipole* dip2; }; void Ropewalk::shoveTheDipoles() { typedef multimap RMap; RMap rapDipoles; // Order the dipoles in maxrapidity and insert new, excited dipole ends double ymin = 0; double ymax = 0; for(int i = 0, N = dipoles.size(); i < N ;++i){ if(dipoles[i].pc->decayed()){ // Set epc to the decay product of pc here dipoles[i].epc = dipoles[i].pc->children()[0]; } else{ // Make a new particle dipoles[i].epc = CurrentGenerator()->getParticle(dipoles[i].pc->id()); dipoles[i].epc->setMomentum(dipoles[i].pc->momentum()); dipoles[i].epc->setVertex(dipoles[i].pc->vertex()); step().addDecayProduct(dipoles[i].pc,dipoles[i].epc, false); } if(dipoles[i].pa->decayed()){ // Set epa to the decay product of pa here dipoles[i].epa = dipoles[i].pa->children()[0]; } else{ dipoles[i].epa = CurrentGenerator()->getParticle(dipoles[i].pa->id()); dipoles[i].epa->setMomentum(dipoles[i].pa->momentum()); dipoles[i].epa->setVertex(dipoles[i].pa->vertex()); step().addDecayProduct(dipoles[i].pa,dipoles[i].epa, false); } // order the dipoles in max rapidity rapDipoles.insert(make_pair(dipoles[i].maxRapidity(), &dipoles[i])); // find maximal and minimal rapidity to sample if(dipoles[i].minRapidity() < ymin) ymin = dipoles[i].minRapidity(); if(dipoles[i].maxRapidity() > ymax) ymax = dipoles[i].maxRapidity(); } vector rapidities; // cout << rCutOff/femtometer << " " << gAmplitude/GeV << " " << gExponent << " " << yCutOff << " " << deltay << " " << tshove/femtometer << " " << deltat/femtometer << endl; // We can let the gluon go off-shell between the pushes, // and on-shell again towards the end // Not fully implemented yet! Do not toggle this on! bool gluonMass = false; // Maybe do better sampling for(double y = ymin; y < ymax; y+=deltay) rapidities.push_back(y); // For each value of ySample, we should have // an excitation pair per dipole we want to treat typedef map > ExMap; ExMap exPairs; Energy gMass = 0.001*keV; for(int i = 0, N = rapidities.size(); i < N; ++i){ // find dipoles that are sampled here // and store them temporarily double ySample = rapidities[i]; vector tmp; for(RMap::iterator rItr = rapDipoles.lower_bound(ySample); rItr != rapDipoles.end(); ++rItr){ if(rItr->second->minRapidity() < ySample) tmp.push_back(rItr->second); } // Construct the excitation particles, one for each sampled dipole vector eParticles; // dipoles to erase vector eraseDipoles; for(int j = 0, M = tmp.size(); j < M; ++j){ // Test if it is possible to add the small excitation to // the dipole Energy pz = gMass*sinh(ySample); Energy e = gMass*cosh(ySample); LorentzMomentum ex = LorentzMomentum(0.0*GeV, 0.0*GeV,pz,e); if(tmp[j]->recoil(ex,true)){ // Add the excitation tmp[j]->recoil(ex); eParticles.push_back(CurrentGenerator()->getParticle(21)); eParticles.back()->setMomentum(ex); } else{ // if it was not possible, remove the dipole from the list // but NOT inside the loop eraseDipoles.push_back(j); } } // Erase dipoles which could not bear an excitation for(int j = 0, M = eraseDipoles.size(); j < M; ++j){ tmp.erase(tmp.begin() + (eraseDipoles[j]-j) ); } // Construct all pairs of possible excitations in this slice exPairs[ySample] = vector(); for(int j = 0, M = tmp.size(); j < M; ++j) for(int k = j + 1; k < M; ++k){ // Both dipoles should span the full rapidity interval //if(tmp[j]->maxRapidity() >= ySample+deltay/2.0 && tmp[k]->maxRapidity() >= ySample+deltay/2.0 // && tmp[j]->minRapidity() <= ySample-deltay/2.0 // && tmp[k]->minRapidity() <= ySample-deltay/2.0){ exPairs[ySample].push_back(Exc(ySample,gMass,eParticles[j],eParticles[k],tmp[j],tmp[k])); //} } } // For all times int ii = 0; for(Length t = 0*femtometer; t < tshove; t+=deltat){ ++ii; // For all slices for(ExMap::iterator slItr = exPairs.begin(); slItr != exPairs.end(); ++slItr) // For all excitation pairs for(int i = 0, N = slItr->second.size(); i < N; ++i){ Exc& ep = slItr->second[i]; LorentzPoint direction = ep.direction(); Length dist = direction.perp(); if(hkPtr && ii < 6){ Histogram* hd = hkPtr->hist("dist"+to_string(ii)); hd->fill(dist/femtometer,hkPtr->currentWeight()); } if(dist < rCutOff){ double gamma = 1.0; if(lorentz){ LorentzMomentum pp = ep.dip1->momentum() + ep.dip2->momentum(); LorentzMomentum p1 = (ep.dip1->R)*pp; LorentzMomentum p2 = (ep.dip2->R)*pp; double gamma1 = sqrt(1 - sqr(p1.perp()/p1.e())); double gamma2 = sqrt(1 - sqr(p2.perp()/p2.e())); gamma = min(gamma1,gamma2); } Energy gain = 0*GeV; if(cylinderShove){ Length R = R0*gExponent; gain = gamma*deltay*deltat*t*gAmplitude*dist/R/R*exp(-dist*dist/2/R/R); /*if(dist > R0*gExponent){ Length R = R0*gExponent; Area area = 2*R*R*acos(dist/2/R) - 0.5*dist*sqrt(4*R*R - dist*dist); Energy inner = gAmplitude*area/M_PI/R/R; gain = gamma*(sqrt((1*GeV + inner)*(1*GeV + inner)) - 1*GeV); } else{ continue; } */ } else{ gain = gamma*deltay*(t+deltat)*deltat*gAmplitude/R0/gExponent/sqrt(2*M_PI)* exp(-dist*dist/2/gExponent/gExponent/R0/R0); } Energy dpx = dist > ZERO ? gain*direction.x()/dist: 0.0*GeV; Energy dpy = dist > ZERO ? gain*direction.y()/dist: 0.0*GeV; // Print stuff if(hkPtr && ii < 6){ Histogram* h = hkPtr->hist("shove"+to_string(ii)); h->fill(gain/MeV,hkPtr->currentWeight()); } ep.shove(dpx,dpy,hkPtr); } } // Propagate the dipoles for(int i = 0, Nd = dipoles.size(); i < Nd; ++i) dipoles[i].propagate(deltat, deltay); } // Now add the excitations to the dipoles for(int i = 0, Nd = dipoles.size(); i < Nd; ++i){ if(dipoles[i].excitations.size() > 0) dipoles[i].excitationsToString(gluonMass, yCutOff, hkPtr); else{ dipoles[i].epc->antiColourConnect(dipoles[i].epa); } } } multimap< double, pair > Ropewalk::Dipole::getDistPairs(){ // Fill a container with the particle pair rapidities in lab system, // ordered by their spans multimap< double, pair > distPairs; DipEx::iterator iplusone = excitations.begin(); ++iplusone; for(DipEx::iterator itr = excitations.begin(); itr != excitations.end(); ++itr,++iplusone){ if(itr == excitations.begin()){ // First -- special LorentzMomentum p1 = (pa->momentum().rapidity() == minRapidity() ? pa->momentum() : pc->momentum()); Energy2 s = (p1 + itr->second->momentum()).m2(); double ysep = (s > sqr(m0) ? 0.5*log(s/sqr(m0)) : 0.0); distPairs.insert(make_pair(ysep, make_pair(minRapidity(),itr->first)) ); } if(iplusone == excitations.end()){ // Last -- special LorentzMomentum p1 = (pa->momentum().rapidity() == maxRapidity() ? pa->momentum() : pc->momentum()); Energy2 s = (p1 + itr->second->momentum()).m2(); double ysep = (s > sqr(m0) ? 0.5*log(s/sqr(m0)) : 0.0); distPairs.insert(make_pair(ysep, make_pair(itr->first,maxRapidity())) ); } else{ Energy2 s = (itr->second->momentum() + iplusone->second->momentum()).m2(); double ysep = (s > sqr(m0) ? 0.5*log(s/sqr(m0)) : 0.0); distPairs.insert(make_pair(ysep, make_pair(itr->first,iplusone->first))); } } return distPairs; } void Ropewalk::Dipole::excitationsToString(bool gluonMass, double yCutOff, HistKeeper* hkPtr){ if(gluonMass){ cout << "Oops! Not implemented yet!" << endl; return; } // Clustering should always be turned on // only turn it off for debugging purposes. // Change the parameters DeltaY or YCutOff instead bool clustering = true; if(clustering){ // Clustering // // For all excitation neighbors in the lab-system, // we calculate their effective rapidity span // and keep only those above yCutOff multimap > distPairs = getDistPairs(); if(hkPtr){ for(auto itr = distPairs.begin(); itr != distPairs.end(); ++itr) hkPtr->hist("yDistBefore")->fill(itr->first,hkPtr->currentWeight()); } if(distPairs.begin() != distPairs.end()) while(distPairs.begin()->first < yCutOff){ // There is only one excitation on the string // The ends will share the momentum provided // that both distances are small multimap >::iterator itr = distPairs.begin(); if(distPairs.size() == 2 && distPairs.rbegin()->first < yCutOff){ DipEx::iterator exP = excitations.find(itr->second.second); if(exP == excitations.end()) exP = excitations.find(itr->second.first); LorentzMomentum panew = epa->momentum() + 0.5*exP->second->momentum(); LorentzMomentum pcnew = epc->momentum() + 0.5*exP->second->momentum(); epa->setMomentum(panew); epc->setMomentum(pcnew); excitations.erase(exP); } // If it is the first piece, put the momentum on the dipole end else if(itr->second.first == minRapidity()){ LorentzMomentum pnew = (epa->momentum().rapidity() == minRapidity() ? epa->momentum() : epc->momentum()); pnew+=excitations[itr->second.second]->momentum(); if(epa->momentum().rapidity() == minRapidity()) epa->setMomentum(pnew); else epc->setMomentum(pnew); // Erase the excitation excitations.erase(itr->second.second); } // If it is the last piece ditto... else if(itr->second.second == maxRapidity()){ LorentzMomentum pnew = (epa->momentum().rapidity() == maxRapidity() ? epa->momentum() : epc->momentum()); pnew+=excitations[itr->second.first]->momentum(); if(epa->momentum().rapidity() == maxRapidity()) epa->setMomentum(pnew); else epc->setMomentum(pnew); // Erase the excitation excitations.erase(itr->second.first); } else{ LorentzMomentum pnew = excitations[itr->second.first]->momentum() + excitations[itr->second.second]->momentum(); // update the first excitations[itr->second.first]->setMomentum(pnew); // erase the second excitations.erase(itr->second.second); } // We can surely do better than recalculating everything // each time - lets see if it kills us if(excitations.empty()) break; distPairs = getDistPairs(); } if(hkPtr){ hkPtr->hist("nExcitations")->fill(excitations.size(),hkPtr->currentWeight()); for(auto itr = distPairs.begin(); itr != distPairs.end(); ++itr) hkPtr->hist("yDistAfter")->fill(itr->first,hkPtr->currentWeight()); } // } // Go through the remaining excitations and sanity check them for(DipEx::iterator itr = excitations.begin(); itr != excitations.end(); ++itr){ // In the dipole rest frame, we should have pz == 0 //LorentzMomentum ppmom = R*itr->second->momentum(); //LorentzMomentum ppn(ppmom.x(),ppmom.y(),0*GeV,ppmom.t()); //itr->second->setMomentum(R.inverse()*ppn); LorentzMomentum ppm = itr->second->momentum(); if(fabs(ppm.x()/GeV) < 1e-5 && fabs(ppm.y()/GeV) < 1e-5 && fabs(ppm.z()/GeV) < 1e-5 ) excitations.erase(itr); else if(hkPtr){ hkPtr->hist("excitationPt")->fill(ppm.perp()/GeV,hkPtr->currentWeight()); LorentzMomentum ppr = R*ppm; hkPtr->hist("excitationRestFramePt")->fill(ppr.perp()/GeV,hkPtr->currentWeight()); } } // If we dont have any excitations, we are done if(excitations.empty()){ epa->colourConnect(epc); return; } // Flag if we start from the colored parton bool color = false; // Link the first excitation to the dipole end with the smallest rapidity if(minRapidity() == epc->momentum().rapidity()){ epa->colourConnect(excitations.rbegin()->second); epc->antiColourConnect(excitations.begin()->second); //excitations.begin()->second->colourConnect(epc); //excitations.rbegin()->second->antiColourConnect(epa); step().addDecayProduct(pc,excitations.begin()->second,false); color = true; } else{ epa->colourConnect(excitations.begin()->second); epc->antiColourConnect(excitations.rbegin()->second); //excitations.begin()->second->antiColourConnect(epa); //excitations.rbegin()->second->colourConnect(epc); step().addDecayProduct(pa,excitations.begin()->second,false); } // If we have only one excitation, we are quickly done if(excitations.size() == 1) return; // We have more... DipEx::iterator iplusone = excitations.begin(); ++iplusone; for(DipEx::iterator itr = excitations.begin(); iplusone != excitations.end(); ++itr,++iplusone){ // Are we halfway through? bool halfway = double(distance(excitations.begin(),iplusone))/double(excitations.size()) > 0.5; if(color){ itr->second->antiColourConnect(iplusone->second); //iplusone->second->colourConnect(itr->second); step().addDecayProduct((halfway ? pa : pc), iplusone->second,false); } else{ itr->second->colourConnect(iplusone->second); //iplusone->second->antiColourConnect(itr->second); step().addDecayProduct((halfway ? pc : pa), iplusone->second,false); } } } void Ropewalk::Dipole::propagate(Length deltat, double deltay){ const Lorentz5Momentum& pcm = epc->momentum(); const Lorentz5Momentum& pam = epa->momentum(); if(pcm.e()/GeV == 0 || pam.e()/GeV == 0) Throw() << "Tried to propagate a dipole end with energy 0; this should never happen." << "This is a serious error - please contact the authors." << Exception::abortnow; // Propagate in dipole rest frame bc = R*(bc + deltat * pcm/(pcm.e())); ba = R*(ba + deltat * pam/(pam.e())); // Propagate in the lab frame LorentzPoint bclab = epc->vertex() + deltat*pcm/pcm.e(); LorentzPoint balab = epa->vertex() + deltat*pam/pam.e(); epc->setVertex(bclab); epa->setVertex(balab); //const_cast((*pc)).setVertex(bclab); //const_cast((*pa)).setVertex(balab); for(DipEx::iterator eItr = excitations.begin(); eItr != excitations.end(); ++eItr){ const Lorentz5Momentum& em = eItr->second->momentum(); // Propagate excitations only in lab frame // only if it has been pushed if(em.e()/GeV > 0){ eItr->second->setVertex(eItr->second->vertex() + deltat*em/(em.perp() + m0*deltay*deltay)); } else{ eItr->second->setVertex(pointAtRap(eItr->first)); } } } void Ropewalk::doBreakups() { typedef multimap BMap; BMap breakups; // First order alldipoles in increasing probability of breaking for ( int i = 0, N = dipoles.size(); i < N; ++i ) breakups.insert(make_pair(dipoles[i].breakupProbability(), &dipoles[i])); // Then start breaking in reverse order mspec.clear(); for ( BMap::reverse_iterator it = breakups.rbegin(); it != breakups.rend(); ++it ) { const Dipole & d = *it->second; if ( d.breakupProbability() < UseRandom::rnd() ) continue; mspec.push_back(d.s()/GeV2); double hinv = 1.0 / d.kappaEnhancement(); // The default parameters should be set // properly from the outside. Now just Pythia8 defaults double rho = pow(0.215, hinv); double x = pow(1.0,hinv); double y = pow(0.027,hinv); Selector qsel; qsel.clear(); if(UseRandom::rnd() < junctionDiquarkProb){ // Take just ordinary quarks qsel.insert(1.0, ParticleID::ubar); qsel.insert(1.0, ParticleID::dbar); qsel.insert(rho, ParticleID::sbar); } else { // Take diquarks // TODO (Pythia 8.2) Set status code of these diquarks to 74 qsel.insert(1.0, ParticleID::ud_0); qsel.insert(3.0*y, ParticleID::dd_1); qsel.insert(3.0*y, ParticleID::ud_1); qsel.insert(3.0*y, ParticleID::uu_1); qsel.insert(rho*x, ParticleID::su_0); qsel.insert(3*x*rho*y, ParticleID::su_1); qsel.insert(3*y*x*x*rho*rho, ParticleID::ss_1); qsel.insert(rho*x, ParticleID::sd_0); qsel.insert(3*x*rho*y, ParticleID::sd_1); } int flav = qsel[UseRandom::rnd()]; PPtr dic = CurrentGenerator()->getParticle(flav); PPtr dia = CurrentGenerator()->getParticle(-flav); //cout << (dic->momentum() + dia->momentum()).m2() / GeV2 << " " << d.s() / GeV2 << endl; if((dic->momentum() + dia->momentum()).m2() > d.s()) continue; // Get boost to dipole rest frame and place the di-quarks there. LorentzRotation R = Utilities::getBoostFromCM(make_pair(d.pc, d.pa)); SimplePhaseSpace::CMS(dic, dia, d.s(), 1.0, 0.0); dia->transform(R); dic->transform(R); // Add the di-quarks as children to the decayed gluons and split // the resulting string. if ( !( step().addDecayProduct(d.pc, dic, true) && step().addDecayProduct(d.pa, dia, true) && step().addDecayProduct(d.pc, dia, false) && step().addDecayProduct(d.pa, dic, false)) ) Throw() << "Adding decay products failed in Ropewalk. " << "This is a serious error - please contact the authors." << Exception::abortnow; tcParticleSet pset(d.string->partons().begin(), d.string->partons().end()); pset.erase(d.pc); pset.erase(d.pa); pset.insert(dic); pset.insert(dia); // Reset the neighboring dipoles to point to the new diquarks. d.broken = true; ColourSinglet * olds = d.string; const vector & oldd = stringdipoles[olds]; for ( int id = 0, Nd = oldd.size(); id < Nd; ++id ) { if ( oldd[id]->pc == d.pa ) oldd[id]->pc = dia; if ( oldd[id]->pa == d.pc ) oldd[id]->pa = dic; } // remove the old string and insert the new ones fixing the // pointers from the dipoles. vector newstrings = ColourSinglet::getSinglets(pset.begin(), pset.end()); for ( int is = 0, Ns = newstrings.size(); is < Ns; ++is ) { ColourSinglet * news = new ColourSinglet(newstrings[is]); tcParticleSet pset(news->partons().begin(), news->partons().end()); vector & newd = stringdipoles[news]; for ( int id = 0, Nd = oldd.size(); id < Nd; ++id ) { if ( !oldd[id]->broken && pset.find(oldd[id]->pc) != pset.end() ) { newd.push_back(oldd[id]); oldd[id]->string = news; } } sort(newd, news->partons()[0]); } stringdipoles.erase(olds); delete olds; } } vector< pair > Ropewalk::getSinglets(double DeltaY) const { vector< pair > ret; ret.reserve(stringdipoles.size()); double toty = 0.0; double totyh = 0.0; double toth = 0.0; double toth2 = 0.0; double minh = 10.0; double maxh = 0.0; // Calculate the kappa-enhancement of the remaining strings. for ( DipoleMap::const_iterator it = stringdipoles.begin(); it != stringdipoles.end(); ++it ) { double sumyh = 0.0; double sumy = 0.0; for ( int id = 0, Nd = it->second.size(); id < Nd; ++id ) { double y = it->second[id]->yspan(m0); double h = it->second[id]->kappaEnhancement(); avh += h*it->second[id]->yspan(1.0*GeV); strl += it->second[id]->yspan(1.0*GeV); if ( DeltaY > 0.0 && abs(it->second[id]->pc->eta()) > DeltaY && abs(it->second[id]->pa->eta()) > DeltaY ) continue; sumy += y; sumyh += y*h; } toty += sumy; totyh += sumyh; if ( sumy > 0.0 ) { double h = 0.1*int(10.0*sumyh/sumy + UseRandom::rnd()); ret.push_back(make_pair(*it->first, h)); toth += sumyh/sumy; toth2 += sqr(sumyh/sumy); minh = min(minh, sumyh/sumy); maxh = max(maxh, sumyh/sumy); } else { ret.push_back(make_pair(*it->first, 1.0)); toth += 1.0; toth2 += 1.0; minh = min(minh, 1.0); maxh = max(maxh, 1.0); } } if ( debug ) CurrentGenerator::log() << ret.size() << "ns " << totyh/max(toty, 0.00001) << "hwa " << minh << "h " << sqrt(toth2/ret.size() - sqr(toth/ret.size())) << "hsd " << endl; if ( avh > 0.0 ) avh /= strl; if ( avp > 0.0 ) avp /= strl0; if ( avq > 0.0 ) avq /= strl0; return ret; } double Ropewalk::averageKappaEnhancement(const vector & dipoles, double DeltaY) const { double sumyh = 0.0; double sumy = 0.0; for ( int id = 0, Nd = dipoles.size(); id < Nd; ++id ) { dipoles[id]->reinit(UseRandom::rnd(), R0, m0); double y = dipoles[id]->yspan(m0); double h = dipoles[id]->firstKappa(); if ( DeltaY > 0.0 && abs(dipoles[id]->pc->eta()) > DeltaY && abs(dipoles[id]->pa->eta()) > DeltaY ) continue; sumy += y; sumyh += y*h; } if ( sumy <= 0.0 ) return 1.0; return sumyh/sumy; } void Ropewalk::sort(vector & dips, tcPPtr p) { if ( p->antiColourLine() ) { map dmap; Dipole * current = 0; for ( int i = 0, N = dips.size(); i < N; ++i ) { if ( dips[i]->pa == p ) current = dips[i]; else dmap[dips[i]->pa] = dips[i]; } dips.clear(); while ( current ) { dips.push_back(current); map::iterator it = dmap.find(current->pc); if ( it == dmap.end() ) current = 0; else { current = it->second; dmap.erase(it); } } if ( !dmap.empty() ) Throw() << "Failed to sort dipoles in Ropewalk. " << "This is a serious error - please contact the authors." << Exception::abortnow; } else { map dmap; Dipole * current = 0; for ( int i = 0, N = dips.size(); i < N; ++i ) { if ( dips[i]->pc == p ) current = dips[i]; else dmap[dips[i]->pc] = dips[i]; } dips.clear(); while ( current ) { dips.push_back(current); map::iterator it = dmap.find(current->pa); if ( it == dmap.end() ) current = 0; else { current = it->second; dmap.erase(it); } } if ( !dmap.empty() ) Throw() << "Failed to sort dipoles in Ropewalk. " << "This is a serious error - please contact the authors." << Exception::abortnow; } } diff --git a/Hadronization/Ropewalk.h b/Hadronization/Ropewalk.h --- a/Hadronization/Ropewalk.h +++ b/Hadronization/Ropewalk.h @@ -1,761 +1,761 @@ // -*- C++ -*- #ifndef TheP8I_Ropewalk_H #define TheP8I_Ropewalk_H // // This is the declaration of the Ropewalk class. // #include "TheP8I/Config/TheP8I.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/EventRecord/ColourSinglet.h" #include "ThePEG/Utilities/UtilityBase.h" #include "ThePEG/Utilities/Throw.h" #include "HistKeeper.h" namespace std { template <> struct less : public binary_function { /** * This is the function called when comparing two pointers to * type_info. */ bool operator()(const ThePEG::ColourSinglet * x, const ThePEG::ColourSinglet * y) const { if ( x && y ) return x->partons() < y->partons(); return !y; } }; } namespace TheP8I { using namespace ThePEG; -typedef Qty<-1,1,0> StringTension; +typedef Tension StringTension; /** * Here is the documentation of the Ropewalk class. */ class Ropewalk { public: double lambdaBefore; vector mspec; /** * Forward declaration of a Dipole. */ struct Dipole; /** * Container class storing information of a dipole overlapping with another; */ struct OverlappingDipole { /** * Standard constructor. */ OverlappingDipole(const Dipole & d, const LorentzRotation R, const Ropewalk * rw); /** * Check if ovelapping at specific rapidity. */ bool overlap(double y, LorentzPoint b1, Length R0) { if ( y < min(ya, yc) || y > max(ya, yc) ) return false; LorentzPoint b2 = ba + (bc - ba)*(y - ya)/(yc - ya); return (b1 - b2).perp() <= R0; } LorentzPoint pointAtRap(double y){ return ba + (bc - ba)*(y - ya)/(yc - ya); } /** * Pointer to the dipole. */ const Dipole * dipole; /** * The boosted rapidities. */ double yc, ya; /** * The propagated and boosted positions. */ LorentzPoint bc, ba; /** * Relative direction. */ int dir; /** * Initially estimated fraction of overlap. */ double yfrac; }; /** * Container class containing information about overlaps for a dipole. * *** ATTENTION *** the meaning of pa and pc is reversed: pa is * always the COLOURED and pc is always the ANTI-COLOURED */ struct Dipole { struct DipoleException : public Exception {}; /** * Constructor. */ Dipole(tcPPtr p1, tcPPtr p2, ColourSinglet & str, Energy m0In) : pc(p2), pa(p1), string(&str), n(0), m(0), p(1), q(0), nb(0), broken(false), hadr(false), m0(m0In) { if(pc->colourLine() != pa->antiColourLine()) swap(pc,pa); LorentzMomentum pcm = pc->momentum(); LorentzMomentum pam = pa->momentum(); //if(pcm.x()/GeV == 0 && pcm.y()/GeV == 0 && pcm.z()/GeV == 0 && pcm.e()/GeV == 0) //cout << pc->id() << " " << pa->id() << endl; //cout << pam.x()/GeV << " " << pam.y()/GeV << " " << pam.z()/GeV << " " << pam.e()/GeV << endl; if ( (pcm + pam).m2() < 0.000001*GeV2 ) Throw() << "Ropewalk found a minisculedipole. This should not happen, " << "but hoping it is just a fluke and throwing away the event." << Exception::eventerror; R = Utilities::boostToCM(make_pair(&pcm,&pam) ); } /** * Return true if this dipole is breakable. Only dipoles between * gluons (not belonging to other dipoles which have broken) and * which has a mass above 2 GeV can break. */ bool breakable() const { return !broken && pc->id() == ParticleID::g && pa->id() == ParticleID::g && !pc->decayed() && !pa->decayed() && s() > 4.0*GeV2; } /** * Calculate the probability that this dipole should break. Taking * into account the number of negative steps in the colour space * and the fraction of overlapping dipoles which are able to * break. */ double breakupProbability() const; void clear(){ overlaps.clear(); n = 0, m = 0, p = 1, q = 0, nb = 0, broken = false, hadr = false; } // Propagate the particles and the excitations void propagate(Length deltat, double deltay); /** * The minimal rapidity in the lab system of the excited dipole ends */ double minRapidity() const { //double y1 = 0.5*log(epc->Pplus()/epc->Pminus()); //double y2 = 0.5*log(epa->Pplus()/epa->Pminus()); double y1 = epc->momentum().rapidity(); double y2 = epa->momentum().rapidity(); return min(y1,y2); } /** * The maximal rapidity in the lab system */ double maxRapidity() const { //double y1 = 0.5*log(epc->Pplus()/epc->Pminus()); //double y2 = 0.5*log(epa->Pplus()/epa->Pminus()); double y1 = epc->momentum().rapidity(); double y2 = epa->momentum().rapidity(); return max(y1,y2); } /* * The b-coordinate (in lab frame) of a point * between the propagated particles, at lab-rapidity yIn */ LorentzPoint pointAtRap(double yIn) const { double ymin = minRapidity(); double ymax = maxRapidity(); LorentzPoint bMin, bMax; if (pc->momentum().rapidity() == ymin ) bMin = pc->vertex(), bMax = pa->vertex(); else bMin = pa->vertex(), bMax = pc->vertex(); // We assume that the rapidity fraction stays the same, // as we cannot calculate yIn in dipole rest frame return bMin + (bMax - bMin)*(yIn - ymin)/(ymax - ymin); } /** * Recalculate overlaps assuming a position a fraction ry from * the coloured parton. */ double reinit(double ry, Length R0, Energy m0); /** * Set and return the effective multiplet. */ void initMultiplet(); void initNewMultiplet(); void initWeightMultiplet(); /** * Calculate the multiplicity given pp and qq; */ static double multiplicity(int pp, int qq) { return ( pp < 0 || qq < 0 || pp + qq == 0 )? 0.0: 0.5*(pp + 1)*(qq + 1)*(pp + qq + 2); } /** * Calculate the average kappa-enhancement. */ double kappaEnhancement() const { // return 1.0; return 0.25*(p + q + 3 - p*q/double(p + q)); } /** * Calculate the kappa-enhancement in the first break-up. */ double firstKappa() const { return 0.25*(2 + 2*p + q); } /** * Calculate boosted kappa-enhancement. */ double firstKappa(double alpha) const { if (alpha > 0) return alpha*firstKappa(); return firstKappa(); } /* LorentzMomentum shoveOverlappers(double dipFrac, double m2Had, Length R0, Lorentz5Momentum pH){ // We first find the b-space point at the breaking point LorentzPoint thisPoint = ba + (bc - ba)*dipFrac; double thisy = pa->momentum().rapidity() + (pc->momentum().rapidity() - pa->momentum().rapidity())*dipFrac; Energy mp = sqrt(pH.mass2()); Energy mtd = sqrt(s()); Energy mSumH = sqrt(m2Had)*GeV; double deltay = mp/(mtd - (mSumH - dipFrac*mtd)/(1 - dipFrac) ); // We then go through all dipoles // overlapping in rapidity Energy ex = 0*GeV, ey = 0*GeV; for(int i = 0, N = overlaps.size(); i < N; ++i){ OverlappingDipole& od = overlaps[i]; // Must overlap at the breaking point if (thisy < min(od.ya, od.yc) || thisy > max(od.ya, od.yc)) continue; LorentzPoint otherPoint = od.pointAtRap(thisy); LorentzPoint direction = thisPoint - otherPoint; // The distance between the two dipoles Length dist = direction.perp(); Energy gain = 2*fabs(deltay)*exp(-dist*dist/4/R0/R0)*GeV; ex+=gain*direction.x()/dist; ey+=gain*direction.y()/dist; od.dipole->recoil(dipFrac,gain,direction); } Lorentz5Momentum ret(ex, ey, 0*GeV, 0*GeV); return ret+pH; } */ /** * Push a dipole */ /*void recoil(double dipFrac, Energy gain, LorentzPoint direction) const { LorentzMomentum pPush(-gain*direction.x()/direction.perp(), -gain*direction.y()/direction.perp(), 0*GeV, 0*GeV); (void) pPush; // We go to the dipole rest frame // and calculate the push on the ends //Lorentz5Momentum pcPush = R*pPush + pc->momentum(); //Lorentz5Momentum paPush = R*pPush + pa->momentum(); // const_cast((*pc)).setMomentum(pcPush); // const_cast((*pa)).setMomentum(paPush); } */ /** * return the momentum of this dipole */ LorentzMomentum momentum() const { return pc->momentum() + pa->momentum(); } /** * Return the squared invariant mass of this dipole. */ Energy2 s() const { return (pc->momentum() + pa->momentum()).m2(); } /** * Return the effective rapidityspan of this dipole. */ double yspan(Energy m0) const { return s() > sqr(m0)? 0.5*log(s()/sqr(m0)): 0.0; } /** * The coloured particle. */ mutable tcPPtr pc; /** * The anti-coloured particle. */ mutable tcPPtr pa; /** * The coloured, excited particle */ PPtr epc; /** * The anti-coloured, excited particle */ PPtr epa; /* * All excitations belonging to this dipole ordered in * rapidity in lab system */ typedef map DipEx; DipEx excitations; void printExcitations() { for(DipEx::iterator itr = excitations.begin(); itr != excitations.end(); ++itr){ cout << itr->first << " " << itr->second->rapidity() << " " << itr->second->momentum().perp()/GeV << endl; } } /* * Recoil the dipole from adding a gluon * this function will not actually add the gluon, only the * recoil, if possible, and return true or false if it was */ bool recoil(const LorentzMomentum& pg, bool dummy = false){ // Keep track of direction int sign = 1; if(epc->rapidity() > epa->rapidity()) sign = -1; // lc momenta after inserting the gluon Energy pplus = epc->Pplus() + epa->Pplus() - pg.plus(); Energy pminus = epc->Pminus() + epa->Pminus() - pg.minus(); // The new lightcone momenta of the dipole ends Energy ppa = 0.0*GeV; Energy ppc = 0.0*GeV; Energy pma = 0.0*GeV; Energy pmc = 0.0*GeV; Energy2 mta2 = epa->mt2(); Energy2 mtc2 = epc->mt2(); Energy mta = sqrt(mta2); Energy mtc = sqrt(mtc2); if ( pplus*pminus <= sqr(mta + mtc) || pplus <= ZERO || pminus <= ZERO ) return false; Energy4 sqarg = sqr(pplus*pminus + epa->mt2() - epc->mt2()) - 4.0*epa->mt2()*pplus*pminus; // Calculate the new momenta if(sqarg <= ZERO ) return false; if(sign > 0){ ppa = (pplus*pminus + epa->mt2() - epc->mt2() + sqrt(sqarg))*0.5/pminus; pma = epa->mt2()/ppa; pmc = pminus - pma; ppc = epc->mt2()/pmc; if ( ppa*mtc < ppc*mta ) return false; // Check rapidity after. } else{ pma = (pplus*pminus + epa->mt2() - epc->mt2() + sqrt(sqarg))*0.5/pplus; ppa = epa->mt2()/pma; ppc = pplus - ppa; pmc = epc->mt2()/ppc; if ( ppa*mtc > ppc*mta ) return false; // Check rapidity after. } LorentzMomentum shifta = LorentzMomentum(epa->momentum().x(),epa->momentum().y(), 0.5*(ppa - pma),0.5*(ppa + pma)); LorentzMomentum shiftc = LorentzMomentum(epc->momentum().x(), epc->momentum().y(), 0.5*(ppc - pmc), 0.5*(ppc + pmc)); // Assign the new momenta if(!dummy){ epa->setMomentum(shifta); epc->setMomentum(shiftc); } return true; } /* * Insert an excitation provided that it is not already there */ void addExcitation(double ylab, PPtr ex){ auto ret = excitations.equal_range(ylab); for(auto itr = ret.first; itr != ret.second; ++itr ) if(ex == itr->second){ return; } excitations.insert(make_pair(ylab,ex)); } /* * Add all excitations to the dipole */ void excitationsToString(bool gluonMass, double yCutOff, HistKeeper* hkPtr = NULL); /* * Extract the distance between all adjecant excitation pairs on a dipole */ multimap< double, pair > getDistPairs(); /** * Return the current step. */ Step & step() const; /** * The propagated and boosted positions. */ LorentzPoint bc, ba; /** * The Lorentz rotation taking us to this dipole's rest system */ LorentzRotation R; /** * The overlapping dipoles with the amount of overlap (negative * means anti overlap). */ vector overlaps; /** * The string to which this Dipole belongs. */ ColourSinglet * string; /** * The summed parallel (n) and anti-parallel (m) overlaps from * other dipoles. */ int n, m; /** * The multiplet. */ int p, q; /** * The number of baryons from junctions */ int nb; /** * Indicate that this dipole has been broken. */ mutable bool broken; mutable bool hadr; /* * The m0 parameter */ Energy m0; }; /** * Convenient typedef */ typedef map > DipoleMap; public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ Ropewalk(const vector & singlets, Length width, Energy scale, double jDC, bool throwaway = true, bool shoving = true, double rcIn = 5.0, StringTension gaIn = 1.0*GeV/femtometer, double geIn = 1.0, double ycIn = 4.0, double dyIn = 0.1, Length tsIn = 1.0*femtometer, Length dtIn = 1.0*femtometer, bool lorentzIn = true, bool cylinderIn = false, Energy ptCutIn = 3.0*GeV, double mstringyIn = 8.0, bool longSoftIn = false, HistKeeper* hk = NULL, bool deb = false); /** * The destructor. */ virtual ~Ropewalk(); //@} /** * Return all ColourSinglet objects with their kappa enhancement. */ vector< pair > getSinglets(double DeltaY = 0.0) const; /** * Calculate the average kappa-enhancement for the given colour * singlet, optionally only taking into account dipoles the * rapidity interval |y| < deltaY. Note that \a cs must point to a * object in the stringdipoles map. @return -1 if cs is not found. */ double averageKappaEnhancement(ColourSinglet * cs, double deltaY = 0.0) const { DipoleMap::const_iterator it = stringdipoles.find(cs); if ( it == stringdipoles.end() ) return -1.0; return averageKappaEnhancement(it, deltaY); } /** * Calculate the average kappa-enhancement for the given iterator * pointing into the stringdipoles map, optionally only taking * into account dipoles the rapidity interval |y| < deltaY. Note * that \a cs must point to a object in the stringdipoles * map. @return -1 if something went wrong. */ double averageKappaEnhancement(DipoleMap::const_iterator it, double deltaY = 0.0) const { return averageKappaEnhancement(it->second, deltaY); } /** * Calculate the average kappa-enhancement the given vector of * dipoles, optionally only taking into account dipoles the * rapidity interval |y| < deltaY. @return -1 if something went * wrong. */ double averageKappaEnhancement(const vector & dipoles, double deltaY = 0.0) const; /** * Propagate a parton a finite time inthe lab-system. */ LorentzPoint propagate(const LorentzPoint & b, const LorentzMomentum & p) const; /** / Get number of juctions */ double getNb(); double getkW(); /** / Get m,n for all strings */ pair getMN(); /** / Get lambda measure for event */ double lambdaSum(double cutoff = 0.1); /** * Calculate the rapidity of a Lorentz vector but limit the * transverse mass to be above m0. */ double limitrap(const LorentzMomentum & p) const; /** * Sort the given vector of dipoles so that the particle \a p is * first and the others follows in sequence. If \a p is a * (anti-)triplet, it will be the pc(pa) of the first dipole, and * the pc(pa) of the second will be the pa(pc) of the first and so * on. If \a p is an octet, the direction will be as if it was a * triplet. */ static void sort(vector & dips, tcPPtr p); protected: /** * Extract all dipoles. */ void getDipoles(); /** * Calculate all overlaps and initialize multiplets in all dipoles. */ void calculateOverlaps(bool doPropagate); /** * Break dipoles (and strings) */ void doBreakups(); /** * Shove the dipoles away from each other */ void shoveTheDipoles(); /** * Return the current step. */ Step & step() const; /** * Make a copy of a ColourSinglet making sure all partons are final. */ ColourSinglet * cloneToFinal(const ColourSinglet & cs); public: DipoleMap::iterator begin() { return stringdipoles.begin(); } DipoleMap::iterator end() { return stringdipoles.end(); } DipoleMap::iterator beginSpectators() { return spectatorDipoles.begin(); } DipoleMap::iterator endSpectators() { return spectatorDipoles.end(); } /** * Exception class. */ struct ColourException: public Exception {}; private: /** * The assumed radius of a string. */ Length R0; /** * The assumed mass for calculating rapidities */ Energy m0; /* * The junction colour fluctuation parameter */ double junctionDiquarkProb; /* * Cutoff distance where we no longer consider * shoving */ Length rCutOff; /* * Amplitude af shoving function */ StringTension gAmplitude; /* * Shoving function (inverse) multiplier */ double gExponent; /* * Cutoff in effective rapidity span when clustering */ double yCutOff; /* * Rapidity slicing for shoving */ double deltay; /* * Shoving time */ Length tshove; /* * Shoving time slicing */ Length deltat; /** * The vector of all Dipoles */ vector dipoles; /** * All participating Dipoles in all string */ DipoleMap stringdipoles; /* * Spectating dipoles */ DipoleMap spectatorDipoles; /* * Lorentz contraction of fields */ bool lorentz; bool cylinderShove; // For keeping only soft long strings Energy stringpTCut; double minStringy; bool longSoft; // For analysis HistKeeper* hkPtr; /** * Debug flag. */ bool debug; public: mutable double strl0; mutable double strl; mutable double avh; mutable double avp; mutable double avq; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ Ropewalk & operator=(const Ropewalk &); double nb; }; } #endif /* TheP8I_Ropewalk_H */ diff --git a/Hadronization/StringFragmentation.cc b/Hadronization/StringFragmentation.cc --- a/Hadronization/StringFragmentation.cc +++ b/Hadronization/StringFragmentation.cc @@ -1,1101 +1,1101 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the StringFragmentation class. // #include "StringFragmentation.h" #include "Ropewalk.h" #include "RandomAverageHandler.h" #include "RandomHandler.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/PDT/RemnantDecayer.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/EventRecord/Step.h" #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DebugItem.h" #include "ThePEG/Utilities/EnumIO.h" #include "ThePEG/Utilities/MaxCmp.h" #include "ThePEG/Utilities/HoldFlag.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include using namespace TheP8I; StringFragmentation::StringFragmentation() : pythia(), fScheme(0), stringR0(0.5*femtometer), stringm0(1.0*GeV), junctionDiquark(1.0), alpha(1.0), average(true), stringyCut(2.0), stringpTCut(6*GeV), fragmentationMass(0.134), baryonSuppression(0.5), window(false), minStringy(8.0), longSoft(false), longStringsOnly(false), throwaway(false), shoving(false), deltay(5.0), tshove(1.0*femtometer), deltat(0.1*femtometer), rCutOff(4.0), gAmplitude(1.0*GeV/femtometer), gExponent(1.0), yCutOff(4.0), lorentz(true), cylinder(false), _eventShapes(0), useThrustAxis(0), opHandler(0), maxTries(2), doAnalysis(0), analysisPath(""), #include "StringFragmentation-init.h" { } StringFragmentation::~StringFragmentation() { if ( opHandler ) delete opHandler; } void StringFragmentation::handle(EventHandler & eh, const tPVector & tagged, const Hint & h) { ++nev; static unsigned long lastEventId = 0; bool secondary = false; // Get the event tPVector all = RemnantDecayer::decayRemnants(tagged, *newStep()); HoldFlag oldFS(fScheme, fScheme); // Prevent funny behavior if hadronizing decays of eg. Upsilon. if ( newStep()->collision()->uniqueId == lastEventId ){ secondary = true; if(fScheme == 4 || fScheme == 5 || fScheme == 1) fScheme = 99; else fScheme = 0; } else lastEventId = newStep()->collision()->uniqueId; if(_eventShapes) _eventShapes->reset(all); vector singlets; singlets.clear(); if ( theCollapser && !secondary ) singlets = theCollapser->collapse(all, newStep()); else singlets = ColourSinglet::getSinglets(all.begin(), all.end()); // Goto correct hadronization scheme switch(fScheme){ case 0: // Pythia { hadronizeSystems(*opHandler->GetPythiaPtr(1.0,nev%10000==0), singlets, all); break; } case 1: // For playing around and printing stuff { // Cartoon of string shoving // if(singlets.size() > 10){ Ropewalk ropewalkInit(singlets, stringR0, stringm0, baryonSuppression, throwaway, false); vector< pair > allStrings = ropewalkInit.getSinglets(stringyCut); tPVector allParticles; for ( vector< pair >::iterator sItr = allStrings.begin(); sItr != allStrings.end(); ++sItr ) for( tcPVector::iterator pItr = sItr->first.partons().begin(); pItr != sItr->first.partons().end(); ++pItr ){ allParticles.push_back(const_ptr_cast(*pItr)); } singlets = theCollapser->collapse(allParticles, newStep()); for(double time = 0.0; time < 1.0; time += 0.1){ Ropewalk ropewalkOne(singlets, stringR0, stringm0, baryonSuppression, false, shoving, rCutOff, gAmplitude, gExponent, yCutOff, deltay, time*femtometer, 0.1*femtometer, lorentz, cylinder, stringpTCut, minStringy, longSoft, NULL, false); auto css = ropewalkOne.getSinglets(1.0); vector cs; for(auto itr = css.begin(); itr != css.end(); ++itr) cs.push_back(itr->first); ofstream out((analysisPath + "../lhcevent.m").c_str(),ios::app); out << PrintStringsToMatlab(cs) << endl; out << "figure;" << endl; out << "tubeplot(a,radius);" << endl; out << "hold on; tubeplot(b,radius);" << endl; out << "hold on; tubeplot(c,radius);" << endl;; out << "hold on; tubeplot(d,radius);" << endl;; out << "hold on; tubeplot(e,radius);" << endl;; out << "hold on; tubeplot(f,radius);" << endl;; out << "hold on; tubeplot(g,radius);" << endl;; out << "hold on; tubeplot(h,radius);" << endl;; out << "hold on; tubeplot(i,radius);" << endl;; out << "hold on; tubeplot(j,radius);" << endl;; out << "hold on; tubeplot(k,radius);" << endl;; out << "hold on; tubeplot(l,radius);" << endl;; out << "hold on; tubeplot(m,radius);" << endl;; } exit(1); } /*if(throwaway){ Ropewalk ropewalk_init(singlets, stringR0, stringm0, baryonSuppression, throwaway, false); vector< pair > allStrings = ropewalk_init.getSinglets(stringyCut); tPVector allParticles; for(vector< pair >::iterator sItr = allStrings.begin(); sItr != allStrings.end(); ++sItr) for(tcPVector::iterator pItr = sItr->first.partons().begin(); pItr != sItr->first.partons().end(); ++pItr) allParticles.push_back(const_ptr_cast(*pItr)); singlets = theCollapser->collapse(allParticles, newStep()); } Ropewalk ropewalk(singlets, stringR0, stringm0, baryonSuppression, false, false); vector< pair > strings = ropewalk.getSinglets(stringyCut); static ofstream out(( generator()->runName() + ".txt").c_str()); if(!std::isnan(ropewalk.lambdaSum()+ropewalk.getkW()+ropewalk.getNb())) out << generator()->currentEvent()->weight() << " " << ropewalk.lambdaSum() << " " << ropewalk.getkW() << " " << ropewalk.getNb() << endl; //vector mspec = ropewalk.mspec; //for(int i = 0, N = mspec.size(); i < N; ++i) out << generator()->currentEvent()->weight() << " " << mspec[i] << endl; */ break; } case 2: // Do shoving and construct parton level histograms. { double weight = generator()->currentEvent()->weight(); if(doAnalysis) hk.setWeight(weight); if (throwaway) { Ropewalk ropewalkInit(singlets, stringR0, stringm0, baryonSuppression, throwaway, false); if(doAnalysis){ hk.hist("stringLengthBefore")->fill(ropewalkInit.lambdaSum(),weight); hk.hist("stringLengthNoCutBefore")->fill(ropewalkInit.lambdaSum(0.0),weight); } vector< pair > allStrings = ropewalkInit.getSinglets(stringyCut); tPVector allParticles; for ( vector< pair >::iterator sItr = allStrings.begin(); sItr != allStrings.end(); ++sItr ) for( tcPVector::iterator pItr = sItr->first.partons().begin(); pItr != sItr->first.partons().end(); ++pItr ){ allParticles.push_back(const_ptr_cast(*pItr)); } singlets = theCollapser->collapse(allParticles, newStep()); } if(doAnalysis){ // string length before shoving Ropewalk rr(singlets,stringR0,stringm0,baryonSuppression,false,false); hk.hist("stringLengthCollapse")->fill(rr.lambdaSum(),weight); hk.hist("stringLengthNoCutCollapse")->fill(rr.lambdaSum(0.0),weight); } typedef multimap OrderedMap; OrderedMap ordered; Ropewalk ropewalk(singlets, stringR0, stringm0, baryonSuppression, false, shoving, rCutOff, gAmplitude, gExponent, yCutOff, deltay, tshove, deltat, lorentz, cylinder, stringpTCut, minStringy, longSoft, (doAnalysis ? &hk : NULL), false); if(doAnalysis){ hk.hist("stringLengthShove")->fill(ropewalk.lambdaSum(),weight); hk.hist("stringLengthNoCutShove")->fill(ropewalk.lambdaSum(0.0),weight); } for(Ropewalk::DipoleMap::iterator itr = ropewalk.begin(); itr != ropewalk.end(); ++itr) ordered.insert(make_pair(maxPT(*itr->first, stringyCut), itr)); for(OrderedMap::iterator it = ordered.begin(); it != ordered.end(); ++it ) { if(doAnalysis){ hk.hist("averageKappa")->fill(ropewalk.averageKappaEnhancement(it->second, stringyCut)); } if(!pythia.getRopeUserHooksPtr()->setDipoles(&(*it->second), stringm0, stringR0, !average, alpha)) { pythia.getRopeUserHooksPtr()->setEnhancement(ropewalk.averageKappaEnhancement(it->second, stringyCut)); for (int i = 0, N = it->second->second.size(); i < N; ++i) it->second->second[i]->hadr = true; } vector toHadronization(1, *(it->second->first)); forcerun = false; if(!hadronizeSystems(pythia, toHadronization, all)) hadronizeSystems(pythia, toHadronization, all); } break; } case 3: // Pipes with average { RandomAverageHandler avg_enhancer(throwaway); // TODO: Implement throwaway functionality in this avg_enhancer.clear(); vector pipes; // Make all the pipes for(vector::iterator sItr = singlets.begin(); sItr!=singlets.end(); ++sItr) pipes.push_back(StringPipe(&(*sItr),stringR0,_eventShapes)); avg_enhancer.SetEvent(pipes); for(vector::iterator it = pipes.begin(); it!=pipes.end(); ++it){ vector toHadronization; double h = avg_enhancer.KappaEnhancement(*it); if(h > 0){ toHadronization.push_back(*(*it).GetSingletPtr()); forcerun = false; if(!hadronizeSystems(*(opHandler->GetPythiaPtr(h,nev%10000==0)),toHadronization,all)) hadronizeSystems(*(opHandler->GetPythiaPtr(1.0,false)),toHadronization,all); } } break; } case 4: // Average kappa over whole strings { Ropewalk ropewalk(singlets, stringR0, stringm0, baryonSuppression,throwaway, false); vector< pair > strings = ropewalk.getSinglets(stringyCut); vector toHadronization; double avge = 0; for(int i = 0, N = strings.size(); i < N; ++i ){ avge += strings[i].second; pythia.getRopeUserHooksPtr()->setEnhancement(strings[i].second); toHadronization.clear(); toHadronization.push_back(strings[i].first); hadronizeSystems(pythia,toHadronization,all); } // ofstream out(( "/scratch/galette/bierlich/work/dipsy/pprange/" + generator()->runName() + "/stringplot.txt").c_str(),ios::app); // out << generator()->currentEvent()->weight() << " " << strings.size() << " " << avge << endl; break; } case 5: // Altering kappa per dipole { double weight = generator()->currentEvent()->weight(); if(doAnalysis) hk.setWeight(weight); if (throwaway) { Ropewalk ropewalkInit(singlets, stringR0, stringm0, baryonSuppression, throwaway, false); if(doAnalysis){ hk.hist("stringLengthBefore")->fill(ropewalkInit.lambdaSum(),weight); hk.hist("stringLengthNoCutBefore")->fill(ropewalkInit.lambdaSum(0.0),weight); } vector< pair > allStrings = ropewalkInit.getSinglets(stringyCut); tPVector allParticles; for ( vector< pair >::iterator sItr = allStrings.begin(); sItr != allStrings.end(); ++sItr ) for( tcPVector::iterator pItr = sItr->first.partons().begin(); pItr != sItr->first.partons().end(); ++pItr ){ allParticles.push_back(const_ptr_cast(*pItr)); } singlets = theCollapser->collapse(allParticles, newStep()); } if(doAnalysis){ // string length before shoving Ropewalk rr(singlets,stringR0,stringm0,baryonSuppression,false,false); hk.hist("stringLengthCollapse")->fill(rr.lambdaSum(),weight); hk.hist("stringLengthNoCutCollapse")->fill(rr.lambdaSum(0.0),weight); } typedef multimap OrderedMap; OrderedMap ordered; Ropewalk ropewalk(singlets, stringR0, stringm0, baryonSuppression, false, shoving, rCutOff, gAmplitude, gExponent, yCutOff, deltay, tshove, deltat, lorentz, cylinder, stringpTCut, minStringy, longSoft, (doAnalysis ? &hk : NULL), false); if(doAnalysis){ hk.hist("stringLengthShove")->fill(ropewalk.lambdaSum(),weight); hk.hist("stringLengthNoCutShove")->fill(ropewalk.lambdaSum(0.0),weight); } for(Ropewalk::DipoleMap::iterator itr = ropewalk.begin(); itr != ropewalk.end(); ++itr) ordered.insert(make_pair(maxPT(*itr->first, stringyCut), itr)); for(OrderedMap::iterator it = ordered.begin(); it != ordered.end(); ++it ) { if(doAnalysis){ hk.hist("averageKappa")->fill(ropewalk.averageKappaEnhancement(it->second, stringyCut)); } else if(!pythia.getRopeUserHooksPtr()->setDipoles(&(*it->second), stringm0, stringR0, !average, alpha)) { pythia.getRopeUserHooksPtr()->setEnhancement(ropewalk.averageKappaEnhancement(it->second, stringyCut)); for (int i = 0, N = it->second->second.size(); i < N; ++i) it->second->second[i]->hadr = true; } vector toHadronization(1, *(it->second->first)); forcerun = false; if(!hadronizeSystems(pythia, toHadronization, all)) hadronizeSystems(pythia, toHadronization, all); } if(!longStringsOnly){ ordered.clear(); for(Ropewalk::DipoleMap::iterator itr = ropewalk.beginSpectators(); itr != ropewalk.endSpectators(); ++itr) ordered.insert(make_pair(maxPT(*itr->first, stringyCut), itr)); for(OrderedMap::iterator it = ordered.begin(); it != ordered.end(); ++it ) { if(doAnalysis){ hk.hist("averageKappa")->fill(ropewalk.averageKappaEnhancement(it->second, stringyCut)); } pythia.getRopeUserHooksPtr()->setEnhancement(1.0); vector toHadronization(1, *(it->second->first)); forcerun = false; if(!hadronizeSystems(pythia, toHadronization, all)) hadronizeSystems(pythia, toHadronization, all); } } break; } case 99: // New Pythia { pythia.getRopeUserHooksPtr()->setEnhancement(1.0); hadronizeSystems(pythia, singlets, all); break; } default: cout << "We should really not be here. This is bad..." << endl; } // fScheme = oldFS; // Do short analysis here: if(doAnalysis){ ofstream out((analysisPath + "analysis.txt").c_str(),ios::app); // out << "" << endl; // out << PrintStringsToMatlab(singlets) << endl; // out << "" << endl; out.close(); } } void StringFragmentation::dofinish(){ if(doAnalysis){ ofstream of((analysisPath + generator()->runName()) + ".histo"); of << hk; of.close(); } //if(doAnalysis!=0){ // YODA::mkWriter("yoda"); // YODA::WriterYODA::write(analysisPath + generator()->runName() + "1D.yoda",_histograms.begin(),_histograms.end()); // YODA::WriterYODA::write(analysisPath + generator()->runName() + "2D.yoda",_histograms2D.begin(),_histograms2D.end()); //} } Energy StringFragmentation::maxPT(const ColourSinglet & cs, double deltaY) { MaxCmp maxpt2; for ( int i = 0, N = cs.partons().size(); i < N; ++i ) { if ( deltaY > 0.0 && abs(cs.partons()[i]->eta()) > deltaY ) continue; maxpt2(cs.partons()[i]->momentum().perp2()); } if ( !maxpt2 ) return ZERO; return sqrt(maxpt2.value()); } bool StringFragmentation:: hadronizeSystems(Pythia8Interface & pyt, const vector & singlets, const tPVector & all) { TheP8EventShapes * _es = NULL; Pythia8::Event & event = pyt.event(); pyt.clearEvent(); for ( int i = 0, N = singlets.size(); i < N; ++i ) { if ( singlets[i].nPieces() == 3 ) { // Simple junction. // Save place where we will store dummy particle. int nsave = event.size(); event.append(22, -21, 0, 0, 0, 0, 0, 0, 0.0, 0.0, 0.0, 0.0); int npart = 0; for ( int ip = 1; ip <= 3; ++ip ) for ( int j = 0, M = singlets[i].piece(ip).size(); j < M; ++j ) { pyt.addParticle(singlets[i].piece(ip)[j], 23, nsave, 0); ++npart; } event[nsave].daughters(nsave + 1, nsave + npart); } else if ( singlets[i].nPieces() == 5 ) { // Double, connected junction // Save place where we will store dummy beam particles. int nb1 = event.size(); int nb2 = nb1 + 1; event.append(2212, -21, 0, 0, nb1 + 2, nb1 + 4, 0, 0, 0.0, 0.0, 0.0, 0.0); event.append(-2212, -21, 0, 0, nb2 + 4, nb2 + 6, 0, 0, 0.0, 0.0, 0.0, 0.0); // Find the string piece connecting the junctions, and the // other loose string pieces. int connector = 0; int q1 = 0; int q2 = 0; int aq1 = 0; int aq2 = 0; for ( int ip = 1; ip <= 5; ++ip ) { if ( singlets[i].sink(ip).first && singlets[i].source(ip).first ) connector = ip; else if ( singlets[i].source(ip).first ) { if ( q1 ) q2 = ip; else q1 = ip; } else if ( singlets[i].sink(ip).first ) { if ( aq1 ) aq2 = ip; else aq1 = ip; } } if ( !connector || !q1 || !q2 || !aq1 || ! aq2 ) Throw() << name() << " found complicated junction string. Although Pythia8 can " << "hadronize junction strings, this one was too complicated." << Exception::runerror; // Insert the partons of the loose triplet ends. int start = event.size(); for ( int j = 0, M = singlets[i].piece(q1).size(); j < M; ++j ) pyt.addParticle(singlets[i].piece(q1)[j], 23, nb1, 0); for ( int j = 0, M = singlets[i].piece(q2).size(); j < M; ++j ) pyt.addParticle(singlets[i].piece(q2)[j], 23, nb1, 0); // Insert dummy triplet incoming parton with correct colour code. int col1 = singlets[i].piece(connector).empty()? event.nextColTag(): pyt.addColourLine(singlets[i].piece(connector).front()->colourLine()); int dum1 = event.size(); event.append(2, -21, nb1, 0, 0, 0, col1, 0, 0.0, 0.0, 0.0, 0.0, 0.0); event[nb1].daughters(start, start + singlets[i].piece(q1).size() + singlets[i].piece(q2).size()); // Insert the partons of the loose anti-triplet ends. start = event.size(); for ( int j = 0, M = singlets[i].piece(aq1).size(); j < M; ++j ) pyt.addParticle(singlets[i].piece(aq1)[j], 23, nb2, 0); for ( int j = 0, M = singlets[i].piece(aq2).size(); j < M; ++j ) pyt.addParticle(singlets[i].piece(aq2)[j], 23, nb2, 0); // Insert dummy anti-triplet incoming parton with correct colour code. int col2 = singlets[i].piece(connector).empty()? col1: pyt.addColourLine(singlets[i].piece(connector).back()->antiColourLine()); int dum2 = event.size(); event.append(-2, -21, nb2, 0, 0, 0, 0, col2, 0.0, 0.0, 0.0, 0.0, 0.0); event[nb2].daughters(start, start + singlets[i].piece(aq1).size() + singlets[i].piece(aq2).size()); // Add the partons from the connecting string piece. for ( int j = 0, M = singlets[i].piece(connector).size(); j < M; ++j ) pyt.addParticle(singlets[i].piece(connector)[j], 23, dum1, dum2); } else if ( singlets[i].nPieces() > 1 ) { // We don't know how to handle other junctions yet. Throw() << name() << " found complicated junction string. Although Pythia8 can " << "hadronize junction strings, that interface is not ready yet." << Exception::runerror; } else { // Normal string for ( int j = 0, M = singlets[i].partons().size(); j < M; ++j ) pyt.addParticle(singlets[i].partons()[j], 23, 0, 0); } } for ( int i = 0, N = all.size(); i < N; ++i ) if ( !all[i]->coloured() ) pyt.addParticle(all[i], 23, 0, 0); int oldsize = event.size(); Pythia8::Event saveEvent = event; int ntry = maxTries; CurrentGenerator::Redirect stdout(cout, false); while ( !pyt.go() && --ntry ) event = saveEvent; //event.list(cout); //event.list(cerr); if ( !ntry ) { static DebugItem printfailed("TheP8I::PrintFailed", 10); if ( printfailed ) { double ymax = -1000.0; double ymin = 1000.0; double sumdy = 0.0; for ( int i = 0, N = singlets.size(); i < N; ++i ) for ( int j = 0, M = singlets[i].partons().size(); j < M; ++j ) { const Particle & p = *singlets[i].partons()[j]; ymax = max(ymax, (_es ? _es->yT(p.momentum()) : p.momentum().rapidity())); //ymax = max(ymax, p.momentum().rapidity()); ymin = min(ymin,(_es ? _es->yT(p.momentum()) : p.momentum().rapidity())); //ymin = min(ymin, p.momentum().rapidity()); cerr << setw(5) << j << setw(14) << p.momentum().rapidity() << setw(14) << p.momentum().perp()/GeV << setw(14) << p.momentum().m()/GeV; if ( j == 0 && p.data().iColour() == PDT::Colour8 ) { cerr << setw(14) << (p.momentum() + singlets[i].partons().back()->momentum()).m()/GeV; sumdy += abs((_es ? _es->yT(p.momentum()) : p.momentum().rapidity()) ) -(_es ? _es->yT(singlets[i].partons().back()->momentum()) : singlets[i].partons().back()->momentum().rapidity()); //sumdy += abs(p.momentum().rapidity() - singlets[i].partons().back()->momentum().rapidity()); } else if ( j > 0 ) { cerr << setw(14) << (p.momentum() + singlets[i].partons()[j-1]->momentum()).m()/GeV; sumdy += abs((_es ? _es->yT(p.momentum()) : p.momentum().rapidity()) ) -(_es ? _es->yT(singlets[i].partons()[j-i]->momentum()) : singlets[i].partons()[j-i]->momentum().rapidity()); //sumdy += abs(p.momentum().rapidity() - singlets[i].partons()[j-1]->momentum().rapidity()); if ( j + 1 < M ) cerr << setw(14) << (p.momentum() + singlets[i].partons()[j-1]->momentum()).m()* (p.momentum() + singlets[i].partons()[j+1]->momentum()).m()/ (p.momentum() + singlets[i].partons()[j+1]->momentum() + singlets[i].partons()[j-1]->momentum()).m()/GeV; } cerr << endl; } cerr << setw(14) << ymax - ymin << setw(14) << sumdy/(ymax - ymin) << endl; } CurrentGenerator::Redirect stdout2(cout, true); //event.list(); pyt.errorlist(); cout << "ThePEG event listing:\n" << *(generator()->currentEvent()); Throw() << "Pythia8 failed to hadronize partonic state:\n" << stdout2.str() << "This event will be discarded!\n" << Exception::warning; throw Veto(); } if(window){ // cout << stringyCut << " " << stringpTCut/GeV << endl; if(forcerun) forcerun = false; else{ for( int i = 1, N = event.size(); i < N; ++i ) { tPPtr p = pyt.getParticle(i); // cout << p->momentum().perp()/GeV << " " << p->rapidity() << " " << p->id() << endl; if (p->momentum().perp() > stringpTCut && abs(p->rapidity()) < stringyCut ){ /* if(abs(p->id()) == 310 || abs(p->id()) == 3122){ static ofstream fout(( generator()->runName() + ".txt").c_str(),ios::app); tcPVector pVector = singlets[0].partons(); fout << p->id() << " " << generator()->currentEvent()->weight() << " "; for(size_t j=0;jid()<20){ fout << pVector[j]->id() << " " << pVector[j]->momentum().perp()/GeV; fout.close(); //} } */ forcerun = true; return false; } } } } map > children; map > parents; for ( int i = 1, N = event.size(); i < N; ++i ) { tPPtr p = pyt.getParticle(i); int d1 = event[i].daughter1(); if ( d1 <= 0 ) continue; children[p].insert(pyt.getParticle(d1)); parents[pyt.getParticle(d1)].insert(p); int d2 = event[i].daughter2(); if ( d2 > 0 ) { children[p].insert(pyt.getParticle(d2)); parents[pyt.getParticle(d2)].insert(p); } for ( int di = d1 + 1; di < d2; ++di ) { children[p].insert(pyt.getParticle(di)); parents[pyt.getParticle(di)].insert(p); } } for ( int i = oldsize, N = event.size(); i < N; ++i ) { PPtr p = pyt.getParticle(i); set & pars = parents[p]; if ( !p ) { - event.list(cout); + event.list(); cout << i << endl; Throw() << "Failed to reconstruct hadronized state from Pythia8:\n" << stdout.str() << "This event will be discarded!\n" << Exception::warning; throw Veto(); } if ( std::isnan((p->momentum().perp()/GeV)) ){ Throw() << "Failed to reconstruct hadronized state from Pythia8:\n" << stdout.str() << "This event will be discarded!\n" << Exception::warning; throw Veto(); } if ( pars.empty() ) { Pythia8::Particle & pyp = event[i]; if ( pyp.mother1() > 0 ) pars.insert(pyt.getParticle(pyp.mother1())); if ( pyp.mother2() > 0 ) pars.insert(pyt.getParticle(pyp.mother2())); for ( int im = pyp.mother1() + 1; im < pyp.mother2(); ++im ) pars.insert(pyt.getParticle(im)); if ( pars.empty() ) { Throw() << "Failed to reconstruct hadronized state from Pythia8:\n" << stdout.str() << "This event will be discarded!\n" << Exception::warning; throw Veto(); } } if ( pars.size() == 1 ) { tPPtr par = *pars.begin(); if ( children[par].size() == 1 && *children[par].begin() == p && par->id() == p->id() ) newStep()->setCopy(par, p); else newStep()->addDecayProduct(par, p); } else { newStep()->addDecayProduct(pars.begin(), pars.end(), p); } } return true; } string StringFragmentation::PrintStringsToMatlab(vector& singlets) { stringstream ss; vector > drawit; vector names; for(int i=0;i<26;++i) names.push_back(char(i+97)); vector::iterator nItr = names.begin(); for(vector::iterator sItr = singlets.begin(); sItr!=singlets.end(); ++sItr){ drawit.clear(); for (tcPVector::iterator pItr = sItr->partons().begin(); pItr!=sItr->partons().end();++pItr) { vector tmp; tmp.clear(); tmp.push_back((*pItr)->rapidity()); tmp.push_back((*pItr)->vertex().x()/femtometer); tmp.push_back((*pItr)->vertex().y()/femtometer); drawit.push_back(tmp); } ss << *nItr << " = ["; for(int i=0, N=drawit.size();i 0 && pythia.version() - 1.234 > 1.0){ // cout << "Pythia version: " << pythia.version() << endl; // cout << "The chosen fragmentation scheeme requires a tweaked " // << "Pythia v. 1.234 for option. I will default you to " // "option 'pythia'." << endl; // fScheme = 0; // } // else{ PytPars p; p.insert(make_pair("StringPT:sigma",theStringPT_sigma)); p.insert(make_pair("StringZ:aLund",theStringZ_aLund)); p.insert(make_pair("StringZ:bLund",theStringZ_bLund)); p.insert(make_pair("StringFlav:probStoUD",theStringFlav_probStoUD)); p.insert(make_pair("StringFlav:probSQtoQQ",theStringFlav_probSQtoQQ)); p.insert(make_pair("StringFlav:probQQ1toQQ0",theStringFlav_probQQ1toQQ0)); p.insert(make_pair("StringFlav:probQQtoQ",theStringFlav_probQQtoQ)); phandler.init(fragmentationMass*fragmentationMass,baryonSuppression,p); pythia.getRopeUserHooksPtr()->setParameterHandler(&phandler); //pythia.getRopeUserHooksPtr()->setWindow(window,stringpTCut); // } } if( !( fScheme == 4 || fScheme == 5 || fScheme == 1) ){ // Not 'else' as it can change above... vector moresettings = theAdditionalP8Settings; moresettings.push_back("StringPT:sigma = " + convert(theStringPT_sigma)); moresettings.push_back("StringZ:aLund = " + convert(theStringZ_aLund)); moresettings.push_back("StringZ:bLund = " + convert(theStringZ_bLund)); moresettings.push_back("StringFlav:probStoUD = " + convert(theStringFlav_probStoUD)); moresettings.push_back("StringFlav:probSQtoQQ = " + convert(theStringFlav_probSQtoQQ)); moresettings.push_back("StringFlav:probQQ1toQQ0 = " + convert(theStringFlav_probQQ1toQQ0)); moresettings.push_back("StringFlav:probQQtoQ = " + convert(theStringFlav_probQQtoQ)); moresettings.push_back("OverlapStrings:fragMass = " + convert(fragmentationMass)); moresettings.push_back("OverlapStrings:baryonSuppression = " + convert(baryonSuppression)); // Initialize the OverlapPythia Handler if ( opHandler ) delete opHandler; opHandler = new OverlapPythiaHandler(this,moresettings); } // Should we do event shapes? if(_eventShapes) delete _eventShapes; _eventShapes = NULL; if(useThrustAxis==1){ _eventShapes = new TheP8EventShapes(); } //if( doAnalysis ) { // _histograms.push_back(new YODA::Histo1D(100,1,15,"/h_lowpT","h_lowpT")); //_histograms.push_back(new YODA::Histo1D(100,1,15,"/h_midpT","h_midpT")); // _histograms.push_back(new YODA::Histo1D(100,1,15,"/h_highpT","h_highpT")); // _histograms2D.push_back(new YODA::Histo2D(10,0.,15,10,0,15,"/m_vs_pT","m_vs_pT")); //} nev = 0; } void StringFragmentation::persistentOutput(PersistentOStream & os) const { os #include "StringFragmentation-output.h" << fScheme << ounit(stringR0,femtometer) << ounit(stringm0,GeV) << junctionDiquark << alpha << average << stringyCut << ounit(stringpTCut,GeV) << fragmentationMass << baryonSuppression << oenum(window) << minStringy << oenum(longSoft) << oenum(longStringsOnly) << oenum(throwaway) << oenum(shoving) << deltay << ounit(tshove,femtometer) << ounit(deltat,femtometer) << rCutOff << ounit(gAmplitude,GeV/femtometer) << gExponent << yCutOff << oenum(lorentz) << oenum(cylinder) << useThrustAxis << maxTries << theCollapser << doAnalysis << analysisPath; } void StringFragmentation::persistentInput(PersistentIStream & is, int) { is #include "StringFragmentation-input.h" >> fScheme >> iunit(stringR0,femtometer) >> iunit(stringm0,GeV) >> junctionDiquark >> alpha >> average >> stringyCut >> iunit(stringpTCut,GeV) >> fragmentationMass >> baryonSuppression >> ienum(window) >> minStringy >> ienum(longSoft) >> ienum(longStringsOnly) >> ienum(throwaway) >> ienum(shoving) >> deltay >> iunit(tshove,femtometer) >> iunit(deltat,femtometer) >> rCutOff >> iunit(gAmplitude,GeV/femtometer) >> gExponent >> yCutOff >> ienum(lorentz) >> ienum(cylinder) >> useThrustAxis >> maxTries >> theCollapser >> doAnalysis >> analysisPath; } ClassDescription StringFragmentation::initStringFragmentation; // Definition of the static class description member. void StringFragmentation::Init() { #include "StringFragmentation-interfaces.h" static Reference interfaceCollapser ("Collapser", "A ThePEG::ClusterCollapser object used to collapse colour singlet " "clusters which are too small to fragment. If no object is given the " "MinistringFragmentetion of Pythia8 is used instead.", &StringFragmentation::theCollapser, true, false, true, true, false); static Switch interfaceFragmentationScheme ("FragmentationScheme", "Different options for how to handle overlapping strings.", &StringFragmentation::fScheme, 0, true, false); static SwitchOption interfaceFragmentationSchemepythia (interfaceFragmentationScheme, "pythia", "Plain old Pythia fragmentation", 0); static SwitchOption interfaceFragmentationSchemedep1 (interfaceFragmentationScheme, "dep1", "Not sure about this.", 1); static SwitchOption interfaceFragmentationSchemedep2 (interfaceFragmentationScheme, "dep2", "Not sure about this.", 2); static SwitchOption interfaceFragmentationSchemenone (interfaceFragmentationScheme, "none", "Plain old Pythia fragmentation", 0); static SwitchOption interfaceFragmentationSchemepipe (interfaceFragmentationScheme, "pipe", "Each string is enclosed by a cylinder. Effective parameters are calculated from the overlap of cylinders.", 3); static SwitchOption interfaceFragmentationSchemeaverage (interfaceFragmentationScheme, "average", "The overlap is calculated for each dipole in all strings, and for each string the effective parameters are obtained from the average overlap.", 4); static SwitchOption interfaceFragmentationSchemedipole (interfaceFragmentationScheme, "dipole", "Effective parameters are calculated for each breakup by determining the overlap of the corresponding dipole with other dipoles.", 5); static Parameter interfaceUseThrustAxis ("UseThrustAxis", "Whether or not to use rapidity wrt. Thrust Axis when counting strings " "(put 1 or 0).", &StringFragmentation::useThrustAxis, 0, 0, 0, true, false, Interface::lowerlim); static Parameter interfaceDoAnalysis ("DoAnalysis", "Can do a short analysis where histograms of the effective parameters are " "plotted, and the strings in (bx,by,y)-space are printed for a couple of events. " "(put 1 or 0).", &StringFragmentation::doAnalysis, 0, 0, 0, true, false, Interface::lowerlim); static Parameter interfaceAnalysisPath ("AnalysisPath", "Set the system path where the (optional) analysis output should appear " " (directory must exist!)" , &StringFragmentation::analysisPath, ""); static Switch interfaceThrowAway ("ThrowAway", "Throw away strings with weight:" "Use 1 - (p + q)/(m + n) for (0,-1) configuration." , &StringFragmentation::throwaway, false, true, false); static SwitchOption interfaceThrowAwayTrue (interfaceThrowAway,"True","enabled.",true); static SwitchOption interfaceThrowAwayFalse (interfaceThrowAway,"False","disabled.",false); static Switch interfaceShoving ("Shoving", "Strings will shove each other around.", &StringFragmentation::shoving, false, true, false); static SwitchOption interfaceShovingTrue (interfaceShoving,"True","enabled.",true); static SwitchOption interfaceShovingFalse (interfaceShoving,"False","disabled.",false); static Parameter interfaceDeltaY ("DeltaY", "The size of rapidity slicing intervals", &StringFragmentation::deltay, 0, 0.1, 0.0, 0, true, false, Interface::lowerlim); static Parameter interfaceTShove ("TShove", "The total shoving -and propagation time", &StringFragmentation::tshove, femtometer, 1.0*femtometer, 0.0*femtometer, 0*femtometer, true, false, Interface::lowerlim); static Parameter interfaceDeltaT ("DeltaT", "Shoving time slicing -- size of a slice", &StringFragmentation::deltat, femtometer, 0.1*femtometer, 0.0*femtometer, 0*femtometer, true, false, Interface::lowerlim); static Parameter interfaceRCutOff ("RCutOff", "Cutoff radius at which we dont let strings shove each" "other around anymore -- given as a multiple of StringR0", &StringFragmentation::rCutOff, 0, 4.0, 0.0, 0, true, false, Interface::lowerlim); static Parameter interfaceGAmplitude ("GAmplitude", "Amplitude of shoving energy density function, should be approx." "the average effective string tension", &StringFragmentation::gAmplitude, GeV/femtometer, 1.0*GeV/femtometer, 0.0*GeV/femtometer, 0*GeV/femtometer, true, false, Interface::lowerlim); static Parameter interfaceGExponent ("GExponent", "Parameter to divide the exponent of energy density function for" "shoving.", &StringFragmentation::gExponent, 0, 10.0, 0.0, 0, true, false, Interface::lowerlim); static Parameter interfaceYCutOff ("YCutOff", "Cutoff rapidity for excitation clustering. This is the minimal" "effective rapidity span for a dipole after clustering", &StringFragmentation::yCutOff, 0, 5.0, 0.0, 0, true, false, Interface::lowerlim); static Switch interfaceLorentz ("Lorentz", "Do Lorentz contraction of string fields moving with a transverse velocity", &StringFragmentation::lorentz, true, true, false); static SwitchOption interfaceLorentzTrue (interfaceLorentz,"True","enabled.",true); static SwitchOption interfaceLorentzFalse (interfaceLorentz,"False","disabled.",false); static Switch interfaceCylinder ("Cylinder", "Do Cylinder Area Shoving", &StringFragmentation::cylinder, false, true, false); static SwitchOption interfaceCylinderTrue (interfaceCylinder,"True","enabled.",true); static SwitchOption interfaceCylinderFalse (interfaceCylinder,"False","disabled.",false); static Switch interfaceWindow ("StringWindow", "Enable the 'window'-cut procedure, off by default." "Parameters will only affect if this is switched on. " , &StringFragmentation::window, false, true, false); static SwitchOption interfaceWindowTrue (interfaceWindow,"True","enabled.",true); static SwitchOption interfaceWindowFalse (interfaceWindow,"False","disabled.",false); static Parameter interfaceStringpTCut ("StringpTCut", "Strings with a constituent pT higher than StringpTCut (GeV), will not participate " " in the ropewalk", &StringFragmentation::stringpTCut, GeV, 3.0*GeV, 0.0*GeV, 0*GeV, true, false, Interface::lowerlim); static Parameter interfaceMinStringy ("MinStringy", "String smaller than MinStringy (in rapidity) will not participate " " in the ropewalk", &StringFragmentation::minStringy, 0, 8.0, 0.0, 0, true, false, Interface::lowerlim); static Switch interfaceLongSoft ("LongSoft", "Only do ropewalk with long and soft strings, set cuts accordingly" , &StringFragmentation::longSoft, false, true, false); static SwitchOption interfaceLongSoftTrue (interfaceLongSoft,"True","enabled.",true); static SwitchOption interfaceLongSoftFalse (interfaceLongSoft,"False","disabled.",false); static Switch interfaceLongStringsOnly ("LongStringsOnly", "Hadronize only long and soft strings, no spectators, set cuts accordingly" , &StringFragmentation::longStringsOnly, false, true, false); static SwitchOption interfaceLongStringsOnlyTrue (interfaceLongStringsOnly,"True","enabled.",true); static SwitchOption interfaceLongStringsOnlyFalse (interfaceLongStringsOnly,"False","disabled.",false); static Parameter interfaceStringm0 ("Stringm0", ".", &StringFragmentation::stringm0, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, true, false, Interface::lowerlim); static Parameter interfaceJunctionDiquark ("JunctionDiquark", "Suppress diquark production in breaking of junctions with this amount", &StringFragmentation::junctionDiquark, 0, 1.0, 0.0, 0, true, false, Interface::lowerlim); static Parameter interfaceAlpha ("Alpha", "Boost the enhanced string tension with additional factor", &StringFragmentation::alpha, 0, 1.0, 1.0, 0, true, false, Interface::lowerlim); static Switch interfaceAverage ("Average", "Disable to try and hadronise strings with individual tensions" "instead of average value." , &StringFragmentation::average, false, true, false); static SwitchOption interfaceAverageTrue (interfaceAverage,"True","enabled.",true); static SwitchOption interfaceAverageFalse (interfaceAverage,"False","disabled.",false); static Parameter interfaceStringyCut ("StringyCut", "No enhancement of strings with a constituent pT higher than StringpTCut (GeV), and within " " StringyCut.", &StringFragmentation::stringyCut, 0, 2.0, 0.0, 0, true, false, Interface::lowerlim); static Parameter interfaceStringR0 ("StringR0", "In the string overlap model the string R0 dictates the minimum radius " " of a string (in fm).", &StringFragmentation::stringR0, femtometer, 0.5*femtometer, 0.0*femtometer, 0*femtometer, true, false, Interface::lowerlim); static Parameter interfaceFragmentationMass ("FragmentationMass", "Set the mass used for the f(z) integral in the overlap string model " " default is pion mass (units GeV).", &StringFragmentation::fragmentationMass, 0, 0.135, 0., 1., true, false, Interface::lowerlim); static Parameter interfaceBaryonSuppression ("BaryonSuppression", "Set the fudge factor used for baryon suppression in the overlap string model. " " One day this should be calculated properly. This day is not today. Probably around 2...", &StringFragmentation::baryonSuppression, 0, 0.5, 0.0, 1.0, true, false, Interface::limited); static Parameter interfaceMaxTries ("MaxTries", "Sometimes Pythia gives up on an event too easily. We therefore allow " "it to re-try a couple of times.", &StringFragmentation::maxTries, 2, 1, 0, true, false, Interface::lowerlim); } string StringFragmentation::convert(double d) { ostringstream os; os << d; return os.str(); } diff --git a/Hadronization/StringFragmentation.h b/Hadronization/StringFragmentation.h --- a/Hadronization/StringFragmentation.h +++ b/Hadronization/StringFragmentation.h @@ -1,390 +1,390 @@ // -*- C++ -*- #ifndef THEP8I_StringFragmentation_H #define THEP8I_StringFragmentation_H // // This is the declaration of the StringFragmentation class. // #include "ThePEG/Handlers/HadronizationHandler.h" #include "ThePEG/Handlers/ClusterCollapser.h" #include "TheP8I/Config/Pythia8Interface.h" #include "OverlapPythiaHandler.h" #include "TheP8EventShapes.h" #include #include "histtools/HistKeeper.h" //#include "YODA/Histo1D.h" //#include "YODA/Histo2D.h" //#include "YODA/Writer.h" //#include "YODA/WriterYODA.h" namespace TheP8I { bool sorting_principle(pair > c1, pair > c2){ return (c1.first->momentum().perp() < c2.first->momentum().perp()); } using namespace ThePEG; /** * Here is the documentation of the StringFragmentation class. * * @see \ref StringFragmentationInterfaces "The interfaces" * defined for StringFragmentation. */ class StringFragmentation: public HadronizationHandler { typedef map PytPars; friend class Ropewalk::Dipole; public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ StringFragmentation(); /** * The destructor. */ virtual ~StringFragmentation(); //@} public: /** * The main function called by the EventHandler class to * perform the Hadronization step. * @param eh the EventHandler in charge of the Event generation. * @param tagged if not empty these are the only particles which should * be considered by the StepHandler. * @param hint a Hint object with possible information from previously * performed steps. * @throws Veto if the StepHandler requires the current step to be * discarded. * @throws Stop if the generation of the current Event should be stopped * after this call. * @throws Exception if something goes wrong. */ virtual void handle(EventHandler & eh, const tPVector & tagged, const Hint & hint); /** @name Functions used by the persistent I/O system. */ //@{ /** * Let the given Pythia8Interface hadronize the ColourSinglet * systems provided. */ bool hadronizeSystems(Pythia8Interface & pyt, const vector & singlets, const tPVector & all); void hadronizeTweak(Pythia8Interface & pyt, const ColourSinglet & singlet, const tPVector & all); /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Helper function to get the maximum transverse momenta of any * parton in a string. Optionally only consider partons within a * rapidity interval |eta| < deltaY. */ static Energy maxPT(const ColourSinglet & cs, double deltaY = 0.0); /** * 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); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ 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; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); //@} private: HistKeeper hk; /** * The interface to the Pythia 8 object, */ Pythia8Interface pythia; /** * Internal switch for the fragmentation scheme */ int fScheme; /** * Object that handles calculation of new Pythia parameters */ ParameterHandler phandler; /** * The intrinsic string radius */ Length stringR0; /** * The assumed mass for calculating rapidities in the reeperbahn scheme */ Energy stringm0; /** * Supression of diquarks emerging from breaking junctions */ double junctionDiquark; double alpha; bool average; /** *????????? If **window** is set, this determines the abs(y) window for which no enhancement will be made */ double stringyCut; /** * Cut on maxpT for a string to go into the ropewalk. REquires **longsoft** to be set */ Energy stringpTCut; /** * Assumed m0 value in calculation of normalization of f(z) the Lund splitting function. */ double fragmentationMass; /** * Parameter surpressing baryonic content further in calculation of effective parameters */ double baryonSuppression; // Force hadronize systems to run to the end, regardless of windows bool forcerun; /** * Can provide a window in y and pT where no enhancement is made, in order to not include * very hard gluons in central rapidity region in a rope */ bool window; /* * Minimal rapidity span for a string to go into the ropewalk. REquires **longsoft** to be set */ double minStringy; /* * Only allow long and soft strings to go into the ropewalk */ bool longSoft; bool longStringsOnly; /** * Throw away some string entirely while making ropes */ bool throwaway; /* * Let the strings shove each other around */ bool shoving; /* * The interval with which we slice in rapidity */ double deltay; /* * Total shoving -and propagation time */ Length tshove; /* * Time slicing */ Length deltat; /* * Cutoff radius at which we don't let dipoles affect each other any more * given as a multiple of R0 */ double rCutOff; - typedef Qty<-1,1,0> StringTension; + typedef Tension StringTension; /* * Amplitude of energy density function for string shoving */ StringTension gAmplitude; /* * Exponent multiple of energy density function for string shoving */ double gExponent; /* * Cutoff rapidity for merging of string excitations */ double yCutOff; /*Lorentz contraction of fields * */ bool lorentz; /* * Cylinder shoving */ bool cylinder; TheP8EventShapes * _eventShapes; /** * Use thrust axis to calculate rapidities; useful for LEP validation */ int useThrustAxis; /** * Additional interfaces to Pythia8 objects in case the string * overlap model is to be used. */ OverlapPythiaHandler * opHandler; /** * Sometimes Pythia gives up on an event too easily. We allow it to * re-try a couple of times. */ int maxTries; // Keep track of how many events we already hadronized using this object, in order to clean up int nev; // Can do a small analysis int doAnalysis; //vector _histograms; //vector _histograms2D; string analysisPath; void dofinish(); // Can print the strings of the event in a matlab friendly format string PrintStringsToMatlab(vector& singlets); string convert(double d); #include "StringFragmentation-var.h" /** * The object used to avoid too small strings in the hadronization. */ ClusterCollapserPtr theCollapser; public: /** * Exception class declaration. */ struct StringFragError: public Exception {}; private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ClassDescription initStringFragmentation; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ StringFragmentation & operator=(const StringFragmentation &); }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of StringFragmentation. */ template <> struct BaseClassTrait { /** Typedef of the first base class of StringFragmentation. */ typedef HadronizationHandler NthBase; }; /** This template specialization informs ThePEG about the name of * the StringFragmentation class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "TheP8I::StringFragmentation"; } /** * The name of a file containing the dynamic library where the class * StringFragmentation is implemented. It may also include several, space-separated, * libraries if the class StringFragmentation 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 "libTheP8I.so"; } }; /** @endcond */ } #endif /* THEP8I_StringFragmentation_H */ diff --git a/Hadronization/StringPipe.cc b/Hadronization/StringPipe.cc --- a/Hadronization/StringPipe.cc +++ b/Hadronization/StringPipe.cc @@ -1,154 +1,154 @@ #include "StringPipe.h" using namespace TheP8I; StringPipe::StringPipe() { } StringPipe::~StringPipe(){ } StringPipe::StringPipe(ColourSinglet* singlet, Length r0, TheP8EventShapes * es) : _externalOverlap(make_pair(0,0)), _es(es), theSinglet(singlet) { tcPPtr pmin = singlet->partons()[0]; tcPPtr pmax = singlet->partons()[0]; for (tcPVector::iterator pItr = singlet->partons().begin(); pItr!=singlet->partons().end();++pItr) { // Find the lowest and the highest rapidities to span the full string between double y = _es ? _es->yT((*pItr)->momentum()) : (*pItr)->rapidity(); if(y<(_es ? _es->yT(pmin->momentum()) : pmin->rapidity())) pmin = *pItr; else if(y>(_es ? _es->yT(pmax->momentum()) : pmax->rapidity())) pmax = *pItr; } _originb.first = (pmin->vertex().x() + pmax->vertex().x()) / 2; _originb.second = (pmin->vertex().y() + pmax->vertex().y()) / 2; _ymin = _es ? _es->yT(pmin->momentum()) : pmin->rapidity(); _ymax = _es ? _es->yT(pmax->momentum()) : pmax->rapidity(); pair Vcs = make_pair(0*femtometer*femtometer,0*femtometer*femtometer); tcPVector::iterator pminus = singlet->partons().begin(); // Calculate pipe radius using distance from parton to mean line // between partons with max and min probability for (tcPVector::iterator pItr = singlet->partons().begin(); pItr!=singlet->partons().end();++pItr) { // Grow the pipe radius Length dx = (*pItr)->vertex().x() - _originb.first; Length dy = (*pItr)->vertex().y() - _originb.second; _radius2 += dx*dx + dy*dy; // Grow the volume of the colour singlet if(pItr!=singlet->partons().begin()){ // Length dbx = (*pItr)->vertex().x() - (*pminus)->vertex().x(); // Length dby = (*pItr)->vertex().y() - (*pminus)->vertex().y(); double dy = _es ? _es->yT((*pItr)->momentum()) - _es->yT((*pminus)->momentum()) : (*pItr)->rapidity() - (*pminus)->rapidity(); // ORIGINAL BELOW -- THIS IS CHANGED PER 20.05.2014, pipe ends should now be parallel // to impact parameter plane. Below they are allowed not to be. if(dy>0) Vcs.first += r0 * r0 * abs(dy); if(dy<0) Vcs.second += r0 * r0 * abs(dy); // Area db2 = dbx * dbx + dby * dby; // if(dy>0) Vcs.first += (r0 * r0 * dy * dy) / (sqrt(db2/(r0*r0) +dy*dy)); // if(dy<0) Vcs.second += (r0 * r0 * dy * dy) / (sqrt(db2/(r0*r0) +dy*dy)); ++pminus; } } // Use radius squared instead of volume to save us from a factor of pi. _radius2 /= (singlet->partons().end() - singlet->partons().begin()); _radius2 += r0*r0; // Overlap is given in both directions _internalOverlap.first = Vcs.first / (_radius2 * (_ymax - _ymin)); _internalOverlap.second = Vcs.second / (_radius2 * (_ymax - _ymin)); } pair StringPipe::ExternalOverlap(StringPipe& other){ if(this==&other || other.GetVolume() == 0*femtometer*femtometer){ return make_pair(0,0); } return make_pair(other._internalOverlap.first * OverlapY(other) * OverlapArea(other) / other.GetVolume(), other._internalOverlap.second * OverlapY(other) * OverlapArea(other) / other.GetVolume()); } Area StringPipe::GetVolume(){ return _radius2*M_PI*(_ymax - _ymin); } Area StringPipe::OverlapArea(StringPipe& other){ // This is the overlap between two cylinders. Purely geometrical. Area d2 = (_originb.first - other._originb.first) * (_originb.first - other._originb.first) + (_originb.second - other._originb.second) * (_originb.second - other._originb.second); Length d = sqrt(d2); Length r = sqrt(_radius2); Length R = sqrt(other._radius2); - Qty<4,0,0,1,1,1> argu = (-d + r + R) * (d + r - R) * (d - r + R) * (d + r + R); + auto argu = (-d + r + R) * (d + r - R) * (d - r + R) * (d + r + R); // Save us from dividing by zero if(argu <= 0*femtometer*femtometer*femtometer*femtometer) return 0*femtometer*femtometer; return _radius2 * acos((d2 + _radius2 - other._radius2)/(2*d*r)) + other._radius2 * acos((d2 + other._radius2 - _radius2)/(2*d*R)) - 0.5 * sqrt(argu); } double StringPipe::IntevalOverlap(double min1, double min2, double max1, double max2){ // Helper function to calculate overlap between two 1D intervals. // This could very well exist somewhere in ThePEG already. if(max1partons().begin(); pItr!=theSinglet->partons().end();++pItr) { if(pT < (*pItr)->momentum().perp()){ pT = (*pItr)->momentum().perp(); y = (*pItr)->rapidity(); } } return y; } Energy StringPipe::MeanpT(){ Energy pT = 0*GeV; for (tcPVector::iterator pItr = theSinglet->partons().begin(); pItr!=theSinglet->partons().end();++pItr) { pT += (*pItr)->momentum().perp(); } return pT/(distance(theSinglet->partons().begin(),theSinglet->partons().end())); } Energy StringPipe::MaxpT(){ Energy pT = 0*GeV; for (tcPVector::iterator pItr = theSinglet->partons().begin(); pItr!=theSinglet->partons().end();++pItr) { if(pT < (*pItr)->momentum().perp()) pT = (*pItr)->momentum().perp(); } return pT; } double StringPipe::OverlapY(StringPipe& other){ return IntevalOverlap(_ymin, other._ymin, _ymax, other._ymax); } ColourSinglet * StringPipe::GetSingletPtr(){ return theSinglet; -} \ No newline at end of file +} diff --git a/src/SimpleLEP.in b/src/SimpleLEP.in --- a/src/SimpleLEP.in +++ b/src/SimpleLEP.in @@ -1,15 +1,24 @@ cd /Defaults/Generators set SimpleLEPGenerator:NumberOfEvents 10000 set SimpleLEPGenerator:DebugLevel 1 set SimpleLEPGenerator:PrintEvent 10 set SimpleLEPGenerator:EventHandler:LuminosityFunction:Energy 91.2 set SimpleLEPGenerator:EventHandler:DecayHandler NULL set SimpleLEPGenerator:DumpPeriod 7000 set SimpleLEPGenerator:EventHandler:HadronizationHandler /TheP8I/Handlers/StringFrag/StringFragmenter create TheP8I::BoseEinsteinHandler BEHandler set SimpleLEPGenerator:EventHandler:DecayHandler BEHandler set /TheP8I/Particles/e+:PDF /Defaults/Partons/SimpleLeptonPDF set /TheP8I/Particles/e-:PDF /Defaults/Partons/SimpleLeptonPDF set /TheP8I/Particles/pi0:Stable Stable set SimpleLEPGenerator:Strategy /TheP8I/StdStrategy +create ThePEG::HepMCFile HepMC2 HepMCAnalysis.so +cp HepMC2 HepMC3 +set HepMC2:Format GenEvent +set HepMC3:Format GenEventHepMC3 +set HepMC2:Filename SimpleLEP.hepmc2 +set HepMC3:Filename SimpleLEP.hepmc3 + +insert SimpleLEPGenerator:AnalysisHandlers[0] HepMC2 +insert SimpleLEPGenerator:AnalysisHandlers[0] HepMC3 saverun SimpleLEP SimpleLEPGenerator