diff --git a/Hadronization/Ropewalk.cc b/Hadronization/Ropewalk.cc --- a/Hadronization/Ropewalk.cc +++ b/Hadronization/Ropewalk.cc @@ -1,1261 +1,1262 @@ // -*- 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)); diff = sel[UseRandom::rnd()]; ++ns; } else { double mac = (pa->momentum() + od->pc->momentum()).m()/GeV; double mbd = (pc->momentum() + od->pa->momentum()).m()/GeV; 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/StringFragmentation.cc b/Hadronization/StringFragmentation.cc --- a/Hadronization/StringFragmentation.cc +++ b/Hadronization/StringFragmentation.cc @@ -1,1100 +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); + 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); 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(); }