Page MenuHomeHEPForge
Files F10881771

No OneTemporary
Actions

View Options

diff --git a/DIPSY/DipoleState.cc b/DIPSY/DipoleState.cc
--- a/DIPSY/DipoleState.cc
+++ b/DIPSY/DipoleState.cc
@@ -1,1461 +1,1461 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DipoleState class.
//
#include "DipoleState.h"
#include "DipoleEventHandler.h"
#include "EventFiller.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/Current.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/MaxCmp.h"
#include "ThePEG/Utilities/ObjectIndexer.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ParticleInfo.h"
#include "CPUTimer.h"
#include <stdio.h>
#include <iostream>
#include <fstream>
using namespace DIPSY;
DipoleState::DipoleState(const DipoleState & x)
: thePlus(x.thePlus), theMinus(x.theMinus), theMinusDeficit(x.theMinusDeficit),
theHandler(x.theHandler), theInitialDipoles(x.theInitialDipoles),
theSwingCandidates(x.theSwingCandidates),
theTouchedSwingCandidates(x.theTouchedSwingCandidates),
theWFInfo(x.theWFInfo),
theWeight(x.theWeight), doTakeHistory(x.doTakeHistory), theYmax(x.theYmax),
theCollidingEnergy(x.theCollidingEnergy), theHistory(x.theHistory),
allDipoles(x.allDipoles), theIncoming(x.theIncoming),
theProducedParticles(x.theProducedParticles), theIncomingShadows(x.theIncomingShadows),
thePrimaryShadow(x.thePrimaryShadow) {
for ( set<DipolePtr>::iterator it = allDipoles.begin();
it != allDipoles.end(); ++it )
(**it).theDipoleState = this;
}
DipoleState::~DipoleState() {}
DipoleStatePtr DipoleState::clone() {
DipoleStatePtr copy = new_ptr(*this);
for ( set<DipolePtr>::iterator it = allDipoles.begin();
it != allDipoles.end(); ++it )
(**it).theDipoleState = this;
Dipole::CloneSet toclone;
for ( set<DipolePtr>::iterator it = allDipoles.begin();
it != allDipoles.end(); ++it ) {
toclone.insert(*it);
(**it).fillReferences(toclone);
}
Dipole::TranslationMap trans;
vector<Ariadne5::ClonePtr> clones;
for ( Dipole::CloneSet::iterator it = toclone.begin();
it != toclone.end(); ++it )
if ( *it ) clones.push_back(trans[*it] = (**it).clone());
for ( int i = 0, N = clones.size(); i < N; ++i ) {
clones[i]->rebind(trans);
if ( tDipolePtr d = dynamic_ptr_cast<tDipolePtr>(clones[i]) ) {
d->theDipoleState = copy;
}
}
copy->allDipoles.clear();
trans.translate(inserter(copy->allDipoles),
allDipoles.begin(), allDipoles.end());
copy->theInitialDipoles.clear();
trans.translate(inserter(copy->theInitialDipoles),
theInitialDipoles.begin(), theInitialDipoles.end());
for ( int i = 0, N = theSwingCandidates.size(); i < N; ++i ) {
copy->theSwingCandidates =
vector< vector<tDipolePtr> >(theSwingCandidates.size());
trans.translate(inserter(copy->theSwingCandidates[i]),
theSwingCandidates[i].begin(), theSwingCandidates[i].end());
}
for ( int i = 0, N = theTouchedSwingCandidates.size(); i < N; ++i ) {
copy->theTouchedSwingCandidates =
vector< vector<tDipolePtr> >(theTouchedSwingCandidates.size());
trans.translate(inserter(copy->theTouchedSwingCandidates[i]),
theTouchedSwingCandidates[i].begin(),
theTouchedSwingCandidates[i].end());
}
return copy;
// *** TODO *** fix shadows here
}
void DipoleState::sortDipolesFS() {
theSwingCandidates.clear();
theTouchedSwingCandidates.clear();
theSwingCandidates.resize(Current<DipoleEventHandler>()->nColours());
theTouchedSwingCandidates.resize(Current<DipoleEventHandler>()->nColours());
list<DipolePtr> dips = getDipoles();
for ( list<DipolePtr>::const_iterator it = dips.begin(); it != dips.end(); it++ ) {
sortDipoleFS(**it);
}
}
void dummycollect(const Dipole & d, int & ndip, Energy & sump, Energy & summ,
TransverseMomentum & sumpt, set<tcDipolePtr> & dips) {
if ( d.children().first && d.children().second ) {
dummycollect(*d.children().first, ndip, sump, summ, sumpt, dips);
dummycollect(*d.children().second, ndip, sump, summ, sumpt, dips);
} else if ( d.children().first ) {
dummycollect(*d.children().first, ndip, sump, summ, sumpt, dips);
} else if ( d.children().second ) { //dead end??
// dummycollect(*d.children().second, ndip, sump, summ, sumpt, dips);
} else {
++ndip;
dips.insert(&d);
if ( d.partons().first->flavour() == ParticleID::g ) {
sumpt += d.partons().first->pT()/2.0;
sump += d.partons().first->plus()/2.0;
summ += d.partons().first->minus()/2.0;
} else {
sumpt += d.partons().first->pT();
sump += d.partons().first->plus();
summ += d.partons().first->minus();
}
if ( d.partons().second->flavour() == ParticleID::g ) {
sumpt += d.partons().second->pT()/2.0;
sump += d.partons().second->plus()/2.0;
summ += d.partons().second->minus()/2.0;
} else {
sumpt += d.partons().second->pT();
sump += d.partons().second->plus();
summ += d.partons().second->minus();
}
}
}
pair<double,int> DipoleState::lambdaMeasure(Energy2 scale,
FactoryBase::tH1DPtr histlength,
FactoryBase::tH1DPtr histmass) const {
pair<double,int> lam (0.0, 0);
for ( set<DipolePtr>::const_iterator it = allDipoles.begin();
it != allDipoles.end(); it++ ) {
Dipole & d = **it;
if ( !d.children().first && !d.children().second && d.participating() ) {
Energy2 m2 = (d.partons().first->momentum() +
d.partons().second->momentum()).m2();
if ( m2 > ZERO ) lam.first += log(m2/scale);
++lam.second;
if ( histmass && m2 >= ZERO ) histmass->fill(sqrt(m2/GeV2));
if ( histlength )
histlength->fill((d.partons().first->position() -
d.partons().second->position()).pt()*GeV);
}
}
return lam;
}
void DipoleState::sortFSDipoles() {
theSwingCandidates.clear();
theSwingCandidates.resize(Current<DipoleEventHandler>()->nColours());
for ( set<DipolePtr>::const_iterator it = allDipoles.begin();
it != allDipoles.end(); it++ ) {
Dipole & d = **it;
if ( !d.children().first && !d.children().second && d.participating() ) {
if ( ( d.neighbors().first &&
d.neighbors().first->colour() == d.colour() ) ||
( d.neighbors().second &&
d.neighbors().second->colour() == d.colour() ) )
generateColourIndex(&d);
forceAt(theSwingCandidates, (**it).colour()).push_back(*it);
}
}
if ( theSwingCandidates.size() != 0 ) return;
set<tcDipolePtr> dipset1;
TransverseMomentum sumpt;
Energy sump = ZERO;
Energy summ = ZERO;
int ndip = 0;
for ( int ic = 0, Nc = theSwingCandidates.size(); ic < Nc; ++ ic )
for ( int i = 0, N = theSwingCandidates[ic].size(); i < N; ++i ) {
const Dipole & d = *theSwingCandidates[ic][i];
dipset1.insert(&d);
++ndip;
if ( d.partons().first->flavour() == ParticleID::g ) {
sumpt += d.partons().first->pT()/2.0;
sump += d.partons().first->plus()/2.0;
summ += d.partons().first->minus()/2.0;
} else {
sumpt += d.partons().first->pT();
sump += d.partons().first->plus();
summ += d.partons().first->minus();
}
if ( d.partons().second->flavour() == ParticleID::g ) {
sumpt += d.partons().second->pT()/2.0;
sump += d.partons().second->plus()/2.0;
summ += d.partons().second->minus()/2.0;
} else {
sumpt += d.partons().second->pT();
sump += d.partons().second->plus();
summ += d.partons().second->minus();
}
}
cerr << "<DIPSY> total momentum in sortFSDipoles." << endl
<< " Plus: " << sump/GeV << endl
<< " Minus: " << summ/GeV << endl
<< " x: " << sumpt.x()/GeV << endl
<< " y: " << sumpt.y()/GeV << endl
<< " Ndip: " << ndip << " (" << dipset1.size() << ")" << endl;
sump = ZERO;
summ = ZERO;
sumpt = TransverseMomentum();
ndip = 0;
set<tcDipolePtr> dipset2;
for ( int i = 0, N = initialDipoles().size(); i < N; ++i )
dummycollect(*initialDipoles()[i], ndip, sump, summ, sumpt, dipset2);
cerr << "<DIPSY> total momentum in getFSSwinger." << endl
<< " Plus: " << sump/GeV << endl
<< " Minus: " << summ/GeV << endl
<< " x: " << sumpt.x()/GeV << endl
<< " y: " << sumpt.y()/GeV << endl
<< " Ndip: " << ndip << " (" << dipset2.size() << ")" << endl;
set<tcDipolePtr> diffset;
set_difference(dipset1.begin(), dipset1.end(), dipset2.begin(), dipset2.end(),
inserter(diffset));
cerr << " " << diffset.size()
<< " dipoles are in sortFSDipoles but not in getFSSwinger" << endl;
diffset.clear();
set_difference(dipset2.begin(), dipset2.end(), dipset1.begin(), dipset1.end(),
inserter(diffset));
cerr << " " << diffset.size()
<< " dipoles are in getFSSwinger but not in sortFSDipoles" << endl;
}
void DipoleState::sortDipoleFS(Dipole & d){
if(d.children().second && d.children().first) {
sortDipole(*(d.children().first));
sortDipole(*(d.children().second)); }
else if(d.children().first)
sortDipole(*(d.children().first));
else if ( d.children().second );
else if ( !(d.participating()) );
else forceAt(theSwingCandidates, d.colour()).push_back(&d);
}
void DipoleState::sortDipoles() {
theSwingCandidates.clear();
theSwingCandidates.resize(Current<DipoleEventHandler>()->nColours());
theTouchedSwingCandidates.clear();
theTouchedSwingCandidates.resize(Current<DipoleEventHandler>()->nColours());
for ( int i = 0, N = theInitialDipoles.size(); i < N; ++i ) {
sortDipole(*(theInitialDipoles[i]));
}
}
void DipoleState::sortDipole(Dipole & d){
if(d.children().second && d.children().first) {
sortDipole(*(d.children().first));
sortDipole(*(d.children().second)); }
else if(d.children().first)
sortDipole(*(d.children().first));
else if ( d.children().second );
else if ( d.interacted() );
else if ( !(d.participating()) );
else if ( !(d.isOn) );
else {
forceAt(theSwingCandidates, d.colour()).push_back(& d);
if ( !d.hasGen() )
forceAt(theTouchedSwingCandidates, d.colour()).push_back(& d);
}
}
void DipoleState::evolve(double ymin, double ymax) {
if ( Current<DipoleEventHandler>()->effectivePartonMode() < 0 && !hasShadows() )
setupShadows();
save();
for ( set<DipolePtr>::const_iterator it = allDipoles.begin(); it != allDipoles.end(); it++ ) {
(*it)->reset();
(*it)->touch();
}
while ( true ) {
sortDipoles();
tDipolePtr sel = tDipolePtr();
for ( int i = 0, N = initialDipoles().size(); i < N; ++i ) {
tDipolePtr si = initialDipoles()[i]->getEmitter(ymin, ymax);
if ( si && si->generatedY() < ymax && si->hasGen() &&
( !sel || si->generatedY() < sel->generatedY() ))
sel = si;
}
if ( sel ) {
ymin = sel->generatedY();
sel->emit();
save();
continue;
}
break;
}
theYmax = ymax;
}
void DipoleState::swingFS(double ymin, double ymax) {
if ( !handler().swingPtr() ) return;
const Swinger & swinger = handler().swinger();
for ( set<DipolePtr>::const_iterator it = allDipoles.begin();
it != allDipoles.end(); it++ ) {
(*it)->reset();
(*it)->touch();
}
sortFSDipoles();
double miny = ymin;
int lastcol = -1;
while ( true ) {
MinCmp<double, Dipole *> sel;
for ( int ic = 0, Nc = theSwingCandidates.size(); ic < Nc; ++ic ) {
const vector<tDipolePtr> & candidates = swingCandidates(ic);
if ( lastcol < 0 || lastcol == ic )
for ( int i = 0, N = candidates.size(); i < N; ++i ) {
Dipole & d = *candidates[i];
InvEnergy2 saver = d.swingCache;
if ( saver < ZERO )
saver = swinger.swingDistanceFS(*d.partons().first,
*d.partons().second, ymin/GeV2);
if ( d.swingDipole() && d.swingDipole()->touched() ) d.reset();
d.swingCache = saver;
}
for ( int i = 0, N = candidates.size(); i < N - 1; ++i ) {
Dipole & d1 = *candidates[i];
bool safe = d1.hasGen();
if ( lastcol < 0 || lastcol == ic )
for ( int j = i + 1; j < N; ++j ) {
Dipole & d2 = *candidates[j];
if ( safe && !d2.touched() ) continue;
InvEnergy2 c =
swinger.swingDistanceFS(*d1.partons().first,
*d2.partons().second, ymin/GeV2);
InvEnergy2 d =
swinger.swingDistanceFS(*d2.partons().first,
*d1.partons().second, ymin/GeV2);
double R = -log( UseRandom::rnd() );
double amp = swinger.swingAmpFS(d1.swingCache, d2.swingCache, c, d);
double yi = Constants::MaxRapidity;
if ( miny*amp + R < Constants::MaxRapidity*amp ) yi = miny + R/amp;
if ( yi < d1.generatedY() || !d1.hasGen() ) {
d1.swingDipole(&d2);
d1.generatedY(yi);
d1.recoilSwing(false);
}
}
sel(d1.generatedY(), &d1);
d1.untouch();
}
if ( candidates.size() ) candidates.back()->untouch();
}
if ( !sel.index() || sel > ymax ) break;
miny = sel;
lastcol = sel.index()->colour();
swinger.recombineFS(*sel.index());
static DebugItem tupleswing("DIPSY::SwingTuple", 6);
static ofstream swingtuple;
static bool isopen = false;
static long lastevent = 0;
if ( tupleswing ) {
if ( !isopen ) {
string filename = CurrentGenerator::current().filename() + "-swing.tuple";
swingtuple.open(filename.c_str());
isopen = true;
}
if ( lastevent != CurrentGenerator::current().currentEventNumber() ) {
lastevent = CurrentGenerator::current().currentEventNumber();
swingtuple << "# " << lastevent << endl;
}
double lrho = -log10(sqrt(miny));
pair<double,int> lam = lambdaMeasure(0.36*GeV2);
swingtuple << lrho << '\t' << lam.first/lam.second << '\t' << lam.first << endl;
}
}
theYmax = ymax;
}
void DipoleState::restoreNonparticipants() {
for ( int i = 0; i < int(theInitialDipoles.size()); i++ ) {
tPartonPtr p1 = theInitialDipoles[i]->partons().first;
tPartonPtr p2 = theInitialDipoles[i]->partons().second;
tPartonPtr p3;
if ( i != int(theInitialDipoles.size())-1 &&
theInitialDipoles[i+1]->partons().first == p2 )
p3 = theInitialDipoles[i+1]->partons().second;
else if ( i != int(theInitialDipoles.size())-1 &&
theInitialDipoles[i+1]->partons().second == p1 )
p3 = theInitialDipoles[i+1]->partons().first;
else if ( i != 0 && theInitialDipoles[i-1]->partons().first == p2 )
p3 = theInitialDipoles[i-1]->partons().second;
else if ( i != 0 && theInitialDipoles[i-1]->partons().second == p1 )
p3 = theInitialDipoles[i-1]->partons().first;
if ( !(p1->interacted()) &&
!(p2->interacted()) &&
!( p3 && p3->interacted() ) ) {
theInitialDipoles[i]->participating(false);
restoreDipole(p1, p2);
}
else
theInitialDipoles[i]->participating(true);
}
}
bool DipoleState::restoreDipole(PartonPtr p1, PartonPtr p2) {
//there must be dipoles to reconnect.
if ( !p1->dipoles().second || !p2->dipoles().first ) return false;
//if already connect, do nothing
if ( p1->dipoles().second == p2->dipoles().first ) return true;
//if p1 and p2 are connected to different partons, just swing.
if ( p1->dipoles().second->partons().second != p2->dipoles().first->partons().first ) {
p1->dipoles().second->swingDipole(p2->dipoles().first);
p1->dipoles().second->recombine();
return true;
}
//if p1 and p2 are separated by only one parton, then first swing away one of
//the dipoles, then repeat.
else {
sortDipoles();
//if no swing can be found, then leave it.
if ( !forceSwing(p1->dipoles().second->partons().second, 0.0, 0.0) ) return false;
//otherwise go ahead and reconnect the dipoles.
p1->dipoles().second->swingDipole(p2->dipoles().first);
p1->dipoles().second->recombine();
return true;
}
}
void DipoleState::normaliseValenceCharge(int mode) {
if ( mode == 0 ) return;
if ( double(int(theInitialDipoles.size())/3) != double(theInitialDipoles.size())/3.0 ) {
return;
}
if ( mode == 1 ) {
set<pair<tPartonPtr, tPartonPtr> > valencePartons;
for ( int i = 0; i < int(theInitialDipoles.size()); i++) {
tPartonPtr p1 = theInitialDipoles[i]->partons().first;
tPartonPtr p2 = theInitialDipoles[i]->partons().second;
//if there are no dipoles to swing, do nothing
if ( !p1->dipoles().second || !p2->dipoles().first )
continue;
//if already connected, do nothing.
if ( p1->dipoles().second->partons().second == p2 )
continue;
//if connected to the same parton, do nothing.
if ( p1->dipoles().second->partons().second == p2->dipoles().first->partons().first )
continue;
//otherwise swing with 50% prob
if ( UseRandom::rnd() > 0.5 ) {
p1->dipoles().second->swingDipole(p2->dipoles().first);
p1->dipoles().second->recombine();
}
}
}
}
void DipoleState::generateColourIndex(tDipolePtr d) {
int isys = -1;
if ( d->neighbors().first )
isys = d->neighbors().first->colourSystem();
if ( d->neighbors().second ) {
int isys2 = d->neighbors().second->colourSystem();
if ( isys >= 0 && isys2 != isys )
cerr << "Dipole bewteen different colour systems!" << endl;
isys = isys2;
}
do {
d->colour(UseRandom::irnd(handler().nColours()), isys);
} while ( ( d->neighbors().first &&
d->colour() == d->neighbors().first->colour() ) ||
( d->neighbors().second &&
d->colour() == d->neighbors().second->colour() ) );
}
const list<PartonPtr> DipoleState::virtualPartons() const {
list<PartonPtr> partons;
for ( set<DipolePtr>::const_iterator it = allDipoles.begin() ; it != allDipoles.end(); it++ ) {
if ( (*it)->children().first || (*it)->children().second ) continue;
if ( !((*it)->neighbors().first) && !((*it)->partons().first->interacted()) )
if ( !((*it)->partons().first->valence()) )
partons.push_back((*it)->partons().first);
if ( !((*it)->partons().second->interacted()) && !((*it)->partons().second->valence()) )
partons.push_back((*it)->partons().second);
}
return partons;
}
const set< list<PartonPtr> > DipoleState::loops() const {
list<PartonPtr> partonsList = getPartons();
set<PartonPtr> partons;
set< list<PartonPtr> > ret;
for( list<PartonPtr>::const_iterator it = partonsList.begin();
it != partonsList.end(); it++ ) {
partons.insert(*it);
}
//go through all partons and add them to colour chains.
while ( !partons.empty() ) {
list<PartonPtr> loop;
PartonPtr head = *partons.begin();
bool forward = true;
PartonPtr p = head;
//follow the colour flow first in forward direction from head until it comes back to the start
//if a quark is hit, then continue in other direction.
do {
//add at end/start of list depending of if we are moving forward or beckwards in colour
if ( forward ) {
loop.push_back(p);
}
else {
loop.push_front(p);
}
//remove from partons, so we know its added.
//could potentially be better to do this without destroy the dipoleState...
if ( !partons.erase(p) ) {
Throw<DipoleConnectionException>()
<< "DIPSY found inconsistent DipoleState when extracting "
"colour singlets!" << Exception::eventerror;
}
//move forward in colour if there are more partons, otherwise turn around
if ( forward ) {
if ( p->dipoles().second ) {
p = p->dipoles().second->partons().second;
if ( Debug::level > 5 ) cout << " continue forward to " << p->y() << endl;
}
else {
if ( head->dipoles().first ) {
//continue at other side, going in other direction
p = head->dipoles().first->partons().first;
}
else {
//want to turn around, but "head" is at the other end, so we are done.
p = head;
}
forward = false;
}
}
else {
if ( p->dipoles().first ) {
p = p->dipoles().first->partons().first;
}
else {
break;
}
}
} while ( p!= head );
ret.insert(loop);
}
return ret;
}
const double DipoleState::highestY() const {
MaxCmp<double> ret;
list< PartonPtr > partons = getPartons();
for ( list<PartonPtr>::const_iterator it = partons.begin();
it != partons.end(); it++ ) {
ret((*it)->y());
ret((*it)->y());
}
return ret;
}
const double DipoleState::lowestY() const {
MinCmp<double> ret;
list< PartonPtr > partons = getPartons();
for ( list<PartonPtr>::const_iterator it = partons.begin();
it != partons.end(); it++ ) {
ret((*it)->y());
ret((*it)->y());
}
return ret;
}
const double DipoleState::ymax() const {
return theYmax;
}
const Energy DipoleState::collidingEnergy() const {
return theCollidingEnergy;
}
void DipoleState::collidingEnergy(Energy E) {
theCollidingEnergy = E;
}
bool DipoleState::diagnosis(bool print) const {
bool ok = true;
int activeDipoles = 0;
int partons = 0;
int interactingDips = 0;
Energy plus = 0.0*GeV;
Energy minus = 0.0*GeV;
Energy plusNeeded = 0.0*GeV;
Energy minusNeeded = 0.0*GeV;
Energy plusMissing = 0.0*GeV;
Energy minusMissing = 0.0*GeV;
bool hasQuark = false;
TransverseMomentum totalpT;
Energy maxpT = 0.0*GeV;
InvEnergy dipSize = ZERO;
for(set< DipolePtr >::iterator it = allDipoles.begin();
it != allDipoles.end(); it++) {
if( ((*it)->children().first || (*it)->children().second) )
continue;
activeDipoles++;
dipSize += (*it)->size();
if ( (*it)->DGLAPsafe() ) {
if ( (*it)->neighbors().first &&
(*it)->neighbors().first->size() > (*it)->size() &&
!((*it)->neighbors().first->DGLAPsafe()) )
Throw<DipoleDGLAPSafeException>()
<< "DIPSY found broken DGLAPchain!" << Exception::eventerror;
if ( (*it)->neighbors().second &&
(*it)->neighbors().second->size() > (*it)->size() &&
!((*it)->neighbors().second->DGLAPsafe()) )
Throw<DipoleDGLAPSafeException>()
<< "DIPSY found broken DGLAPchain!" << Exception::eventerror;
}
if ( (*it)->interacted() ) interactingDips++;
if ( ((*it)->neighbors().first &&
(*it)->neighbors().first->children().first) ||
((*it)->neighbors().first &&
(*it)->neighbors().first->children().second) )
Throw<DipoleConnectionException>()
<< "DIPSY found dipole where negiboring dipole has decayed!"
<< Exception::eventerror;
if( (*it)->partons().first->dipoles().second != (*it) )
Throw<DipoleConnectionException>()
<< "DIPSY found dipole where first parton "
<< "doesn't point back at the dipole!" << Exception::eventerror;
if( (*it)->partons().second->dipoles().first != (*it) )
Throw<DipoleConnectionException>()
<< "DIPSY found dipole where second parton "
<< "doesn't point back at the dipole!" << Exception::eventerror;
if ( (*it)->neighbors().first &&
(*it)->neighbors().first->neighbors().second != (*it) )
Throw<DipoleConnectionException>()
<< "DIPSY found dipole where neighbouring dipole "
<< "doesn't point back at the dipole!" << Exception::eventerror;
if ( (*it)->neighbors().second &&
(*it)->neighbors().second->neighbors().first != (*it) )
Throw<DipoleConnectionException>()
<< "DIPSY found dipole where neighbouring dipole "
<< "doesn't point back at the dipole!" << Exception::eventerror;
if ( (*it)->partons().first == (*it)->partons().second )
Throw<DipoleConnectionException>()
<< "DIPSY found dipole looping back to itself!" << Exception::eventerror;
PartonPtr p = (*it)->partons().first;
partons++;
if ( p->dipoles().first && p->dipoles().first->partons().second != p )
Throw<DipoleConnectionException>()
<< "DIPSY found parton where first dipole "
<< "doesn't point back at the parton!" << Exception::eventerror;
if ( p->dipoles().second && p->dipoles().second->partons().first != p )
Throw<DipoleConnectionException>()
<< "DIPSY found parton where second dipole "
<< "doesn't point back at the parton!" << Exception::eventerror;
if ( p->dipoles().first && p->dipoles().second &&
p->dipoles().first == p->dipoles().second )
Throw<DipoleConnectionException>()
<< "DIPSY found parton colour connected to itself!"
<< Exception::eventerror;
if ( p->dipoles().first && (p->dipoles().first->children().first ||
p->dipoles().first->children().second) )
Throw<DipoleConnectionException>()
<< "DIPSY found parton where first dipole has decayed!"
<< Exception::eventerror;
if ( p->dipoles().second && (p->dipoles().second->children().first ||
p->dipoles().second->children().second) )
Throw<DipoleConnectionException>()
<< "DIPSY found parton where second dipole has decayed!"
<< Exception::eventerror;
if ( p->y() < -100.0 || p->y() > 100.0 )
Throw<DipoleKinematicsException>()
<< "DIPSY found parton impossible rapidity!" << Exception::eventerror;
if ( p->plus() < ZERO )
Throw<DipoleKinematicsException>()
<< "DIPSY found parton at oY " << p->oY() << " has negative plus "
<< p->plus()/GeV << "!" << Exception::eventerror;
if ( p->minus() < ZERO )
Throw<DipoleKinematicsException>()
<< "DIPSY found parton at oY " << p->oY() << " has negative minus "
<< p->minus()/GeV << "!" << Exception::eventerror;
totalpT += p->pT();
if ( p->pT().pt() > maxpT ) {
maxpT = p->pT().pt();
if ( maxpT > 1000.0*GeV )
Throw<DipoleKinematicsException>()
<< "DIPSY found parton with very large transverse momentum: "
<< maxpT/GeV << " (p+ = " << p->plus()/GeV << ", p- = "
<< p->minus()/GeV << ", y = " << p->y() << ")" << Exception::warning;
}
if( abs((p->plus()*p->minus() - sqr(p->pT().pt()))/sqr(GeV)) >
0.001*sqr(p->pT().pt())/sqr(GeV) ) {
Throw<DipoleKinematicsException>()
<< "DIPSY found off-shell parton with invariant mass square "
<< (p->plus()*p->minus() - sqr(p->pT().pt()))/sqr(GeV)
<< "!" << Exception::eventerror;
ok = false;
}
if( p->rightMoving() ) {
plus += p->plus();
minusMissing += p->minus();
minusNeeded += p->pT().pt()*exp( p->y() );
}
else {
minus += p->minus();
plusMissing += p->plus();
plusNeeded += p->pT().pt()*exp( -p->y() );
}
if( !((*it)->neighbors().second) ) {
hasQuark = true;
p = (*it)->partons().second;
partons++;
totalpT += p->pT();
if( p->rightMoving() ) {
plus += p->plus();
minusMissing += p->minus();
minusNeeded += p->pT().pt()*exp( p->y() );
}
else {
minus += p->minus();
plusMissing += p->plus();
plusNeeded += p->pT().pt()*exp( -p->y() );
}
}
}
if( totalpT.pt() > 0.00000001*GeV ) {
Throw<DipoleKinematicsException>()
<< "DIPSY found transverse momentum imbalance: " << totalpT.pt()/GeV
<< Exception::warning;
ok = false;
}
Energy acceptableError = 0.1*GeV;
if ( handler().eventFiller().mode() == 1 )
acceptableError = 0.000001*GeV;
if ( abs(plus - thePlus) > 2.0*acceptableError &&
abs(plus + minus + minusMissing + plusMissing - thePlus - theMinus) >
2.0*acceptableError ) {
Throw<DipoleKinematicsException>()
<< "DIPSY found energy non-conservation: "
<< (plus + minus + minusMissing + plusMissing
- thePlus - theMinus)/GeV/2.0 << " GeV" << Exception::warning;
ok = false;
}
if ( isnan(plus/GeV) || isnan(minus/GeV) )
Throw<DipoleKinematicsException>()
<< "DIPSY found parton with nan kinematics!" << Exception::runerror;
if ( print ) {
cout << "------------------ state " << this << " ------------------------------\n";
if (!ok )
cout << "| NOT OK!!!" << endl;
cout << setprecision(10);
cout <<
"| original plus = " << thePlus/GeV << endl <<
"| originl minus = " << theMinus/GeV << endl <<
"| plus = " << plus/GeV << endl <<
"| minus = " << minus/GeV << endl <<
"| missing minus = " << minusMissing/GeV << endl <<
"| needed minus = " << minusNeeded/GeV << endl <<
"| missing plus = " << plusMissing/GeV << endl <<
"| needed plus = " << plusNeeded/GeV << endl <<
"| total plus = " << (plus + plusNeeded - plusMissing)/GeV << endl <<
"| total minus = " << (minus + minusNeeded - minusMissing)/GeV << endl <<
"| total pT = " << "( "<<totalpT.x()/GeV<<", "<<totalpT.y()/GeV<<" )" << endl <<
"| originl enrgy = " << (thePlus + theMinus)/GeV/2.0 << endl <<
"| total energy = " << (plus + minus + minusMissing + plusMissing)/GeV/2.0 << endl <<
"| missing enrgy = " << (plus+minus+minusMissing+plusMissing-thePlus-theMinus)/GeV/2.0 << endl <<
"| allDipoles has " << allDipoles.size() << " dipoles." << endl <<
"| number of active dipoles: " << activeDipoles << endl <<
"| number of partons " << partons << endl <<
"| hasQuark: " << hasQuark << endl <<
"| number of interacting dipoles: " << interactingDips << endl <<
"| max pT : " << maxpT/GeV << endl <<
"| average dipole size: " << GeV*dipSize/double(activeDipoles) << endl;
cout << "----------------------------------------------------------------------\n";
if (ok ) cout << "| Found no errors! :)\n";
else cout << "| STATE NOT OK!! :o\n";
cout << "----------------------------------------------------------------------\n";
}
return ok;
}
bool DipoleState::forceSwing(PartonPtr p, double ymax1, double ymax2) {
DipolePtr d1 = p->dipoles().first;
DipolePtr d2 = p->dipoles().second;
if ( !d1 || !d2 )
cout << "forcing a swing for a parton that doesnt have 2 dipoles!" << endl;
//check for swings of same colour up to ymax1
d1->reset();
d1->touch();
bool foundFirstDipole = d1->forceGenerateRec(ymax1);
d2->reset();
d2->touch();
bool foundSecondDipole = d2->forceGenerateRec(ymax1);
if( foundFirstDipole || foundSecondDipole ) {
if ( !foundFirstDipole )
d2->recombine();
else if ( !foundSecondDipole )
d1->recombine();
else {
if ( d1->generatedY() < d2->generatedY() )
d1->recombine();
else
d2->recombine();
}
return true;
}
//check for swings of other colours up to ymax2
int trueColour = d1->colour();
int trueSys = d1->colourSystem();
for ( int c=0;c < Current<DipoleEventHandler>()->nColours();c++ ) {
d1->colour(c, trueSys);
d1->forceGenerateRec(ymax2);
}
d1->colour(trueColour);
foundFirstDipole = d1->swingDipole();
trueColour = d2->colour();
trueSys = d2->colourSystem();
for ( int c=0;c < Current<DipoleEventHandler>()->nColours();c++ ) {
d2->colour(c, trueSys);
d2->forceGenerateRec(ymax2);
}
d2->colour(trueColour);
foundSecondDipole = d2->swingDipole();
if( foundFirstDipole || foundSecondDipole ) {
if ( !foundFirstDipole )
d2->recombine();
else if ( !foundSecondDipole )
d1->recombine();
else {
if ( d1->generatedY() < d2->generatedY() )
d1->recombine();
else
d2->recombine();
return true;
}
}
return false;
}
DipoleStatePtr DipoleState::merge(DipoleStatePtr otherState) {
//reasign all other dipoles to belong to this state.
for( set< DipolePtr >::iterator it = otherState->allDipoles.begin();
it != otherState->allDipoles.end(); it++ ) {
(*it)->dipoleState(this);
}
//fill up initialDipoles and allDipoles with the other state.
theInitialDipoles.insert(theInitialDipoles.begin(),
otherState->initialDipoles().begin(),
otherState->initialDipoles().end());
theIncoming.insert(otherState->incoming().begin(), otherState->incoming().end());
allDipoles.insert(otherState->allDipoles.begin(),otherState->allDipoles.end());
theProducedParticles.insert(otherState->theProducedParticles.begin(),
otherState->theProducedParticles.end());
theIncomingShadows.insert(theIncomingShadows.end(),
otherState->theIncomingShadows.begin(),
otherState->theIncomingShadows.end());
//add up the original particles momenta
thePlus += otherState->thePlus;
theMinus += otherState->theMinus;
return this;
}
DipoleStatePtr DipoleState::collide( DipoleStatePtr otherState,
const vector<FList::const_iterator> & sel,
const ImpactParameters & b ) {
//mirror the other state in y
otherState->mirror(0.0);
//move the other state according to the impact parameter
otherState->translate(b);
//reasign all other dipoles to belong to this state.
for( set< DipolePtr >::iterator it = otherState->allDipoles.begin();
it != otherState->allDipoles.end(); it++ ) {
(*it)->dipoleState(this);
}
//fill up initialDipoles and allDipoles with the other state.
theInitialDipoles.insert(theInitialDipoles.begin(),
otherState->initialDipoles().begin(),
otherState->initialDipoles().end());
allDipoles.insert(otherState->allDipoles.begin(),otherState->allDipoles.end());
//decide if each collision in sel is elastic or not, and recouple the non-elastic ones
for( vector< FList::const_iterator>::const_iterator it = sel.begin();
it != sel.end(); it++ ) {
if( (*it)->first.second < 2.0*UseRandom::rnd() ) {
(*it)->second.first->swingDipole( (*it)->second.second );
(*it)->second.first->recombine();
}
}
return this;
}
void DipoleState::mirror(double y0) {
for ( set< DipolePtr >::iterator it = allDipoles.begin();
it != allDipoles.end(); it++) {
if ( !((*it)->children().first || (*it)->children().second) ) {
(*it)->partons().first->mirror(y0);
if ( !((*it)->neighbors().second) ) (*it)->partons().second->mirror(y0);
}
}
swap(thePlus, theMinus);
for ( map<PPtr, vector<PartonPtr> >::iterator it = theIncoming.begin();
it != theIncoming.end(); ++it )
it->first->setMomentum(lightCone(it->first->momentum().minus(),
it->first->momentum().plus()));
for ( int i = 0, N = theIncomingShadows.size(); i < N; ++i )
theIncomingShadows[i]->mirror(y0);
}
void DipoleState::translate(const ImpactParameters & b) {
for(set< DipolePtr >::iterator it = allDipoles.begin();
it != allDipoles.end(); it++) {
if( !((*it)->children().first || (*it)->children().second) ) {
b.translate((*it)->partons().first);
if( !((*it)->neighbors().second) )
b.translate((*it)->partons().second);
}
}
for ( int i = 0, N = theIncomingShadows.size(); i < N; ++i )
theIncomingShadows[i]->translate(b);
}
void DipoleState::fixValence(Step & step) const {
for ( map<PPtr, vector<PartonPtr> >::const_iterator it = theIncoming.begin();
it != theIncoming.end(); ++it ) {
tcWFInfoPtr wfi = WFInfo::getWFInfo(*it->first);
if ( !wfi ) continue;
PVector valence;
for ( int i = 0, N = it->second.size(); i < N; ++i ) {
valence.push_back(getParticle(it->second[i], true));
if ( !valence.back() ) valence.pop_back();
}
wfi->wf().fixValence(step, it->first, valence);
}
LorentzRotation Rshift = Utilities::boostToCM(step.collision()->incoming());
Utilities::transform(step.particles(), Rshift);
}
LorentzMomentum DipoleState::changeMass(LorentzMomentum p, Energy m, LorentzMomentum ref,
LorentzRotation * Rshift) {
Energy2 s = (p + ref).m2();
LorentzRotation R(0.0, 0.0, -(p.z() + ref.z())/(p.e() + ref.e()));
Energy z = Math::sign(SimplePhaseSpace::getMagnitude(s, ref.mt(), sqrt(sqr(m) + p.perp2())), (R*p).z());
LorentzMomentum pnew(p.x(), p.y(), z, sqrt(sqr(m) + sqr(z) + p.perp2()));
LorentzMomentum refnew(ref.x(), ref.y(), -z, sqrt(sqr(z) + ref.mt2()));
s = (pnew + refnew).m2();
LorentzRotation shift = LorentzRotation(0.0, 0.0, ref.z()/ref.e())*
LorentzRotation(0.0, 0.0, -refnew.z()/refnew.e());
pnew *= shift;
refnew *= shift;
if ( Rshift ) *Rshift = R.boostZ((pnew.z() + ref.z())/(pnew.e() + ref.e()));
return pnew;
}
vector<pair<Parton::Point, InvEnergy> > DipoleState::points() {
list< PartonPtr > partons = getPartons();
vector<pair<Parton::Point, InvEnergy> > ret;
for ( list<PartonPtr>::iterator it = partons.begin(); it != partons.end(); it++) {
InvEnergy r = max((*it)->dipoles().first->size(), (*it)->dipoles().first->size());
ret.push_back(make_pair((*it)->position(), r));
}
return ret;
}
list< PartonPtr > DipoleState::getPartons() const {
list<PartonPtr> partons;
for ( set<DipolePtr>::const_iterator it = allDipoles.begin() ; it != allDipoles.end(); it++ ) {
if ( (*it)->children().first || (*it)->children().second ) continue;
partons.push_back((*it)->partons().first);
if ( !((*it)->neighbors().second) ) //end of chain --> add also right parton
partons.push_back((*it)->partons().second);
}
return partons;
}
list< DipolePtr > DipoleState::getDipoles() const {
list<DipolePtr> dipoles;
for ( set<DipolePtr>::const_iterator it = allDipoles.begin() ; it != allDipoles.end(); it++ ) {
if ( (*it)->children().first || (*it)->children().second ) continue;
dipoles.push_back(*it);
}
return dipoles;
}
void DipoleState::makeOriginal() {
//clear the original dipoles
theInitialDipoles.clear(); //owns the dipoles
theSwingCandidates.clear(); //trancient
theTouchedSwingCandidates.clear(); //trancient
//go through allDipoles and set them to on or off.
for ( set<DipolePtr>::iterator it = allDipoles.begin();
it != allDipoles.end(); it++ ) {
DipolePtr d = *it;
if ( d->children().first || d->children().second ) d->isOn = false;
else d->isOn = true;
}
//go through allDipoles again
for ( set<DipolePtr>::iterator it = allDipoles.begin();
it != allDipoles.end(); it++ ) {
DipolePtr d = *it;
//remove neighbors that are off, both from dipole and parton
if ( d->isOn ) {
if ( d->neighbors().first && !d->neighbors().first->isOn ) {
d->neighbors(make_pair(DipolePtr(), d->neighbors().second));
d->partons().first->dipoles(make_pair(DipolePtr(), d->partons().first->dipoles().second) );
}
if ( d->neighbors().second && !d->neighbors().second->isOn ) {
d->neighbors(make_pair(d->neighbors().first, DipolePtr()));
d->partons().second->dipoles(make_pair(d->partons().second->dipoles().first, DipolePtr()) );
}
}
}
//go through allDipoles again
for ( set<DipolePtr>::iterator it = allDipoles.begin();
it != allDipoles.end(); ) {
DipolePtr d = *it;
//remove all pointers to off dipoles (not needed, as transient pointers)
//insert in original dipoles if the dipole is on.
if ( d->isOn ) {
theInitialDipoles.push_back(d);
d->partons().first->parents() = make_pair(PartonPtr(), PartonPtr());
d->partons().second->parents() = make_pair(PartonPtr(), PartonPtr());
d->partons().first->removeAllChildren();
d->partons().second->removeAllChildren();
d->reset();
it++;
}
//remove the dipole from allDipoles if it is off.
if ( !d->isOn ) {
set<DipolePtr>::iterator jt = it;
it++;
allDipoles.erase(jt);
}
}
//set parton data as if they were valence partons
list<PartonPtr> partons = getPartons();
int i = 0;
for( list<PartonPtr>::iterator it = partons.begin();
it != partons.end(); it++ ) {
PartonPtr p = *it;
p->number(i);
p->ordered(true);
p->valencePlus(p->plus());
p->valencePT(p->pT());
p->oY(p->y());
p->parents(make_pair(PartonPtr(), PartonPtr()));
i++;
}
}
vector<DipoleState::String> DipoleState::strings() {
set< list<PartonPtr> > colourLoops = loops();
vector<DipoleState::String> ret;
for ( set< list<PartonPtr> >::const_iterator loop = colourLoops.begin();
loop != colourLoops.end(); loop++ ) {
DipoleState::String s;
for ( list<PartonPtr>::const_iterator p = (*loop).begin(); p != (*loop).end(); p++ ) {
s.push_back( *p );
}
if ( Debug::level > 5 ) cout << "new string of size " << s.size() << endl;
ret.push_back( s );
}
return ret;
}
//Assume incoming particles massless for now. Valid as long higest pT
//is much larger than the mass.
void DipoleState::balanceMomenta() {
//Sum up how much p-/p+ is needed on the left/right side.
//Assumes that no particles got pushed over to the other side during reabsorbtion.
Energy PlmMissing = 0.0*GeV;
Energy PlmNeeded = 0.0*GeV;
Energy PrpMissing = 0.0*GeV;
Energy PrpNeeded = 0.0*GeV;
Energy Plp = 0.0*GeV;
Energy Prm = 0.0*GeV;
PartonPtr p;
for(set< DipolePtr >::iterator it = allDipoles.begin();
it != allDipoles.end(); it++) {
if( ((*it)->children().first || (*it)->children().second) )
continue;
//Do first parton always
p = (*it)->partons().first;
if( p->rightMoving() ) {
PlmMissing += p->minus();
PlmNeeded += p->pT().pt()*exp( p->y() );
Plp += p->plus();
}
else {
Prm += p->minus();
PrpMissing += p->plus();
PrpNeeded += p->pT().pt()*exp( -p->y() );
}
if( !((*it)->neighbors().second) ) {
//do second parton as well if at the end of a chain
p = (*it)->partons().second;
if( p->rightMoving() ) {
PlmMissing += p->minus();
PlmNeeded += p->pT().pt()*exp( p->y() );
Plp += p->plus();
}
else {
Prm += p->minus();
PrpMissing += p->plus();
PrpNeeded += p->pT().pt()*exp( -p->y() );
}
}
}
Energy PpShuffled = PrpNeeded - PrpMissing;
Energy PmShuffled = PlmNeeded - PlmMissing;
double a = -(1 + PmShuffled/Prm);
double b = -PlmNeeded/Prm;
double c = -(1 + PpShuffled/Plp);
double d = -PrpNeeded/Plp;
double e = (b - d)/c + a;
//how much p+/p- has to be scaled down on left/right side.
double q2 = - e/2 + sqrt( sqr(e)/4 + d*a/c );
double q1 = (d - c*q2)/q2;
if( max(abs(log(q1)),abs(log(q2))) > 0.1 ) {
cout << "yl shifted " << log(q2) << endl <<
"yr shifted " << log(q1) << endl;
}
if ( PlmMissing > Prm ) {
cout << "not enough p- in right state! VETO!" << endl;
}
else if ( PrpMissing > Plp ) {
cout << "not enough p+ in left state! VETO!" << endl;
}
//now shift all the partons in rapidity with ln q1 or ln q2 on left/right side.
for(set< DipolePtr >::iterator it = allDipoles.begin();
it != allDipoles.end(); it++) {
if( ((*it)->children().first || (*it)->children().second) )
continue;
//Do first parton always
p = (*it)->partons().first;
if( p->rightMoving() ) {
p->plus( p->plus()*q1 );
p->y( log(p->pT().pt()/p->plus() ) );
p->minus( p->pT().pt()*exp( p->y() ) );
}
else {
p->minus( p->minus()*q2 );
p->y( log( p->minus()/p->pT().pt() ) );
p->plus( p->pT().pt()*exp( -p->y() ) );
}
if( !((*it)->neighbors().second) ) {
//do second parton as well if at the end of a chain
p = (*it)->partons().second;
if( p->rightMoving() ) {
p->plus( p->plus()*q1 );
p->y( log(p->pT().pt()/p->plus() ) );
p->minus( p->pT().pt()*exp( p->y() ) );
}
else {
p->minus( p->minus()*q2 );
p->y( log( p->minus()/p->pT().pt() ) );
p->plus( p->pT().pt()*exp( -p->y() ) );
}
}
}
}
void DipoleState::saveGluonsToFile(double weight) const {
EventPtr event = handler().currentEvent();
string filename = handler().generator()->filename();
//set up outfile
ostringstream os;
PPair incoming = event->incoming();
LorentzPoint p1 = incoming.first->vertex();
LorentzPoint p2 = incoming.second->vertex();
double bx = (p1-p2).x()/femtometer;
double by = (p1-p2).y()/femtometer;
//TODO: interface!
os << "/home/christoffer/DIPSYevents/" << filename << "/event" << event->number() << ".dat";
// os << "/scratch/parton/christof/DIPSY/events/" << filename << "/event" << weight << ".dat";
weight = 2.0*int(sqrt(sqr(bx)+sqr(by))) + 1.0;
ofstream outfile (os.str().c_str());
//general information
outfile << "Generated with DIPSY 1.1" << endl;
outfile << "These are the real gluons only. No FSR or hadronisation has been done." << endl;
//print weight and impact parameter.
outfile << endl << "eventweight: " << weight << endl << endl;
outfile << "impact_parameter(fm):" << " " << bx << " " << by << endl;
vector<PartonPtr> partonVector;
int n = 0;
for(set< DipolePtr >::iterator it = allDipoles.begin();
it != allDipoles.end(); it++) {
if( ((*it)->children().first || (*it)->children().second) )
continue;
PartonPtr p = (*it)->partons().first;
if ( p->y() < -3 || p->y() > 3 ) continue;
partonVector.push_back(p);
p->number(n);
n++;
}
outfile << "number_of_particles:" << " " << partonVector.size() << endl;
outfile << "number of spectating nucleons: " << numberOfSpectators() << endl << endl;
//print gluons
outfile << "particle_number" << " "
<< "transverse_position_x(fm)" << " "
<< "transverse_position_y(fm)" << " "
<< "transverse_momentum_x(GeV)" << " "
<< "transverse_momentum_y(GeV)" << " "
<< "rapidity" << " "
<< "colour_neighbour_number" << " "
<< "anticolour_neighbour_number" << endl << endl;
for(int i = 0; i < int(partonVector.size()); i++) {
PartonPtr p = partonVector[i];
outfile << p->number() << " "
<< p->position().x()*hbarc/femtometer << " "
<< p->position().y()*hbarc/femtometer << " "
<< p->pT().x()/GeV << " "
<< p->pT().y()/GeV << " "
<< p->y() << " ";
if ( p->dipoles().first )
outfile << p->dipoles().first->partons().first->number() << " ";
else outfile << -1 << " ";
if ( p->dipoles().second )
outfile << p->dipoles().second->partons().second->number() << " ";
else outfile << -1 << " ";
outfile << endl;
}
outfile.close();
cout << "printed gluons to file " << os.str().c_str() << endl;
}
int DipoleState::numberOfSpectators() const {
int nSpectators = 0;
for ( int i = 0; i < int(theInitialDipoles.size()); i++ ) {
tPartonPtr p1 = theInitialDipoles[i]->partons().first;
tPartonPtr p2 = theInitialDipoles[i]->partons().second;
tPartonPtr p3;
if ( i != int(theInitialDipoles.size())-1 &&
theInitialDipoles[i+1]->partons().first == p2 )
p3 = theInitialDipoles[i+1]->partons().second;
else if ( i != int(theInitialDipoles.size())-1 &&
theInitialDipoles[i+1]->partons().second == p1 )
p3 = theInitialDipoles[i+1]->partons().first;
else if ( i != 0 && theInitialDipoles[i-1]->partons().first == p2 )
p3 = theInitialDipoles[i-1]->partons().second;
else if ( i != 0 && theInitialDipoles[i-1]->partons().second == p1 )
p3 = theInitialDipoles[i-1]->partons().first;
if ( !(p1->interacted()) &&
!(p2->interacted()) &&
!( p3 && p3->interacted() ) ) {
nSpectators++;
}
}
return nSpectators/3;
}
void DipoleState::unifyColourSystems(int sys) {
for ( set<DipolePtr>::iterator it = allDipoles.begin();
it != allDipoles.end(); ++it ) (**it).colourSystem(sys);
}
tPPtr DipoleState::getParticle(tcPartonPtr parton, bool nocreate) const {
map<tcPartonPtr,PPtr>::iterator pit = theProducedParticles.find(parton);
if ( pit != theProducedParticles.end() ) return pit->second;
if ( nocreate ) return tPPtr();
if( parton->plus() <= 0.0*GeV ) cout << "p+ is not positive :(" << endl <<
"p+ = " << parton->plus()/GeV << endl <<
"p- = " << parton->minus()/GeV << endl <<
"pt = " << parton->pT().pt()/GeV << endl <<
"y = " << parton->y() << endl;
PPtr p = CurrentGenerator::current().getParticle(parton->flavour());
p->setMomentum(parton->momentum());
p->setVertex(LorentzPoint(parton->position().x()*hbarc, parton->position().y()*hbarc,
ZERO, ZERO));
p->getInfo().push_back(new_ptr(ParticleInfo(parton)));
theProducedParticles[parton] = p;
return p;
}
void DipoleState::setupShadows() {
thePrimaryShadow = tcSPartonPtr();
set<tPartonPtr> valence;
for ( int i = 0, N = initialDipoles().size(); i < N ; ++i ) {
valence.insert(initialDipoles()[i]->partons().first);
valence.insert(initialDipoles()[i]->partons().second);
}
for ( set<tPartonPtr>::iterator it = valence.begin(); it != valence.end(); ++it )
theIncomingShadows.push_back(ShadowParton::createValence(**it));
}
void DipoleState::flagShadowsOnShell(tcSPartonPtr valence) {
for ( int i = 0, N = theIncomingShadows.size(); i < N; ++i )
if ( theIncomingShadows[i] != valence )
theIncomingShadows[i]->onShell(true);
}
void DipoleState::propagateShadowMomenta() {
for ( int i = 0, N = theIncomingShadows.size(); i < N; ++i )
theIncomingShadows[i]->propagateShadowMomenta();
}
void DipoleState::resetShadows() {
for ( int i = 0, N = theIncomingShadows.size(); i < N; ++i )
theIncomingShadows[i]->reset();
}
void DipoleState::resetInteractedShadows() {
for ( int i = 0, N = theIncomingShadows.size(); i < N; ++i )
theIncomingShadows[i]->resetInteracted();
}
void DipoleState::performShadowInteraction() {
for ( int i = 0, N = theIncomingShadows.size(); i < N; ++i )
if ( theIncomingShadows[i] != thePrimaryShadow )
theIncomingShadows[i]->setOnShell();
}
LorentzMomentum DipoleState::incomingMomentum(tcSPartonPtr valence, bool final) {
- if ( thePrimaryShadow != valence ) return valence->momentum();
+ if ( thePrimaryShadow && thePrimaryShadow != valence ) return valence->momentum();
LorentzMomentum ret = lightCone(plus(), minus(), TransverseMomentum());
for ( int i = 0, N = theIncomingShadows.size(); i < N; ++i )
if ( theIncomingShadows[i] != valence ) {
ret -= theIncomingShadows[i]->momentum();
if ( final) theIncomingShadows[i]->setOnShell();
}
return ret;
}
void DipoleState::acceptShadowInteraction(tcSPartonPtr valence) {
if ( thePrimaryShadow ) return;
thePrimaryShadow = valence;
}
void DipoleState::debugShadowTree() {
cerr << "=== Shadow Tree Structure ===" << endl;
for ( int i = 0, N = theIncomingShadows.size(); i < N; ++i )
theIncomingShadows[i]->debugTree("");
}
void DipoleState::persistentOutput(PersistentOStream & os) const {
os << ounit(thePlus, GeV) << ounit(theMinus, GeV)
<< ounit(theMinusDeficit, GeV) << theHandler << theInitialDipoles
<< theSwingCandidates << theTouchedSwingCandidates
<< theWFInfo << theWeight << doTakeHistory
<< theYmax << ounit(theCollidingEnergy, GeV) << allDipoles << theIncoming
<< theProducedParticles << theIncomingShadows;
}
void DipoleState::persistentInput(PersistentIStream & is, int) {
is >> iunit(thePlus, GeV)>> iunit(theMinus, GeV)
>> iunit(theMinusDeficit, GeV) >> theHandler >> theInitialDipoles
>> theSwingCandidates >> theTouchedSwingCandidates
>> theWFInfo >> theWeight >> doTakeHistory
>> theYmax >> iunit(theCollidingEnergy, GeV) >> allDipoles >> theIncoming
>> theProducedParticles >> theIncomingShadows;
}
// Static variable needed for the type description system in ThePEG.
#include "ThePEG/Utilities/DescribeClass.h"
DescribeClass<DipoleState,PersistentBase>
describeDIPSYDipoleState("DIPSY::DipoleState", "libAriadne5.so libDIPSY.so");
void DipoleState::Init() {}
diff --git a/DIPSY/DipoleXSec.cc b/DIPSY/DipoleXSec.cc
--- a/DIPSY/DipoleXSec.cc
+++ b/DIPSY/DipoleXSec.cc
@@ -1,1855 +1,1856 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DipoleXSec class.
//
#include "DipoleXSec.h"
#include "Dipole.h"
#include "DipoleState.h"
#include "Parton.h"
#include "DipoleEventHandler.h"
#include "RealParton.h"
#include "RealPartonState.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/Current.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/DebugItem.h"
#include "ThePEG/Utilities/Throw.h"
#ifdef ThePEG_TEMPLATES_IN_CC_FILE
// #include "DipoleXSec.tcc"
#endif
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "gsl/gsl_sf_bessel.h"
using namespace ThePEG;
using namespace DIPSY;
DipoleXSec::~DipoleXSec() {}
DipoleInteraction::DipoleInteraction(Dipole & dlin, Dipole & drin,
const ImpactParameters & bin)
: dips(&dlin, &drin), dnext(dlin.neighbors().first, drin.neighbors().first),
dprev(dlin.neighbors().second, drin.neighbors().second),
b(&bin), sints(tSPartonPtr(), tSPartonPtr()),
- f2(0.0), uf2(0.0), norec(false), kt(ZERO) {
+ f2(0.0), uf2(0.0), norec(false), kt(ZERO), id(0) {
// If different colours interaction is either between c-c or
// cbar-cbar. If the same colour we only haver c-cbar.
ints = make_pair(dlin.partons().first, drin.partons().first);
spec = make_pair(dlin.partons().second, drin.partons().second);
if ( dlin.colour() == drin.colour() ) swap(ints.second, spec.second);
if ( UseRandom::rndbool() ) {
swap(ints.first, spec.first);
swap(ints.second, spec.second);
}
Parton::Point r = bin.difference(ints.first->position(), ints.second->position());
d2 = min(min(r.pt2(), dlin.size2()), drin.size2());
kt = 1.0/sqrt(d2);
rec = -Current<DipoleEventHandler>()->emitter().pTScale()*r/d2;
if ( ints.first->shadow() ) {
sints = make_pair(ints.first->shadow()->resolve(d2, spec.first),
ints.second->shadow()->resolve(d2, spec.second));
}
}
void DipoleInteraction::prepare() const {
sints = make_pair(ints.first->shadow()->resolve(d2, spec.first),
ints.second->shadow()->resolve(d2, spec.second));
sints.first->insertInteraction(this);
sints.second->insertInteraction(this);
}
bool DipoleInteraction::check(bool final) const {
sints = make_pair(ints.first->shadow()->resolve(d2, spec.first),
ints.second->shadow()->resolve(d2, spec.second));
ShadowParton::Propagator ppi = ints.first->shadow()->propagator(d2, spec.first, this, final);
ShadowParton::Propagator ppj = ints.second->shadow()->propagator(d2, spec.second, this, final);
if ( ppi.fail || ppj.fail ) return false;
LorentzMomentum pi = ppi.p;
LorentzMomentum pj = ppj.p;
TransverseMomentum rrec;
if ( !norec ) rrec = rec;
TransverseMomentum pti = TransverseMomentum(pi.x(), pi.y()) + rrec;
Energy2 mti2 = pti.pt2() + sqr(sints.first->mass());
TransverseMomentum ptj = TransverseMomentum(pj.x(), pj.y()) + b->invRotatePT(-rrec);
Energy2 mtj2 = ptj.pt2() + sqr(sints.second->mass());
Energy PP = pi.plus() + pj.minus();
Energy PM = pi.minus() + pj.plus();
- if ( PP <= min(ppi.colpos, ppi.acopos) ||
- PM <= min(ppj.colpos, ppj.acopos) ) return false;
+ if ( PP > min(ppi.colpos, ppi.acopos) ||
+ PM > min(ppj.colpos, ppj.acopos) ) return false;
if ( PP*PM + mti2 - mtj2 < ZERO ||
PP*PM + mtj2 - mti2 < ZERO ) return false;
Energy2 sqri = (sqr(PP*PM + mti2 - mtj2) - 4.0*mti2*PM*PP)/sqr(PM);
Energy2 sqrj = (sqr(PP*PM + mtj2 - mti2) - 4.0*mtj2*PM*PP)/sqr(PP);
if ( sqri < ZERO || sqrj < ZERO ) return false;
+ if ( !id ) return true;
sints.first->pTplus(pti, 0.5*(PP + mti2/PM - mtj2/PM + sqrt(sqri)));
- sints.first->interacting(true);
+ sints.first->interacting(id);
+ sints.first->acceptInteraction(this);
sints.second->pTplus(ptj, 0.5*(PM + mtj2/PP - mti2/PP + sqrt(sqrj)));
- sints.second->interacting(true);
+ sints.second->interacting(id);
+ sints.second->acceptInteraction(this);
if ( final ) {
sints.first->setOnShell(true);
- ints.first->interact(true);
+ sints.first->original()->interact(true);
sints.second->setOnShell(true);
- ints.second->interact(true);
+ sints.second->original()->interact(true);
}
return true;
}
void DipoleInteraction::reject() const {
sints.first->rejectInteraction(this);
sints.second->rejectInteraction(this);
}
void DipoleInteraction::accept() const {
sints = make_pair(ints.first->shadow()->resolve(d2, spec.first),
ints.second->shadow()->resolve(d2, spec.second));
- sints.first->acceptInteraction(this);
- sints.second->acceptInteraction(this);
}
InvEnergy DipoleXSec::rMax() const {
return theRMax > 0.0*InvGeV? theRMax: Current<DipoleEventHandler>()->rMax();
}
double DipoleXSec::fSinFn(const Parton::Point & rho1, const Parton::Point & rho2,
const TransverseMomentum & pt) const {
if ( sinFunction == 1 ) {
double r1 = rho1.pt2()*pt.pt2();
double r2 = rho2.pt2()*pt.pt2();
return r1*r2/(4.0*(r1 + 1.0)*(r2 + 1.0));
}
double s1 = pt.x()*rho1.x() + pt.y()*rho1.y();
double s2 = pt.x()*rho2.x() + pt.y()*rho2.y();
return sqr(sin(s1)*sin(s2));
}
double DipoleXSec::
fij(const pair<tPartonPtr, tPartonPtr> left,
const pair<tPartonPtr, tPartonPtr> right,
const ImpactParameters & b, bool veto) const {
Nfij++;
tcPartonPtr p11 = left.first;
tcPartonPtr p12 = left.second;
tcPartonPtr p21 = right.first;
tcPartonPtr p22 = right.second;
//TODO: keep only interaction 0, as that is the only one supported anyway
if ( theInteraction == 0 || theInteraction == 1 || theInteraction == 3 || theInteraction == 4 ) {
InvEnergy2 rr11 = b.dist2(*p11, *p21);
InvEnergy2 rr21 = b.dist2(*p12, *p21);
InvEnergy2 rr12 = b.dist2(*p11, *p22);
InvEnergy2 rr22 = b.dist2(*p12, *p22);
if ( veto && kinematicsVeto(left, right, b) ) {
return 0.0;
}
TransverseMomentum pt;
double pTScale = Current<DipoleEventHandler>()->emitter().pTScale();
InvEnergy rscale = sqrt(min(min(min(rr12, rr21),min(rr11, rr22)),
min(p11->dist2(*p12), p21->dist2(*p22))));
double fudgeME = 1.0;
if ( theInteraction == 0 ) {
pair<bool, bool> ints;
ints = int0Partons(p11, p12, p21, p22, b);
Parton::Point r1 = b.difference((ints.first ? p11->position():p12->position()),
(ints.second ? p21->position():p22->position()));
if ( r1.pt2() < min(p11->dist2(*p12), p21->dist2(*p22)) &&
Current<DipoleEventHandler>()->fudgeME() ) {
double deltay = ( ints.first? p11->y(): p12->y() ) +
( ints.second? p21->y(): p22->y() );
fudgeME = 1.0 - 1.0/(1.0 + cosh(deltay));
fudgeME *= Current<DipoleEventHandler>()->fudgeFactorME();
}
pt = pTScale*r1/r1.pt2();
rscale = sqrt(min(min(r1.pt2(),p11->dist2(*p12)), p21->dist2(*p22)));
}
if ( theInteraction == 1 ) { //4 parton distance
Parton::Point r1 = b.difference(p11->position(), p21->position());
Parton::Point r2 = b.difference(p11->position(), p22->position());
Parton::Point r3 = b.difference(p12->position(), p21->position());
Parton::Point r4 = b.difference(p12->position(), p22->position());
pt = pTScale*(r1/sqr(r1.pt()) + r2/sqr(r2.pt()) +
r3/sqr(r3.pt()) + r4/sqr(r4.pt()))*0.35;
}
if ( theInteraction == 3 ) { //2 parton distance
Parton::Point r1 = b.difference(p11->position(), p22->position());
Parton::Point r2 = b.difference(p12->position(), p21->position());
pt = pTScale*(r1/sqr(r1.pt()) + r2/sqr(r2.pt()))*0.6;
}
Parton::Point rho1 = (p12->position() - p11->position())/2.0;
Parton::Point rho2 = b.rotate(p22->position() - p21->position())/2.0;
return 8.0*fudgeME*sqr(Current<DipoleEventHandler>()->alphaS(rscale))*
fSinFn(rho1, rho2, pt)*
sqr(sqr(pt.pt()))/sqr(sqr(pt.pt()) + sqr(pTScale/rMax()));
}
return 0.0;
}
double DipoleXSec::
fij(const Dipole & dleft, const Dipole & dright,
const ImpactParameters & b, bool veto) const {
pair<tPartonPtr, tPartonPtr> left = dleft.partons();
pair<tPartonPtr, tPartonPtr> right = dright.partons();
pair<bool, bool> ints;
Nfij++;
double fudgeME = 1.0;
tcPartonPtr p11 = left.first;
tcPartonPtr p12 = left.second;
tcPartonPtr p21 = right.first;
tcPartonPtr p22 = right.second;
ints = int0Partons(p11, p12, p21, p22, b);
//TODO: keep only interaction 0, as that is the only one supported anyway
if ( theInteraction == 0 || theInteraction == 1 || theInteraction == 3 || theInteraction == 4 ) {
InvEnergy2 rr11 = b.dist2(*p11, *p21);
InvEnergy2 rr21 = b.dist2(*p12, *p21);
InvEnergy2 rr12 = b.dist2(*p11, *p22);
InvEnergy2 rr22 = b.dist2(*p12, *p22);
if ( veto && kinematicsVeto(dleft, dright, b, ints) ) {
return 0.0;
}
TransverseMomentum pt;
double pTScale = Current<DipoleEventHandler>()->emitter().pTScale();
InvEnergy rscale = sqrt(min(min(min(rr12, rr21),min(rr11, rr22)),
min(p11->dist2(*p12), p21->dist2(*p22))));
if ( theInteraction == 0 ) {
Parton::Point r1 = b.difference((ints.first ? p11->position():p12->position()),
(ints.second ? p21->position():p22->position()));
if ( r1.pt2() < min(p11->dist2(*p12), p21->dist2(*p22)) &&
Current<DipoleEventHandler>()->fudgeME() ) {
double deltay = ( ints.first? p11->y(): p12->y() ) -
( ints.second? p21->y(): p22->y() );
fudgeME = 1.0 - 1.0/(1.0 + cosh(deltay));
}
pt = pTScale*r1/sqr(r1.pt());
rscale = sqrt(min(min(sqr(r1.pt()),p11->dist2(*p12)), p21->dist2(*p22)));
}
if ( theInteraction == 1 ) { //4 parton distance
Parton::Point r1 = b.difference(p11->position(), p21->position());
Parton::Point r2 = b.difference(p11->position(), p22->position());
Parton::Point r3 = b.difference(p12->position(), p21->position());
Parton::Point r4 = b.difference(p12->position(), p22->position());
pt = pTScale*(r1/sqr(r1.pt()) + r2/sqr(r2.pt()) +
r3/sqr(r3.pt()) + r4/sqr(r4.pt()))*0.35;
}
if ( theInteraction == 3 ) { //2 parton distance
Parton::Point r1 = b.difference(p11->position(), p22->position());
Parton::Point r2 = b.difference(p12->position(), p21->position());
pt = pTScale*(r1/sqr(r1.pt()) + r2/sqr(r2.pt()))*0.6;
}
Parton::Point rho1 = (p12->position() - p11->position())/2.0;
Parton::Point rho2 = b.rotate(p22->position() - p21->position())/2.0;
return 8.0*fudgeME*sqr(Current<DipoleEventHandler>()->alphaS(rscale))*
fSinFn(rho1, rho2, pt)*
sqr(sqr(pt.pt()))/sqr(sqr(pt.pt()) + sqr(pTScale/rMax()));
}
return 0.0;
}
DipoleInteraction DipoleXSec::
fij(const ImpactParameters & b, Dipole & dleft, Dipole & dright,
bool veto) const {
DipoleInteraction di(dleft, dright, b);
Nfij++;
double fudgeME = 1.0;
if ( veto && kinematicsVeto(di) ) return di;
Energy m0 = Current<DipoleEventHandler>()->emitter().pTScale()/rMax();
Parton::Point r1 = b.difference(di.ints.first->position(), di.ints.second->position());
if ( r1.pt2() < min(dleft.size2(), dright.size2()) &&
Current<DipoleEventHandler>()->fudgeME() ) {
double deltay = di.ints.first->y() + di.ints.second->y();
fudgeME = 1.0 - 1.0/(1.0 + cosh(deltay));
fudgeME *= Current<DipoleEventHandler>()->fudgeFactorME();
}
if ( Current<DipoleEventHandler>()->fudgeME() > 1 ) {
if ( dleft.colour() == dright.colour() ) fudgeME *= 9.0/2.0;
else fudgeME *= 9.0/16.0;
}
Parton::Point rho1 = di.dips.first->vSize()/2.0;
Parton::Point rho2 = b.rotate(di.dips.second->vSize())/2.0;
di.f2 = 16.0*fudgeME*sqr(Current<DipoleEventHandler>()->alphaS(sqrt(di.d2)))*
fSinFn(rho1, rho2, di.rec)*sqr(di.rec.pt2())/sqr(di.rec.pt2() + sqr(m0));
di.uf2 = unitarize(di.f2);
return di;
}
bool DipoleXSec::kinematicsVeto(const pair<tPartonPtr, tPartonPtr> left,
const pair<tPartonPtr, tPartonPtr> right,
const ImpactParameters & b) const {
pair<pair<bool, bool>, pair<bool, bool> > ints = doesInt(left, right, b);
InteractionRecoil recs = recoil(left, right, b, ints);
if ( theInteraction == 0 ) {
//the key partons for f_ij
pair<bool, bool> ints0 = int0Partons(left.first, left.second, right.first, right.second, b);
//first set up effective partons with range etc
tPartonPtr p1 = (ints0.first ? left.first:left.second);
tPartonPtr p2 = (ints0.second ? right.first:right.second);
//there MAY be secondary partons in the interaction, decided by ints.
tPartonPtr p1sec, p2sec;
if ( ints.first.first && ints.first.second )
p1sec = (ints0.first ? left.second:left.first);
if ( ints.second.first && ints.second.second )
p2sec = (ints0.second ? right.second:right.first);
//only distance between key partons will affect range (reasonable?)
InvEnergy range1 = sqrt(min(left.first->dist2(*left.second)/4.0, b.dist2(*p1, *p2)));
InvEnergy range2 = sqrt(min(right.first->dist2(*right.second)/4.0, b.dist2(*p1, *p2)));
if ( Current<DipoleEventHandler>()->emitter().rangeMode() == 1 ) {
range1 = sqrt(left.first->dist2(*left.second)/4.0);
range2 = sqrt(right.first->dist2(*right.second)/4.0);
}
EffectivePartonPtr ep1 = EffectiveParton::create(*p1, range1);
EffectivePartonPtr ep2 = EffectiveParton::create(*p2, range2);
EffectivePartonPtr ep1sec, ep2sec;
if ( p1sec ) ep1sec = EffectiveParton::create(*p1sec, range1);
if ( p2sec ) ep2sec = EffectiveParton::create(*p2sec, range2);
TransverseMomentum rec1 = (ints0.first ? recs.first.first:recs.first.second);
TransverseMomentum rec2 = (ints0.second ? recs.second.first:recs.second.second);
Energy pt1 = (ep1->pT() + rec1).pt();
Energy minus1 = sqr(pt1)/ep1->plus();
Energy pt2 = (ep2->pT() + rec2).pt();
Energy minus2 = sqr(pt2)/ep2->plus();
//sum up total supplied and needed LC momentum from each side
Energy leftPlus = ep1->plus() + (p1sec ? ep1sec->plus():ZERO);
Energy leftMinus = minus1;
if ( p1sec ) {
TransverseMomentum rec1sec = (ints0.first ? recs.first.second:recs.first.first);
Energy pt1sec = (ep1sec->pT() + rec1sec).pt();
Energy minus1sec = sqr(pt1sec)/ep1sec->plus();
leftMinus += minus1sec;
}
Energy rightPlus = ep2->plus() + (p2sec ? ep2sec->plus():ZERO);
Energy rightMinus = minus2 + (p2sec ? ep2sec->minus():ZERO);
if ( p2sec ) {
TransverseMomentum rec2sec = (ints0.second ? recs.second.second:recs.second.first);
Energy pt2sec = (ep2sec->pT() + rec2sec).pt();
Energy minus2sec = sqr(pt2sec)/ep2sec->plus();
rightMinus += minus2sec;
}
if ( theIntOrdering == 2 ) {
Energy maxRec = max(max(recs.first.first.pt() ,recs.first.second.pt()),
max(recs.second.first.pt() ,recs.second.second.pt()));
if ( leftPlus*rightPlus < 16.0*sqr(maxRec) ) return true;
else return false;
}
//check enough energy to set all on shell
if ( leftPlus*rightPlus < 16.0*leftMinus*rightMinus ) {
return true;
}
//take LC momentum transfers in account
Energy plus1 = ep1->plus()*(1.0 - rightMinus/leftPlus);
minus1 /= 1.0 - rightMinus/leftPlus;
Energy plus2 = ep2->plus()*(1.0 - leftMinus/rightPlus);
minus2 /= 1.0 - leftMinus/rightPlus;
//check ordering of the key partons (secondaries can be unordered if they want)
double PSInf = Current<DipoleEventHandler>()->emitter().PSInflation();
if ( theIntOrdering == 0 ) {
//just check plus and minus ordering after recoils
if ( plus1*PSInf < minus2 ) return true;
if ( plus2*PSInf < minus1 ) return true;
}
else if ( theIntOrdering == 1 ) {
//check p- ordering from both sides, as if it was evolution.
//that is, as if no future recoils.
double PMOrd = Current<DipoleEventHandler>()->emitter().PMinusOrdering();
TransverseMomentum rec1 = (ints0.first ? recs.first.first:recs.first.second);
TransverseMomentum rec2 = (ints0.second ? recs.second.first:recs.second.second);
if ( sqr(max(rec1.pt(), rec2.pt())*PSInf) < minus1*minus2*PMOrd )
return true;
}
//check that the partons stay on their side, to avoid doublecounting
double yInt = 0;
if ( p1->dipoles().first ) yInt = p1->dipoles().first->dipoleState().ymax();
else if ( p1->dipoles().second ) yInt = p1->dipoles().second->dipoleState().ymax();
double y1 = 0.5*log(minus1/plus1);
if ( y1 > yInt ) return true;
double y2 = 0.5*log(plus2/minus2);
if ( y2 < yInt ) return true;
//should we check secondaries as well here?
//if no veto triggered, it should not be vetoed
return false;
}
else {
cerr << "only interaction 0 is supported in kinematics veto atm, sorry" << endl;
}
return false;
}
bool DipoleXSec::kinematicsVeto(const Dipole & dleft,
const Dipole & dright,
const ImpactParameters & b,
const pair<bool,bool> & ints0) const {
pair<tPartonPtr, tPartonPtr> left = dleft.partons();
pair<tPartonPtr, tPartonPtr> right = dright.partons();
pair<pair<bool, bool>, pair<bool, bool> > ints = doesInt(left, right, b);
InteractionRecoil recs = recoil(left, right, b, ints, ints0);
if ( theInteraction == 0 ) {
//first set up effective partons with range etc
tPartonPtr p1 = (ints0.first ? left.first:left.second);
tPartonPtr p2 = (ints0.second ? right.first:right.second);
//there MAY be secondary partons in the interaction, decided by ints.
tPartonPtr p1sec, p2sec;
if ( ints.first.first && ints.first.second )
p1sec = (ints0.first ? left.second:left.first);
if ( ints.second.first && ints.second.second )
p2sec = (ints0.second ? right.second:right.first);
//only distance between key partons will affect range (reasonable?)
InvEnergy range1 = sqrt(min(left.first->dist2(*left.second)/4.0, b.dist2(*p1, *p2)));
InvEnergy range2 = sqrt(min(right.first->dist2(*right.second)/4.0, b.dist2(*p1, *p2)));
if ( Current<DipoleEventHandler>()->emitter().rangeMode() == 1 ) {
range1 = sqrt(left.first->dist2(*left.second)/4.0);
range2 = sqrt(right.first->dist2(*right.second)/4.0);
}
EffectivePartonPtr ep1 = dleft.getEff(p1, range1);
EffectivePartonPtr ep2 = dright.getEff(p2, range2);
EffectivePartonPtr ep1sec, ep2sec;
if ( p1sec ) ep1sec = dleft.getEff(p1sec, range1);
if ( p2sec ) ep2sec = dright.getEff(p2sec, range2);
TransverseMomentum rec1 = (ints0.first ? recs.first.first:recs.first.second);
TransverseMomentum rec2 = (ints0.second ? recs.second.first:recs.second.second);
Energy pt1 = (ep1->pT() + rec1).pt();
Energy minus1 = sqr(pt1)/ep1->plus();
Energy pt2 = (ep2->pT() + rec2).pt();
Energy minus2 = sqr(pt2)/ep2->plus();
//sum up total supplied and needed LC momentum from each side
Energy leftPlus = ep1->plus() + (p1sec ? ep1sec->plus():ZERO);
Energy leftMinus = minus1;
if ( p1sec ) {
TransverseMomentum rec1sec = (ints0.first ? recs.first.second:recs.first.first);
Energy pt1sec = (ep1sec->pT() + rec1sec).pt();
Energy minus1sec = sqr(pt1sec)/ep1sec->plus();
leftMinus += minus1sec;
}
Energy rightPlus = ep2->plus() + (p2sec ? ep2sec->plus():ZERO);
Energy rightMinus = minus2 + (p2sec ? ep2sec->minus():ZERO);
if ( p2sec ) {
TransverseMomentum rec2sec = (ints0.second ? recs.second.second:recs.second.first);
Energy pt2sec = (ep2sec->pT() + rec2sec).pt();
Energy minus2sec = sqr(pt2sec)/ep2sec->plus();
rightMinus += minus2sec;
}
if ( theIntOrdering == 2 ) {
Energy maxRec = max(max(recs.first.first.pt() ,recs.first.second.pt()),
max(recs.second.first.pt() ,recs.second.second.pt()));
if ( leftPlus*rightPlus < 16.0*sqr(maxRec) ) return true;
else return false;
}
//check enough energy to set all on shell
if ( leftPlus*rightPlus < 16.0*leftMinus*rightMinus ) {
return true;
}
//take LC momentum transfers in account
Energy plus1 = ep1->plus()*(1.0 - rightMinus/leftPlus);
minus1 /= 1.0 - rightMinus/leftPlus;
Energy plus2 = ep2->plus()*(1.0 - leftMinus/rightPlus);
minus2 /= 1.0 - leftMinus/rightPlus;
//check ordering of the key partons (secondaries can be unordered if they want)
double PSInf = Current<DipoleEventHandler>()->emitter().PSInflation();
if ( theIntOrdering == 0 ) {
//just check plus and minus ordering after recoils
if ( plus1*PSInf < minus2 ) return true;
if ( plus2*PSInf < minus1 ) return true;
}
else if ( theIntOrdering == 1 ) {
//check p- ordering from both sides, as if it was evolution.
//that is, as if no future recoils.
double PMOrd = Current<DipoleEventHandler>()->emitter().PMinusOrdering();
TransverseMomentum rec1 = (ints0.first ? recs.first.first:recs.first.second);
TransverseMomentum rec2 = (ints0.second ? recs.second.first:recs.second.second);
if ( sqr(max(rec1.pt(), rec2.pt())*PSInf) < minus1*minus2*PMOrd )
return true;
}
//check that the partons stay on their side, to avoid doublecounting
double yInt = 0;
if ( p1->dipoles().first ) yInt = p1->dipoles().first->dipoleState().ymax();
else if ( p1->dipoles().second ) yInt = p1->dipoles().second->dipoleState().ymax();
double y1 = 0.5*log(minus1/plus1);
if ( y1 > yInt ) return true;
double y2 = 0.5*log(plus2/minus2);
if ( y2 < yInt ) return true;
//should we check secondaries as well here?
//if no veto triggered, it should not be vetoed
return false;
}
else {
cerr << "only interaction 0 is supported in kinematics veto atm, sorry" << endl;
}
return false;
}
bool DipoleXSec::kinematicsVeto(const DipoleInteraction & di) const {
if ( di.sints.first ) return !di.check(); // Hey! We're using shadows!
ImpactParameters b = *di.b;
const Dipole & dleft = *di.dips.first;
const Dipole & dright = *di.dips.second;
//first set up effective partons with range etc
tcPartonPtr p1 = di.ints.first;
tcPartonPtr p2 = di.ints.second;
//only distance between key partons will affect range (reasonable?)
InvEnergy range1 = sqrt(min(dleft.size2()/4.0, b.dist2(*p1, *p2)));
InvEnergy range2 = sqrt(min(dright.size2()/4.0, b.dist2(*p1, *p2)));
EffectivePartonPtr ep1 = dleft.getEff(p1, range1);
EffectivePartonPtr ep2 = dright.getEff(p2, range2);
TransverseMomentum rec1 = di.rec;
TransverseMomentum rec2 = di.b->invRotatePT(-di.rec);
Energy pt1 = (ep1->pT() + rec1).pt();
Energy minus1 = sqr(pt1)/ep1->plus();
Energy pt2 = (ep2->pT() + rec2).pt();
Energy minus2 = sqr(pt2)/ep2->plus();
//sum up total supplied and needed LC momentum from each side
Energy leftPlus = ep1->plus();
Energy leftMinus = minus1;
Energy rightPlus = ep2->plus();
Energy rightMinus = minus2;
if ( theIntOrdering == 2 ) {
if ( leftPlus*rightPlus < 16.0*rec1.pt2() ) return true;
else return false;
}
//check enough energy to set all on shell
if ( leftPlus*rightPlus < 16.0*leftMinus*rightMinus ) {
return true;
}
//take LC momentum transfers in account
Energy plus1 = ep1->plus()*(1.0 - rightMinus/leftPlus);
minus1 /= 1.0 - rightMinus/leftPlus;
Energy plus2 = ep2->plus()*(1.0 - leftMinus/rightPlus);
minus2 /= 1.0 - leftMinus/rightPlus;
//check ordering of the key partons (secondaries can be unordered if they want)
double PSInf = Current<DipoleEventHandler>()->emitter().PSInflation();
if ( theIntOrdering == 0 ) {
//just check plus and minus ordering after recoils
if ( plus1*PSInf < minus2 ) return true;
if ( plus2*PSInf < minus1 ) return true;
}
//check that the partons stay on their side, to avoid doublecounting
double yInt = 0;
if ( p1->dipoles().first ) yInt = p1->dipoles().first->dipoleState().ymax();
else if ( p1->dipoles().second ) yInt = p1->dipoles().second->dipoleState().ymax();
double y1 = 0.5*log(minus1/plus1);
if ( y1 > yInt ) return true;
double y2 = 0.5*log(plus2/minus2);
if ( y2 < yInt ) return true;
//should we check secondaries as well here?
//if no veto triggered, it should not be vetoed
return false;
return false;
}
double DipoleXSec::
sumf(const DipoleState & sl, const DipoleState & sr,
const ImpactParameters & b) const {
Nfij = 0;
NBVeto = 0;
NScalVeto = 0;
scalVeto = 0.0;
bVeto = 0.0;
vector<tDipolePtr> dl;
sl.extract(back_inserter(dl));
vector<tDipolePtr> dr;
sr.extract(back_inserter(dr));
double sum = 0.0;
for ( int i = 0, N = dl.size(); i < N; ++i )
for ( int j = 0, M = dr.size(); j < M; ++j )
sum += fij(*dl[i], *dr[j], b);
return sum;
}
double DipoleXSec::
sumf(const ImpactParameters & b, const DipoleState & sl, const DipoleState & sr) const {
Nfij = 0;
NBVeto = 0;
NScalVeto = 0;
scalVeto = 0.0;
bVeto = 0.0;
vector<tDipolePtr> dl;
sl.extract(back_inserter(dl));
vector<tDipolePtr> dr;
sr.extract(back_inserter(dr));
double sum = 0.0;
for ( int i = 0, N = dl.size(); i < N; ++i )
for ( int j = 0, M = dr.size(); j < M; ++j ) {
DipoleInteraction di = fij(b, *dl[i], *dr[j]);
sum += di.f2/2.0; /*** TODO: this is because we need proper f here. FIX CONFUSION */
}
return sum;
}
DipoleXSec::FList
DipoleXSec::flist(const DipoleState & sl, const DipoleState & sr,
const ImpactParameters & b) const {
FList ret;
vector<tDipolePtr> dl;
sl.extract(back_inserter(dl));
vector<tDipolePtr> dr;
sr.extract(back_inserter(dr));
//dont save the lowest fij. Otherwise the FList for LHC PbPb will take too much memory.
double cutoff = min(0.00000001, 1.0/double(dl.size()*dr.size()));
if ( dl.size()*dr.size() < 1000000 ) cutoff = 0.0000000001;
int total = 0;
int skipped = 0;
for ( int i = 0, N = dl.size(); i < N; ++i ) {
for ( int j = 0, M = dr.size(); j < M; ++j ) {
double f = fij(dl[i]->partons(), dr[j]->partons(), b, false)*2.0;
//extra 2 added from the non-diffractive interaction probability.
total++;
if ( f > cutoff &&
!kinematicsVeto(dl[i]->partons(), dr[j]->partons(), b) )
ret.insert(make_pair(make_pair(f, unitarize(f)), make_pair(dl[i], dr[j])));
else skipped++;
}
}
return ret;
}
DipoleInteraction::List
DipoleXSec::flist(const ImpactParameters & b,
const DipoleState & sl, const DipoleState & sr) const {
DipoleInteraction::List ret;
vector<tDipolePtr> dl;
sl.extract(back_inserter(dl));
vector<tDipolePtr> dr;
sr.extract(back_inserter(dr));
//dont save the lowest fij. Otherwise the FList for LHC PbPb will take too much memory.
double cutoff = min(0.00000001, 1.0/double(dl.size()*dr.size()));
if ( dl.size()*dr.size() < 1000000 ) cutoff = 0.0000000001;
int total = 0;
int skipped = 0;
for ( int i = 0, N = dl.size(); i < N; ++i ) {
for ( int j = 0, M = dr.size(); j < M; ++j ) {
/*** TODO: WHY IS VETO FALSE - because we don't want to check
kinematics if below cut ***/
DipoleInteraction di = fij(b, *dl[i], *dr[j], false);
total++;
if ( di.f2 > cutoff && !kinematicsVeto(di) )
ret.insert(di);
else
skipped++;
}
}
return ret;
}
DipoleXSec::InteractionRecoil
DipoleXSec::recoil(const DipoleInteraction & di) const {
InteractionRecoil ret;
if ( di.norec ) return ret;
if ( di.ints.first == di.dips.first->partons().first )
ret.first.first = di.rec;
else
ret.first.second = di.rec;
if ( di.ints.second == di.dips.second->partons().first )
ret.second.first = di.b->invRotatePT(-di.rec);
else
ret.second.second = di.b->invRotatePT(-di.rec);
return ret;
}
DipoleXSec::InteractionRecoil
DipoleXSec::recoil(const pair<tPartonPtr, tPartonPtr> left,
const pair<tPartonPtr, tPartonPtr> right,
const ImpactParameters & b,
pair<pair<bool, bool>, pair<bool, bool> > doesInt) const {
return recoil(left, right, b, doesInt, int0Partons(left.first, left.second,
right.first, right.second, b));
}
DipoleXSec::InteractionRecoil
DipoleXSec::recoil(const pair<tPartonPtr, tPartonPtr> left,
const pair<tPartonPtr, tPartonPtr> right,
const ImpactParameters & b,
pair<pair<bool, bool>, pair<bool, bool> > doesInt,
pair<bool, bool> ints) const {
tPartonPtr p1 = left.first;
tPartonPtr p2 = left.second;
tPartonPtr p3 = right.first;
tPartonPtr p4 = right.second;
double pTScale = Current<DipoleEventHandler>()->emitter().pTScale();
if ( theInteraction == 2 ) { //swing recoil, only new dips
TransverseMomentum rec14 = -pTScale*b.difference(p1->position(), p4->position())/
( b.dist2(*p1,*p4) );
TransverseMomentum rec23 = -pTScale*b.difference(p2->position(), p3->position())/
( b.dist2(*p2,*p3) );
return make_pair(make_pair(rec14, rec23),
make_pair(-rec23, -rec14));
}
//4 parton-parton recoils (full diagonals)
if ( theInteraction == 1 ) {
TransverseMomentum rec13 = -pTScale*b.difference(p1->position(), p3->position())/
(b.dist2(*p1,*p3) );
TransverseMomentum rec14 = -pTScale*b.difference(p1->position(), p4->position())/
(b.dist2(*p1,*p4) );
TransverseMomentum rec23 = -pTScale*b.difference(p2->position(), p3->position())/
(b.dist2(*p2,*p3) );
TransverseMomentum rec24 = -pTScale*b.difference(p2->position(), p4->position())/
(b.dist2(*p2,*p4) );
return make_pair(make_pair(rec14 + rec13, rec23 + rec24),
make_pair(b.invRotatePT(-rec23 - rec13),
b.invRotatePT(-rec14 - rec24)));
}
//2 parton-parton recoils (no diagonals)
if ( theInteraction == 3 ) {
TransverseMomentum rec13 = -pTScale*b.difference(p1->position(), p3->position())/
(b.dist2(*p1,*p3) );
TransverseMomentum rec14 = -pTScale*b.difference(p1->position(), p4->position())/
(b.dist2(*p1,*p4) );
TransverseMomentum rec23 = -pTScale*b.difference(p2->position(), p3->position())/
(b.dist2(*p2,*p3) );
TransverseMomentum rec24 = -pTScale*b.difference(p2->position(), p4->position())/
(b.dist2(*p2,*p4) );
if ( doesInt.first.first && doesInt.first.second &&
doesInt.second.first && doesInt.second.second )
return make_pair(make_pair(rec14, rec23),
make_pair(b.invRotatePT(-rec23), b.invRotatePT(-rec14)));
else {
TransverseMomentum rec11 = TransverseMomentum();
if ( doesInt.first.first && doesInt.second.first ) rec11 += rec13;
if ( doesInt.first.first && doesInt.second.second ) rec11 += rec14;
TransverseMomentum rec12 = TransverseMomentum();
if ( doesInt.first.second && doesInt.second.first ) rec12 += rec23;
if ( doesInt.first.second && doesInt.second.second ) rec12 += rec24;
TransverseMomentum rec21 = TransverseMomentum();
if ( doesInt.second.first && doesInt.first.first ) rec21 -= rec13;
if ( doesInt.second.first && doesInt.first.second ) rec21 -= rec23;
TransverseMomentum rec22 = TransverseMomentum();
if ( doesInt.second.second && doesInt.first.first ) rec22 -= rec14;
if ( doesInt.second.second && doesInt.first.second ) rec22 -= rec24;
return make_pair(make_pair(rec11, rec12),
make_pair(b.invRotatePT(rec21), b.invRotatePT(rec22)));
}
}
if ( theInteraction == 0 ) { //dip-dip recoil
TransverseMomentum rec13 = -pTScale*b.difference(p1->position(), p3->position())/
(b.dist2(*p1,*p3) );
TransverseMomentum rec14 = -pTScale*b.difference(p1->position(), p4->position())/
(b.dist2(*p1,*p4) );
TransverseMomentum rec23 = -pTScale*b.difference(p2->position(), p3->position())/
(b.dist2(*p2,*p3) );
TransverseMomentum rec24 = -pTScale*b.difference(p2->position(), p4->position())/
(b.dist2(*p2,*p4) );
/*** TODO: WHAT IS THIS? ***/
if ( p1->oY() + p3->oY() > p2->oY() + p4->oY() ) rec24 = TransverseMomentum();
else rec13 = TransverseMomentum();
// if ( p1->flavour() == ParticleID::g ) {
// rec13 /= 2.0;
// rec14 /= 2.0;
// }
// if ( p2->flavour() == ParticleID::g ) {
// rec23 /= 2.0;
// rec24 /= 2.0;
// }
// if ( p3->flavour() == ParticleID::g ) {
// rec13 /= 2.0;
// rec23 /= 2.0;
// }
// if ( p4->flavour() == ParticleID::g ) {
// rec14 /= 2.0;
// rec24 /= 2.0;
// }
TransverseMomentum rec11 = TransverseMomentum();
if ( doesInt.first.first && doesInt.second.first ) rec11 += rec13;
if ( doesInt.first.first && doesInt.second.second ) rec11 += rec14;
TransverseMomentum rec12 = TransverseMomentum();
if ( doesInt.first.second && doesInt.second.first ) rec12 += rec23;
if ( doesInt.first.second && doesInt.second.second ) rec12 += rec24;
TransverseMomentum rec21 = TransverseMomentum();
if ( doesInt.second.first && doesInt.first.first ) rec21 -= rec13;
if ( doesInt.second.first && doesInt.first.second ) rec21 -= rec23;
TransverseMomentum rec22 = TransverseMomentum();
if ( doesInt.second.second && doesInt.first.first ) rec22 -= rec14;
if ( doesInt.second.second && doesInt.first.second ) rec22 -= rec24;
return make_pair(make_pair(rec11, rec12),
make_pair(b.invRotatePT(rec21), b.invRotatePT(rec22)));
}
return make_pair(make_pair(TransverseMomentum(), TransverseMomentum()),
make_pair(TransverseMomentum(), TransverseMomentum()));
}
DipoleXSec::RealInteraction
DipoleXSec::initialiseInteraction(const pair<DipolePtr, DipolePtr> inter,
RealPartonStatePtr lrs, RealPartonStatePtr rrs,
pair<pair<bool, bool>, pair<bool, bool> > doesInt,
const ImpactParameters & b) const {
RealInteraction ret;
ret.lrs = lrs;
ret.rrs = rrs;
ret.d1 = inter.first;
ret.d2 = inter.second;
if ( doesInt.first.first ) ret.p11 = lrs->getReal(ret.d1->partons().first);
else ret.p11 = RealPartonPtr();
if ( doesInt.first.second ) ret.p12 = lrs->getReal(ret.d1->partons().second);
else ret.p12 = RealPartonPtr();
if ( doesInt.second.first ) ret.p21 = rrs->getReal(inter.second->partons().first);
else ret.p21 = RealPartonPtr();
if ( doesInt.second.second ) ret.p22 = rrs->getReal(inter.second->partons().second);
else ret.p22 = RealPartonPtr();
if ( ret.p11 && ret.p11->fluct != -1 && ret.p12 && ret.p12->fluct == ret.p11->fluct )
lrs->splitFluct(ret.p11, ret.p12);
if ( ret.p21 && ret.p21->fluct != -1 && ret.p22 && ret.p22->fluct == ret.p21->fluct )
rrs->splitFluct(ret.p21, ret.p22);
ret.range11 = sqrt(min(min(inter.first->partons().first->dist2
(*inter.second->partons().first),
inter.first->partons().first->dist2
(*inter.second->partons().second)),
inter.first->partons().first->dist2
(*inter.first->partons().second)/4.0));
ret.range12 = sqrt(min(min(inter.first->partons().second->dist2
(*inter.second->partons().first),
inter.first->partons().second->dist2
(*inter.second->partons().second)),
inter.first->partons().second->dist2
(*inter.first->partons().first)/4.0));
ret.range21 = sqrt(min(min(inter.second->partons().first->dist2
(*inter.first->partons().first),
inter.second->partons().first->dist2
(*inter.first->partons().second)),
inter.second->partons().first->dist2
(*inter.second->partons().second)/4.0));
ret.range22 = sqrt(min(min(inter.second->partons().second->dist2
(*inter.first->partons().first),
inter.second->partons().second->dist2
(*inter.first->partons().second)),
inter.second->partons().second->dist2
(*inter.second->partons().first)/4.0));
if ( theInteraction == 3 ) {
ret.range11 = sqrt(min(inter.first->partons().first->dist2
(*inter.second->partons().second),
inter.first->partons().first->dist2
(*inter.first->partons().second)/4.0));
ret.range12 = sqrt(min(inter.first->partons().second->dist2
(*inter.second->partons().first),
inter.first->partons().second->dist2
(*inter.first->partons().first)/4.0));
ret.range21 = sqrt(min(inter.second->partons().first->dist2
(*inter.first->partons().second),
inter.second->partons().first->dist2
(*inter.second->partons().second)/4.0));
ret.range22 = sqrt(min(inter.second->partons().second->dist2
(*inter.first->partons().first),
inter.second->partons().second->dist2
(*inter.second->partons().first)/4.0));
}
if ( theInteraction == 0 ) {
pair<bool, bool> ints;
ints = int0Partons(ret.d1->partons().first, ret.d1->partons().second,
inter.second->partons().first, inter.second->partons().second, b);
InvEnergy2 range2;
if ( ints.first && ints.second ) range2 = b.dist2(*ret.p11->theParton,*ret.p21->theParton);
if ( !ints.first && ints.second ) range2 = b.dist2(*ret.p12->theParton,*ret.p21->theParton);
if ( ints.first && !ints.second ) range2 = b.dist2(*ret.p11->theParton,*ret.p22->theParton);
if ( !ints.first && !ints.second ) range2 = b.dist2(*ret.p12->theParton,*ret.p22->theParton);
InvEnergy range = sqrt(range2);
ret.range11 = min(range, ret.d1->size()/2.0);
ret.range12 = min(range, ret.d1->size()/2.0);
ret.range21 = min(range, ret.d2->size()/2.0);
ret.range22 = min(range, ret.d2->size()/2.0);
if ( Current<DipoleEventHandler>()->emitter().rangeMode() == 1 ) {
ret.range11 = ret.d1->size()/2.0;
ret.range12 = ret.d1->size()/2.0;
ret.range21 = ret.d2->size()/2.0;
ret.range22 = ret.d2->size()/2.0;
}
}
InvEnergy2 rr11 = (ret.p11 && ret.p21) ? ret.p11->theParton->dist2(*ret.p21->theParton):ZERO;
InvEnergy2 rr12 = (ret.p11 && ret.p22) ? ret.p11->theParton->dist2(*ret.p22->theParton):ZERO;
InvEnergy2 rr21 = (ret.p12 && ret.p21) ? ret.p12->theParton->dist2(*ret.p21->theParton):ZERO;
InvEnergy2 rr22 = (ret.p12 && ret.p22) ? ret.p12->theParton->dist2(*ret.p22->theParton):ZERO;
InvEnergy2 rm11 = (ret.p11 && ret.p11->mothers.first ?
ret.p11->theParton->dist2(*ret.p11->mothers.first->theParton):ZERO);
InvEnergy2 rm12 = (ret.p12 && ret.p12->mothers.first ?
ret.p12->theParton->dist2(*ret.p12->mothers.first->theParton):ZERO);
InvEnergy2 rm21 = (ret.p21 && ret.p21->mothers.first ?
ret.p21->theParton->dist2(*ret.p21->mothers.first->theParton):ZERO);
InvEnergy2 rm22 = (ret.p22 && ret.p22->mothers.first ?
ret.p22->theParton->dist2(*ret.p22->mothers.first->theParton):ZERO);
ret.max11 = rr11 < rm11 && rr11 < rm21;
ret.max12 = rr12 < rm11 && rr12 < rm22;
ret.max21 = rr21 < rm12 && rr21 < rm21;
ret.max22 = rr22 < rm12 && rr22 < rm22;
Energy2 den = ((ret.p11 ? 1./rr11:ZERO) + (ret.p12 ? 1./rr12:ZERO)
+ (ret.p21 ? 1./rr21:ZERO) + (ret.p22 ? 1./rr22:ZERO));
ret.P11 = ret.p11 ? (1./rr11)/den:ZERO;
ret.P12 = ret.p12 ? (1./rr12)/den:ZERO;
ret.P21 = ret.p21 ? (1./rr21)/den:ZERO;
ret.P22 = ret.p22 ? (1./rr22)/den:ZERO;
//look only at the new dipole pairs, the cross pairs get zero weight.
if ( theInteraction == 3 ) {
den = (1./rr12 + 1./rr21);
ret.P12 = (1./rr12)/den;
ret.P21 = (1./rr21)/den;
ret.P11 = 0.0;
ret.P22 = 0.0;
}
//set only the selected pair to weight 1, the others to 0.
if ( theInteraction == 0 ) {
pair<bool, bool> ints;
ints = int0Partons(inter.first->partons().first, inter.first->partons().second,
inter.second->partons().first, inter.second->partons().second, b);
ret.P11 = 0.0;
ret.P22 = 0.0;
ret.P21 = 0.0;
ret.P12 = 0.0;
if ( ints.first && ints.second ) ret.P11 = 1.0;
if ( ints.first && !ints.second ) ret.P12 = 1.0;
if ( !ints.first && ints.second ) ret.P21 = 1.0;
if ( !ints.first && !ints.second ) ret.P22 = 1.0;
}
return ret;
}
DipoleXSec::RealInteraction
DipoleXSec::initialiseInteraction(const DipoleInteraction & di,
RealPartonStatePtr lrs, RealPartonStatePtr rrs) const {
RealInteraction ret;
const ImpactParameters & b = *di.b;
ret.lrs = lrs;
ret.rrs = rrs;
ret.d1 = di.dips.first;
ret.d2 = di.dips.second;
if ( ret.d1->partons().first == di.ints.first ) ret.p11 = lrs->getReal(di.ints.first);
if ( ret.d1->partons().second == di.ints.first ) ret.p12 = lrs->getReal(di.ints.first);
if ( ret.d2->partons().first == di.ints.second ) ret.p21 = rrs->getReal(di.ints.second);
if ( ret.d2->partons().second == di.ints.second ) ret.p22 = rrs->getReal(di.ints.second);
if ( ret.p11 && ret.p11->fluct != -1 && ret.p12 && ret.p12->fluct == ret.p11->fluct )
lrs->splitFluct(ret.p11, ret.p12);
if ( ret.p21 && ret.p21->fluct != -1 && ret.p22 && ret.p22->fluct == ret.p21->fluct )
rrs->splitFluct(ret.p21, ret.p22);
ret.range11 = sqrt(min(min(ret.d1->partons().first->dist2
(*ret.d2->partons().first),
ret.d1->partons().first->dist2
(*ret.d2->partons().second)),
ret.d1->partons().first->dist2
(*ret.d1->partons().second)/4.0));
ret.range12 = sqrt(min(min(ret.d1->partons().second->dist2
(*ret.d2->partons().first),
ret.d1->partons().second->dist2
(*ret.d2->partons().second)),
ret.d1->partons().second->dist2
(*ret.d1->partons().first)/4.0));
ret.range21 = sqrt(min(min(ret.d2->partons().first->dist2
(*ret.d1->partons().first),
ret.d2->partons().first->dist2
(*ret.d1->partons().second)),
ret.d2->partons().first->dist2
(*ret.d2->partons().second)/4.0));
ret.range22 = sqrt(min(min(ret.d2->partons().second->dist2
(*ret.d1->partons().first),
ret.d2->partons().second->dist2
(*ret.d1->partons().second)),
ret.d2->partons().second->dist2
(*ret.d2->partons().first)/4.0));
pair<bool, bool> ints(ret.d1->partons().first == di.ints.first,
ret.d2->partons().first == di.ints.second);
InvEnergy2 range2 = b.dist2(*di.ints.first, *di.ints.second);
InvEnergy range = sqrt(range2);
ret.range11 = min(range, ret.d1->size()/2.0);
ret.range12 = min(range, ret.d1->size()/2.0);
ret.range21 = min(range, ret.d2->size()/2.0);
ret.range22 = min(range, ret.d2->size()/2.0);
if ( Current<DipoleEventHandler>()->emitter().rangeMode() == 1 ) {
ret.range11 = ret.d1->size()/2.0;
ret.range12 = ret.d1->size()/2.0;
ret.range21 = ret.d2->size()/2.0;
ret.range22 = ret.d2->size()/2.0;
}
InvEnergy2 rr11 = (ret.p11 && ret.p21) ? ret.p11->theParton->dist2(*ret.p21->theParton):ZERO;
InvEnergy2 rr12 = (ret.p11 && ret.p22) ? ret.p11->theParton->dist2(*ret.p22->theParton):ZERO;
InvEnergy2 rr21 = (ret.p12 && ret.p21) ? ret.p12->theParton->dist2(*ret.p21->theParton):ZERO;
InvEnergy2 rr22 = (ret.p12 && ret.p22) ? ret.p12->theParton->dist2(*ret.p22->theParton):ZERO;
InvEnergy2 rm11 = (ret.p11 && ret.p11->mothers.first ?
ret.p11->theParton->dist2(*ret.p11->mothers.first->theParton):ZERO);
InvEnergy2 rm12 = (ret.p12 && ret.p12->mothers.first ?
ret.p12->theParton->dist2(*ret.p12->mothers.first->theParton):ZERO);
InvEnergy2 rm21 = (ret.p21 && ret.p21->mothers.first ?
ret.p21->theParton->dist2(*ret.p21->mothers.first->theParton):ZERO);
InvEnergy2 rm22 = (ret.p22 && ret.p22->mothers.first ?
ret.p22->theParton->dist2(*ret.p22->mothers.first->theParton):ZERO);
ret.max11 = rr11 < rm11 && rr11 < rm21;
ret.max12 = rr12 < rm11 && rr12 < rm22;
ret.max21 = rr21 < rm12 && rr21 < rm21;
ret.max22 = rr22 < rm12 && rr22 < rm22;
Energy2 den = ((ret.p11 ? 1./rr11:ZERO) + (ret.p12 ? 1./rr12:ZERO)
+ (ret.p21 ? 1./rr21:ZERO) + (ret.p22 ? 1./rr22:ZERO));
ret.P11 = ret.p11 ? (1./rr11)/den:ZERO;
ret.P12 = ret.p12 ? (1./rr12)/den:ZERO;
ret.P21 = ret.p21 ? (1./rr21)/den:ZERO;
ret.P22 = ret.p22 ? (1./rr22)/den:ZERO;
//look only at the new dipole pairs, the cross pairs get zero weight.
if ( theInteraction == 3 ) {
den = (1./rr12 + 1./rr21);
ret.P12 = (1./rr12)/den;
ret.P21 = (1./rr21)/den;
ret.P11 = 0.0;
ret.P22 = 0.0;
}
//set only the selected pair to weight 1, the others to 0.
ret.P11 = 0.0;
ret.P22 = 0.0;
ret.P21 = 0.0;
ret.P12 = 0.0;
if ( ints.first && ints.second ) ret.P11 = 1.0;
if ( ints.first && !ints.second ) ret.P12 = 1.0;
if ( !ints.first && ints.second ) ret.P21 = 1.0;
if ( !ints.first && !ints.second ) ret.P22 = 1.0;
return ret;
}
pair<pair<bool, bool>, pair<bool, bool> >
DipoleXSec::doesInt(const pair<tPartonPtr, tPartonPtr> left,
const pair<tPartonPtr, tPartonPtr> right,
const ImpactParameters & b) const {
//decide which of the 4 partons will end up on shell
return make_pair(make_pair(true, true), make_pair(true, true));
if ( Current<DipoleEventHandler>()->eventFiller().singleMother() ) {
//when single mother not all are chosen
pair<pair<bool, bool>, pair<bool, bool> > ret =
make_pair(make_pair(false, false), make_pair(false, false));
if ( theInteraction == 0 ) {
//if interaction 0, then always keep the same partons used in fij
pair<bool, bool> ints;
ints = int0Partons(left.first, left.second, right.first, right.second, b);
if ( ints.first ) ret.first.first = true;
else ret.first.second = true;
if ( ints.second ) ret.second.first =true;
else ret.second.second = true;
}
else {
//select which two partons the gluons are exchanged between,
//weighted by 1/distance^2
InvEnergy2 r11 = b.dist2(*left.first, *right.first);
InvEnergy2 r12 = b.dist2(*left.first, *right.second);
InvEnergy2 r21 = b.dist2(*left.second, *right.first);
InvEnergy2 r22 = b.dist2(*left.second, *right.second);
Selector<int, double> sel;
sel.insert(1./r11/(1./r11 + 1./r12 + 1./r21 + 1./r22), 1);
sel.insert(1./r12/(1./r11 + 1./r12 + 1./r21 + 1./r22), 2);
sel.insert(1./r21/(1./r11 + 1./r12 + 1./r21 + 1./r22), 3);
sel.insert(1./r22/(1./r11 + 1./r12 + 1./r21 + 1./r22), 4);
double dummy = (r11 + r22 + r12 + r21)*sqr(GeV);
double pseudoRnd = 10000.0*dummy - int(10000.0*dummy);
switch ( sel[pseudoRnd] ) {
case 1:
ret.first.first = true;
ret.second.first = true;
break;
case 2:
ret.first.first = true;
ret.second.second = true;
break;
case 3:
ret.first.second = true;
ret.second.first = true;
break;
case 4:
ret.first.second = true;
ret.second.second = true;
break;
}
}
//check for swinged emission, in which case both partons are real
if ( !(left.first->parents().second == left.second ||
left.second->parents().first == left.first) ) {
ret.first.first = true;
ret.first.second = true;
}
if ( !(right.first->parents().second == right.second ||
right.second->parents().first == right.first) ) {
ret.second.first = true;
ret.second.second = true;
}
return ret;
}
//if not single mother, all partons are on shell
return make_pair(make_pair(true, true), make_pair(true, true));
}
pair<bool, bool> DipoleXSec::int0Partons(tcPartonPtr p11, tcPartonPtr p12,
tcPartonPtr p21, tcPartonPtr p22,
const ImpactParameters & b) const {
//chose partons (pseudo-)randomly
InvEnergy2 rr11 = b.dist2(*p11, *p21);
InvEnergy2 rr22 = b.dist2(*p12, *p22);
double dummy = (rr11 + rr22)*GeV2;
while ( dummy > 100 ) dummy /= 10;
bool firstDipole = (1000000.0*dummy - long(1000000.0*dummy) > 0.5 );
bool secondDipole = (10000000.0*dummy - long(10000000.0*dummy) > 0.5 );
return make_pair(firstDipole, secondDipole);
}
void DipoleXSec::updateMomenta(RealInteraction * i) const {
if ( i->p11 ) {
pair<Energy, Energy> pm11 = i->p11->effectivePlusMinus(i->range11, true);
i->effectivePlus11 = pm11.first;
i->effectiveMinus11 = pm11.second;
}
if ( i->p12 ) {
pair<Energy, Energy> pm12 = i->p12->effectivePlusMinus(i->range12, false);
i->effectivePlus12 = pm12.first;
i->effectiveMinus12 = pm12.second;
}
if ( i->p21 ) {
pair<Energy, Energy> pm21 = i->p21->effectivePlusMinus(i->range21, true);
i->effectivePlus21 = pm21.first;
i->effectiveMinus21 = pm21.second;
}
if ( i->p22 ) {
pair<Energy, Energy> pm22 = i->p22->effectivePlusMinus(i->range22, false);
i->effectivePlus22 = pm22.first;
i->effectiveMinus22 = pm22.second;
}
}
void DipoleXSec::doTransverseRecoils(RealInteraction i, InteractionRecoil recs) const {
if ( Current<DipoleEventHandler>()->eventFiller().mode() != 1 ) {
i.lrs->totalRecoil += recs.first.first + recs.first.second;
i.rrs->totalRecoil += recs.second.first + recs.second.second;
if ( i.p11 ) i.p11->intRecoil += recs.first.first;
if ( i.p12 ) i.p12->intRecoil += recs.first.second;
if ( i.p21 ) i.p21->intRecoil += recs.second.first;
if ( i.p22 ) i.p22->intRecoil += recs.second.second;
}
if ( i.p11 ) {
if ( i.p22 ) i.p22->doEffectiveRecoil(i.p11, i.range11, true, ZERO, -recs.first.first);
else i.p21->doEffectiveRecoil(i.p11, i.range11, true, ZERO, -recs.first.first);
}
if ( i.p12 ) {
if ( i.p22 ) i.p22->doEffectiveRecoil(i.p12, i.range12, false, ZERO, -recs.first.second);
else i.p21->doEffectiveRecoil(i.p12, i.range12, false, ZERO, -recs.first.second);
}
if ( i.p21 ) {
if ( i.p12 ) i.p12->doEffectiveRecoil(i.p21, i.range21, true, ZERO, -recs.second.first);
else i.p11->doEffectiveRecoil(i.p21, i.range21, true, ZERO, -recs.second.first);
}
if ( i.p22 ) {
if ( i.p11 ) i.p11->doEffectiveRecoil(i.p22, i.range22, false, ZERO, -recs.second.second);
else i.p12->doEffectiveRecoil(i.p22, i.range22, false, ZERO, -recs.second.second);
}
}
pair<double, double> DipoleXSec::findBoosts(RealInteraction i) const {
Energy neededMinus = ZERO;
if ( checkOffShell ) neededMinus += i.lrs->neededValenceMinus();
if ( i.p11 ) neededMinus += i.p11->effectiveGiveMinus(i.range11, true);
if ( i.p12 ) neededMinus += i.p12->effectiveGiveMinus(i.range12, false);
Energy neededPlus = ZERO;
if ( checkOffShell ) neededPlus += i.rrs->neededValenceMinus();
if ( i.p21 ) neededPlus += i.p21->effectiveGiveMinus(i.range21, true);
if ( i.p22 ) neededPlus += i.p22->effectiveGiveMinus(i.range22, false);
Energy intPlus1 = ZERO;
if ( i.p11 ) intPlus1 += i.effectivePlus11;
if ( i.p12 ) intPlus1 += i.effectivePlus12;
Energy intPlus2 = ZERO;
if ( i.p21 ) intPlus2 += i.p21->inRangeMinus(i.range21, true);
if ( i.p22 ) intPlus2 += i.p22->inRangeMinus(i.range22, false);
Energy intMinus1 = ZERO;
if ( i.p11 ) intMinus1 += i.p11->inRangeMinus(i.range11, true);
if ( i.p12 ) intMinus1 += i.p12->inRangeMinus(i.range12, false);
Energy intMinus2 = ZERO;
if ( i.p21 ) intMinus2 += i.effectivePlus21;
if ( i.p22 ) intMinus2 += i.effectivePlus22;
Energy evoPlus2 = neededPlus - intPlus2;
Energy evoMinus1 = neededMinus - intMinus1;
pair<double, double> boosts = findBoosts(intPlus1, intPlus2, intMinus1, intMinus2, evoPlus2, evoMinus1);
if ( boosts.first == 0.0 || isnan(boosts.first) || isnan(boosts.second) ) {
return make_pair(0.0,0.0);
}
if ( boosts.second < 0.0 || boosts.first < 0.0 ) {
//can be over 1.0 if the valenceminus is more than needed
return make_pair(0.0,0.0);;
}
return boosts;
}
void DipoleXSec::doBoosts(RealInteraction i, pair<double, double> boosts) const {
doBoost(i.p11, i.range11, i.p12, i.range12, boosts.first);
doBoost(i.p21, i.range21, i.p22, i.range22, boosts.second);
}
void DipoleXSec::reduceRecoil(RealInteraction RI, InteractionRecoil recs) const {
//check if the interacting partons have passed each other
while ( max((RI.p11 ? RI.p11->y:RI.p12->y), (RI.p12 ? RI.p12->y:RI.p11->y)) >
-max((RI.p21 ? RI.p21->y:RI.p22->y), (RI.p22 ? RI.p22->y:RI.p21->y)) ) {
//and if the recoil is larger than 1 GeV
if ( max(max(recs.first.first.pt(), recs.first.second.pt()),
max(recs.second.first.pt(), recs.second.second.pt())) < 1*GeV ) {
break;
}
//then undo the recoil and redo 90% of it.
recs = make_pair(make_pair(-recs.first.first, -recs.first.second),
make_pair(-recs.second.first, -recs.second.second));
doTransverseRecoils(RI, recs);
recs = make_pair(make_pair(-recs.first.first*0.9, -recs.first.second*0.9),
make_pair(-recs.second.first*0.9, -recs.second.second*0.9));
doTransverseRecoils(RI, recs);
}
}
bool DipoleXSec::doInteraction(InteractionRecoil recs, const FList::const_iterator inter,
RealPartonStatePtr lrs, RealPartonStatePtr rrs,
pair<pair<bool, bool>, pair<bool, bool> > doesInt,
const ImpactParameters & b) const {
static DebugItem checkkinematics("DIPSY::CheckKinematics", 6);
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found at start of interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found at start of interaction!" << Exception::warning;
RealInteraction RI = initialiseInteraction(inter->second, lrs, rrs, doesInt, b);
doTransverseRecoils(RI, recs);
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after transverse recoil in interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after transverse recoil in interaction!" << Exception::warning;
if ( theRecoilReduction ) {
reduceRecoil(RI, recs);
}
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after recoil reduction in interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after recoil reduction in interaction!" << Exception::warning;
updateMomenta(& RI);
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after update in interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after update in interaction!" << Exception::warning;
pair<double, double> boosts = findBoosts(RI);
if ( boosts.first == 0.0 ) {
return false;
}
doBoosts(RI, boosts);
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after boost in interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after boost in interaction!" << Exception::warning;
if ( ordered(RI, recs, b) ) {
setOnShell(RI);
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found at end of interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found at end of interaction!" << Exception::warning;
return true;
}
return false;
}
bool DipoleXSec::doInteraction(InteractionRecoil recs,
const DipoleInteraction::List::const_iterator inter,
RealPartonStatePtr lrs, RealPartonStatePtr rrs,
pair<pair<bool, bool>, pair<bool, bool> > doesInt) const {
const ImpactParameters & b = *(inter->b);
static DebugItem checkkinematics("DIPSY::CheckKinematics", 6);
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found at start of interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found at start of interaction!" << Exception::warning;
RealInteraction RI = initialiseInteraction(*inter, lrs, rrs);
doTransverseRecoils(RI, recs);
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after transverse recoil in interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after transverse recoil in interaction!" << Exception::warning;
if ( theRecoilReduction ) {
reduceRecoil(RI, recs);
}
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after recoil reduction in interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after recoil reduction in interaction!" << Exception::warning;
updateMomenta(& RI);
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after update in interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after update in interaction!" << Exception::warning;
pair<double, double> boosts = findBoosts(RI);
if ( boosts.first == 0.0 ) {
return false;
}
doBoosts(RI, boosts);
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after boost in interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found after boost in interaction!" << Exception::warning;
if ( ordered(RI, recs, b) ) {
setOnShell(RI);
if ( checkkinematics && lrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found at end of interaction!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<InteractionKinematicException>()
<< "negatives found at end of interaction!" << Exception::warning;
return true;
}
return false;
}
bool DipoleXSec::ordered(RealInteraction i, InteractionRecoil recs,
const ImpactParameters & b) const {
//this options means only momentum conservation, and that is already checked by
//finding a boost. no extra ordering wanted.
if ( theIntOrdering == 2 ) {
return true;
}
double PSInf = Current<DipoleEventHandler>()->emitter().PSInflation();
bool ordered = true;
bool both = Current<DipoleEventHandler>()->emitter().bothOrderedFS();
//if not a local pt max, check p+- ordering
if ( theIntOrdering == 0 ) {
if ( theInteraction == 0 ) {
//check ordering only for the key partons
pair<bool, bool> ints0 = int0Partons(i.d1->partons().first, i.d1->partons().second,
i.d2->partons().first, i.d2->partons().second, b);
Energy effMinus1 = (ints0.first ? i.effectiveMinus11:i.effectiveMinus12);
Energy effMinus2 = (ints0.second ? i.effectiveMinus21:i.effectiveMinus22);
Energy effPlus1 = (ints0.first ? i.effectivePlus11:i.effectivePlus12);
Energy effPlus2 = (ints0.second ? i.effectivePlus21:i.effectivePlus22);
if ( effPlus1*PSInf < effMinus2 ) ordered = false;
if ( effPlus2*PSInf < effMinus1 ) ordered = false;
}
else {
if ( !i.max11 && i.p11 && i.p21 ) {
if ( i.effectivePlus11*PSInf < i.effectiveMinus21*(both ? 1.0:i.P11) ) ordered = false;
if ( (both ? 1.0:i.P11)*i.effectiveMinus11/PSInf > i.effectivePlus21 ) ordered = false;
}
if ( !i.max12 && i.p11 && i.p22 ) {
if ( i.effectivePlus11*PSInf < i.effectiveMinus22*(both ? 1.0:i.P12) ) ordered = false;
if ( (both ? 1.0:i.P12)*i.effectiveMinus11/PSInf > i.effectivePlus22 ) ordered = false;
}
if ( !i.max21 && i.p12 && i.p21 ) {
if ( i.effectivePlus12*PSInf < i.effectiveMinus21*(both ? 1.0:i.P21) ) ordered = false;
if ( (both ? 1.0:i.P21)*i.effectiveMinus12/PSInf > i.effectivePlus21 ) ordered = false;
}
if ( !i.max22 && i.p12 && i.p22 ) {
if ( i.effectivePlus12*PSInf < i.effectiveMinus22*(both ? 1.0:i.P22) ) ordered = false;
if ( (both ? 1.0:i.P22)*i.effectiveMinus12/PSInf > i.effectivePlus22 ) ordered = false;
}
}
}
if ( theIntOrdering == 1 && theInteraction == 0 ) {
double PMOrd = Current<DipoleEventHandler>()->emitter().PMinusOrdering();
//check ordering only for the key partons
pair<bool, bool> ints0 = int0Partons(i.d1->partons().first, i.d1->partons().second,
i.d2->partons().first, i.d2->partons().second, b);
Energy effMinus1 = (ints0.first ? i.effectiveMinus11:i.effectiveMinus12);
Energy effMinus2 = (ints0.second ? i.effectiveMinus21:i.effectiveMinus22);
Energy rec1 = (ints0.first ? recs.first.first.pt():recs.first.second.pt());
Energy rec2 = (ints0.second ? recs.second.first.pt():recs.second.second.pt());
if ( sqr(max(rec1, rec2)*PSInf) < effMinus1*effMinus2*PMOrd ) ordered = false;
}
if ( (i.p11 && i.p11->searchNegative(i.range11, true)) ||
(i.p12 && i.p12->searchNegative(i.range12, false)) ||
(i.p21 && i.p21->searchNegative(i.range21, true)) ||
(i.p22 && i.p22->searchNegative(i.range22, false)) )
ordered = false;
return ordered;
}
void DipoleXSec::setOnShell(RealInteraction i) const {
if ( i.p11 ) i.p11->effectiveGiveMinus(i.range11, true);
if ( i.p12 ) i.p12->effectiveGiveMinus(i.range12, false);
if ( i.p21 ) i.p21->effectiveGiveMinus(i.range21, true);
if ( i.p22 ) i.p22->effectiveGiveMinus(i.range22, false);
}
bool DipoleXSec::reconnect(tDipolePtr d1, tDipolePtr d2) const {
if ( d1->children().second && !d1->children().first ) {
return reconnect(d1->children().second, d2);
}
if ( d2->children().second && !d2->children().first ) {
return reconnect(d1, d2->children().second);
}
if ( !d1->children().first && !d2->children().first ) { //none has rescattered
d1->swingDipole(d2);
Current<DipoleEventHandler>()->swinger().recombine(*d1);
interact(*d1->children().first, *d2->children().first);
// d1->children().first->interact(*d2->children().first);
// d2->children().first->interact(*d1->children().first);
return true;
}
tPartonPtr p11 = d1->partons().first;
tPartonPtr p12 = d1->partons().second;
tPartonPtr p21 = d2->partons().first;
tPartonPtr p22 = d2->partons().second;
tDipolePtr d11 = p11->dipoles().second;
tDipolePtr d12 = p12->dipoles().first;
tDipolePtr d21 = p21->dipoles().second;
tDipolePtr d22 = p22->dipoles().first;
tDipolePtr swing1;
tDipolePtr swing2;
if ( d11 == d12 ) d1 = d11; //dipole has swinged back
if ( d21 == d22 ) d2 = d21;
if ( !d1->children().first && !d2->children().first ) { //both are original dipole
swing1 = d1;
swing2 = d2;
}
else if ( d1->children().first && !d2->children().first ) { //one dip d1 is rescattering
tDipolePtr temp = d1;
d1 = d2;
d2 = temp;
p11 = d1->partons().first;
p12 = d1->partons().second;
p21 = d2->partons().first;
p22 = d2->partons().second;
d11 = p11->dipoles().second;
d12 = p12->dipoles().first;
d21 = p21->dipoles().second;
d22 = p22->dipoles().first;
}
if ( !d1->children().first && d2->children().first ) { //one dip d2 is rescattering
if ( p11 != d21->partons().second && p12 != d22->partons().first ) {
//no connections, swing with one of the two randomly
swing1 = d1;
if ( UseRandom::rnd() > 0.5 ) swing2 = d21;
else swing2 = d22;
}
else if ( p11 == d21->partons().second && p12 != d22->partons().first ) {
//share one swinged dip, swing with the other
swing1 = d1;
swing2 = d22;
}
else if ( p11 != d21->partons().second && p12 == d22->partons().first ) {
//share other swinged dip, swing with the first
swing1 = d1;
swing2 = d21;
}
else if ( p11 == d21->partons().second && p12 == d22->partons().first ) {
//share both swinged dip, swing back to original
swing1 = d21;
swing2 = d22;
}
}
else if ( d1->children().first && d2->children().first ) { //both dips are rescattering
bool found = false;
int i = 0;
while ( !found && i < 1000 ) {
if ( UseRandom::rnd() > 0.5 ) swing1 = d11;
else swing1 = d12;
if ( UseRandom::rnd() > 0.5 ) swing2 = d21;
else swing2 = d22;
if ( swing1 != swing2 && swing1->partons().first != swing2->partons().second &&
swing2->partons().first != swing1->partons().second )
found = true;
}
}
if ( !swing1 || !swing2 ) {
d1->dipoleState().diagnosis(true);
return false;
}
swing1->swingDipole(swing2);
Current<DipoleEventHandler>()->swinger().recombine(*swing1);
interact(*swing1->children().first, *swing2->children().first);
// swing1->children().first->interact(*swing2->children().first);
// swing2->children().first->interact(*swing1->children().first);
return true;
}
void DipoleXSec::interact(Dipole & d1, Dipole & d2) const {
d1.interact(d2, partonicInteraction());
d2.interact(d1, partonicInteraction());
}
vector<pair<DipolePtr, DipolePtr> >
DipoleXSec::getColourExchanges(tRealPartonStatePtr lrs, tRealPartonStatePtr rrs) const {
vector<pair<DipolePtr, DipolePtr> > ret;
while ( ret.empty() ) {
list<tDipolePtr>::iterator rightDip = rrs->interactions.begin();
for ( list<tDipolePtr>::iterator leftDip = lrs->interactions.begin();
leftDip != lrs->interactions.end(); leftDip++, rightDip++ ) {
if ( unitarize(fij((*leftDip)->partons(), (*rightDip)->partons(),
ImpactParameters(), false))*0.99
< UseRandom::rnd() || true )
ret.push_back(make_pair(*leftDip, *rightDip));
}
}
return ret;
}
pair< double, double> DipoleXSec::findBoosts(Energy intPlus1, Energy intPlus2,
Energy intMinus1, Energy intMinus2,
Energy evoPlus2, Energy evoMinus1) const {
Energy2 A = (intPlus1 - evoPlus2)*intMinus2;
Energy2 B = intPlus1*(evoMinus1 + intMinus1 - intMinus2) -
evoPlus2*(evoMinus1 - intMinus2) - intPlus2*intMinus2;
Energy2 C = -intPlus2*(evoMinus1 - intMinus2);
if ( sqr(B/(2*A)) - C/A < 0.0 ) return make_pair(0.0, 0.0);
double y = -B/(2*A) + sqrt(sqr(B/(2*A)) - C/A); //is the + sqrt() always the right solution?
double x = 1.0 - intPlus2/(y*intPlus1) - evoPlus2/intPlus1;
return make_pair(x, y);
}
void DipoleXSec::doBoost(tRealPartonPtr p1, InvEnergy range1,
tRealPartonPtr p2, InvEnergy range2, double x) const {
//we here have to use the plusweighted recoil to fit the x, y above.
//Other weights get a lot more complicated equations for the boosts.
if ( p1 ) {
pair<Energy, Energy> pm1 = p1->effectivePlusMinus(range1, true);
p1->doPlusWeightedRecoil(p1, range1, true, (1.0 - x)*pm1.first, TransverseMomentum());
}
if ( p2 ) {
pair<Energy, Energy> pm2 = p2->effectivePlusMinus(range2, false);
p2->doPlusWeightedRecoil(p2, range2, false, (1.0 - x)*pm2.first, TransverseMomentum());
}
}
double DipoleXSec::unitarize(double f) const {
return Math::exp1m(-f);
}
void DipoleXSec::persistentOutput(PersistentOStream & os) const {
os << ounit(theRMax, InvGeV) << theInteraction << sinFunction << usePartonicInteraction
<< theIntOrdering << theRecoilReduction << checkOffShell;
}
void DipoleXSec::persistentInput(PersistentIStream & is, int) {
is >> iunit(theRMax, InvGeV) >> theInteraction >> sinFunction >> usePartonicInteraction
>> theIntOrdering >> theRecoilReduction >> checkOffShell;
}
// Static variable needed for the type description system in ThePEG.
#include "ThePEG/Utilities/DescribeClass.h"
DescribeClass<DipoleXSec,HandlerBase>
describeDIPSYDipoleXSec("DIPSY::DipoleXSec", "libAriadne5.so libDIPSY.so");
void DipoleXSec::Init() {
static ClassDocumentation<DipoleXSec> documentation
("There is no documentation for the DipoleXSec class");
static Parameter<DipoleXSec,InvEnergy> interfaceRMax
("RMax",
"The confinement scale (in iverse GeV). If set to zero, "
"the value of <interface>DipoleEventHandler::RMax</interface> of the "
"controlling event handler will be used.",
&DipoleXSec::theRMax, InvGeV, 0.0*InvGeV, 0.0*InvGeV, 0*InvGeV,
true, false, Interface::lowerlim);
static Parameter<DipoleXSec, int> interfaceInteraction
("Interaction",
"Which interaction to be used. Determines f_{ij} and recoils. "
"0 is the sinus interaction with 4 identical recoils"
"1 is the sinus interaction with recoils between all 4 parton pairs"
"2 is the swing interaction"
"3 is the sinus interaction with recoils between the 2 parton pairs"
" that get the new dipoles",
&DipoleXSec::theInteraction, 1, 1, 0, 0,
true, false, Interface::lowerlim);
static Switch<DipoleXSec,int> interfaceSinFunction
("SinFunction",
"Determine what approximation to the sine functions in the interaction strength "
"to use.",
&DipoleXSec::sinFunction, 0, true, false);
static SwitchOption interfaceSinFunctionExact
(interfaceSinFunction,
"Exact",
"The actual function.",
0);
static SwitchOption interfaceSinFunctionAverage
(interfaceSinFunction,
"Average",
"Average over angle of dipole and of the resulting bessel fuction.",
1);
static Switch<DipoleXSec,int> interfaceIntOrdering
("IntOrdering",
"How the real state is found from the virtual cascade. Speed versus consistency.",
&DipoleXSec::theIntOrdering, 0, true, false);
static SwitchOption interfaceIntOrderingDefault
(interfaceIntOrdering,
"Default",
"plus and minus reuired to be ordered after all recoils.",
0);
static SwitchOption interfaceIntOrderingAsEvo
(interfaceIntOrdering,
"AsEvo",
"Try to emulate the ordering in the evolution by requesting ordering"
" on the colliding particle if it would've been part of the cascade."
" That is, ignore recoils from the colliding cascade.",
1);
static SwitchOption interfaceIntOrderingVeryOpen
(interfaceIntOrdering,
"VeryOpen",
"Only checks that there is enough energy to put the interaction recoil on shell."
" Does not care about ordering, or setting evolution pt on shell. Same as Emils code.",
2);
static Switch<DipoleXSec,int> interfaceRecoilReduction
("RecoilReduction",
"What to do with large recoils",
&DipoleXSec::theRecoilReduction, 0, true, false);
static SwitchOption interfaceRecoilReductionOff
(interfaceRecoilReduction,
"Off",
"Do nothing special.",
0);
static SwitchOption interfaceRecoilReductionRapidityOrdered
(interfaceRecoilReduction,
"RapidityOrdered",
"Reduce the recoil until all interacting partons are rapidity ordered with each other.",
1);
static Switch<DipoleXSec,bool> interfaceCheckOffShell
("CheckOffShell",
"Make sure there is energy available to put incoming particles on-shell. Only necessary for virtual photons.",
&DipoleXSec::checkOffShell, true, true, false);
static SwitchOption interfaceCheckOffShellYes
(interfaceCheckOffShell,
"Yes",
"Do the check",
true);
static SwitchOption interfaceCheckOffShellNo
(interfaceCheckOffShell,
"No",
"No check is made.",
false);
static Switch<DipoleXSec,bool> interfacePartonicInteraction
("PartonicInteraction",
"Flag determining if only one parton in each dipole is considered interacting or both.",
&DipoleXSec::usePartonicInteraction, false, true, false);
static SwitchOption interfacePartonicInteractionDipole
(interfacePartonicInteraction,
"Dipole",
"The both partons in a dipole interact.",
false);
static SwitchOption interfacePartonicInteractionParton
(interfacePartonicInteraction,
"Parton",
"Only one parton in a dipole interacts.",
true);
}
diff --git a/DIPSY/DipoleXSec.h b/DIPSY/DipoleXSec.h
--- a/DIPSY/DipoleXSec.h
+++ b/DIPSY/DipoleXSec.h
@@ -1,646 +1,651 @@
// -*- C++ -*-
#ifndef DIPSY_DipoleXSec_H
#define DIPSY_DipoleXSec_H
//
// This is the declaration of the DipoleXSec class.
//
#include "ThePEG/Handlers/HandlerBase.h"
#include "DipoleXSec.fh"
#include "Dipole.h"
#include "DipoleState.fh"
#include "ImpactParameters.h"
#include "RealPartonState.fh"
#include "RealParton.fh"
#include "EffectiveParton.h"
namespace DIPSY {
using namespace ThePEG;
/**
* A dipole-dipole interaction used for book keeping.
*/
struct DipoleInteraction {
/**
* Constructor taking two dipoles.
*/
DipoleInteraction(Dipole & dlin, Dipole & drin,
const ImpactParameters & bin);
/**
* Prepare an interactions to be checked for the first time.
*/
void prepare() const;
/**
* Check that the interaction can be performed.
*/
bool check(bool final = false) const;
/**
* Accept the interaction.
*/
void accept() const;
/**
* Reject this interaction.
*/
void reject() const;
/**
* Perform the interaction.
*/
bool perform() const {
return check(true);
}
/**
* The two dipoles.
*/
pair<tDipolePtr,tDipolePtr> dips;
/**
* The two dipoles next to these.
*/
pair<tDipolePtr,tDipolePtr> dnext;
/**
* The two dipoles previous to these.
*/
pair<tDipolePtr,tDipolePtr> dprev;
/**
* A pointer to the impact parameter object.
*/
const ImpactParameters * b;
/**
* Which partons are actually interacting?
*/
pair<tPartonPtr,tPartonPtr> ints;
/**
* Which partons are merely spectators
*/
pair<tPartonPtr,tPartonPtr> spec;
/**
* Which partons are actually interacting?
*/
mutable pair<tSPartonPtr,tSPartonPtr> sints;
/**
* The distance between the interacting partons.
*/
InvEnergy2 d2;
/**
* The interaction strength (times two for total xsec).
*/
double f2;
/**
* The unitarized interaction strength.
*/
double uf2;
/**
* The transverse momentum recoil associated with the interaction.
*/
TransverseMomentum rec;
/**
* Set to true if the two partons which interacts has already received a recoil.
*/
mutable bool norec;
/**
* The transverse momentum scale associated with the interactions.
*/
Energy kt;
+ /**
+ * The number of this interaction in the order of tried interactions.
+ */
+ mutable int id;
+
/** Struct for ordering by dipoles. */
struct DipoleOrder {
/** Main operator */
bool operator()(const DipoleInteraction & i1, const DipoleInteraction & i2) const {
return i1.dips.first < i2.dips.first ||
( i1.dips.first == i2.dips.first && i1.dips.second < i2.dips.second );
}
};
/** Struct for ordering by dipoles. */
struct StrengthOrder {
/** Main operator */
bool operator()(const DipoleInteraction & i1, const DipoleInteraction & i2) const {
return i1.f2 > i2.f2;
}
};
/**
* A settype used to store an ordered list of potential interactions
* between dipoles in colliding dipole systems. Ordered in
* decreasing interaction strength.
*/
typedef multiset<DipoleInteraction,DipoleInteraction::StrengthOrder> List;
/** Struct for ordering by scale. */
struct PTOrder {
/** Main operator */
bool operator()(const List::const_iterator & i1,
const List::const_iterator & i2) const {
return i1->kt > i2->kt;
}
};
/**
* A settype used to store an ordered list of potential interactions
* between dipoles in colliding dipole systems. Ordered in
* decreasing interaction strength.
*/
typedef multiset<List::const_iterator,PTOrder> PTSet;
/**
* A settype used to store a list of potential interactions
* between dipoles in colliding dipole systems.
*/
typedef set<DipoleInteraction,DipoleInteraction::DipoleOrder> DipList;
};
/**
* DipoleXSec is the base class of all objects capable of calculating
* the scattering probability for two colliding dipoles. The main
* virtual functions to be overridden are called fij, sumf and flist.
*
* @see \ref DipoleXSecInterfaces "The interfaces"
* defined for DipoleXSec.
*/
class DipoleXSec: public HandlerBase {
public:
/**
* A maptype used to store an ordered list of scattering
* probabilities (the bare one paired with the unitarized one)
* between dipoles in colliding dipole systems.
*/
typedef multimap<pair<double,double>, pair<tDipolePtr,tDipolePtr>,
std::greater< pair<double,double> > > FList;
/**
* Interaction pT of the four partons involved in an interaction. Default order is
* the first and second parton in colour order from the left state first,
* then the first and second parton from right state.
*/
typedef pair< pair<TransverseMomentum, TransverseMomentum>,
pair<TransverseMomentum, TransverseMomentum> > InteractionRecoil;
/**
* A vector of dipole pairs represented as iterators into an FList.
*/
typedef vector<FList::const_iterator> DipolePairVector;
/**
* An ordered map of dipole pairs represented as iterators into an
* FList.
*/
typedef multimap<double,FList::const_iterator,std::greater<double> >
DipolePairMap;
/**
* A interaction of RealPartons, with all needed information about the interaction.
*/
struct RealInteraction {
RealPartonStatePtr lrs, rrs;
DipolePtr d1, d2;
RealPartonPtr p11, p12, p21, p22;
InvEnergy range11, range12, range21, range22;
Energy effectivePlus11, effectivePlus12, effectivePlus21, effectivePlus22;
Energy effectiveMinus11, effectiveMinus12, effectiveMinus21, effectiveMinus22;
bool max11, max12, max21, max22;
double P11, P12, P21, P22;
};
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
inline DipoleXSec()
: theRMax(0.0*InvGeV), theInteraction(0), sinFunction(0), usePartonicInteraction(false),
theIntOrdering(0), theRecoilReduction(0), checkOffShell(true) {}
/**
* The copy constructor.
*/
inline DipoleXSec(const DipoleXSec & x)
: HandlerBase(x), theRMax(x.theRMax), theInteraction(x.theInteraction),
sinFunction(x.sinFunction), usePartonicInteraction(x.usePartonicInteraction),
theIntOrdering(x.theIntOrdering),theRecoilReduction(x.theRecoilReduction),
checkOffShell(x.checkOffShell) {}
/**
* The destructor.
*/
virtual ~DipoleXSec();
//@}
public:
/** @name Virtual functions to be overridden by subclasses. */
//@{
/**
* Calculate the scattering probability for the two given dipoles
* using the given ImpactParameters.
*/
virtual double fij(const pair<tPartonPtr, tPartonPtr>,
const pair<tPartonPtr, tPartonPtr>,
const ImpactParameters & b,
bool veto = true) const;
/**
* Calculate the scattering probability for the two given dipoles
* using the given ImpactParameters.
*/
virtual double fij(const Dipole &,
const Dipole &,
const ImpactParameters & b,
bool veto = true) const;
/**
* Calculate the scattering probability for the two given dipoles
* using the given ImpactParameters.
*/
virtual DipoleInteraction fij(const ImpactParameters & b, Dipole &, Dipole &,
bool veto = true) const;
/**
* Return false if an interaction would be kinematically forbidden.
*/
virtual bool kinematicsVeto(const pair<tPartonPtr, tPartonPtr>,
const pair<tPartonPtr, tPartonPtr>,
const ImpactParameters & b) const;
/**
* Return false if an interaction would be kinematically forbidden.
*/
virtual bool kinematicsVeto(const Dipole &,
const Dipole &,
const ImpactParameters & b,
const pair<bool,bool> & ints) const;
/**
* Return false if an interaction would be kinematically forbidden.
*/
virtual bool kinematicsVeto(const DipoleInteraction &) const;
/**
* Calculate the total scattering probability for the two given
* dipole systems using the given ImpactParameters.
*/
virtual double sumf(const DipoleState &, const DipoleState &,
const ImpactParameters &) const;
/**
* Calculate the total scattering probability for the two given
* dipole systems using the given ImpactParameters.
*/
virtual double sumf(const ImpactParameters &,
const DipoleState &, const DipoleState &) const;
/**
* Calculate the total scattering probability all dipole pairs
* the two given dipole systems using the given ImpactParameters.
*/
virtual FList flist(const DipoleState &, const DipoleState &,
const ImpactParameters &) const;
/**
* Calculate the total scattering probability all dipole pairs
* the two given dipole systems using the given ImpactParameters.
*/
virtual DipoleInteraction::List flist(const ImpactParameters &,
const DipoleState &, const DipoleState &) const;
/**
* Return a unitarized scattering probability, given the
* ununitarized one.
*/
virtual double unitarize(double f) const;
/**
* Returns the recoils the interaction would like to give the 4 involved partons.
* Assumes the states are not yet rotated.
*/
virtual InteractionRecoil
recoil(const pair<tPartonPtr, tPartonPtr> left,
const pair<tPartonPtr, tPartonPtr> right,
const ImpactParameters & b,
pair<pair<bool, bool>, pair<bool, bool> > doesInt,
pair<bool,bool> ints0) const;
/**
* Returns the recoils the interaction would like to give the 4 involved partons.
* Assumes the states are not yet rotated.
*/
virtual InteractionRecoil
recoil(const pair<tPartonPtr, tPartonPtr> left,
const pair<tPartonPtr, tPartonPtr> right,
const ImpactParameters & b,
pair<pair<bool, bool>, pair<bool, bool> > doesInt =
make_pair(make_pair(true, true), make_pair(true, true))) const;
/**
* Returns the recoils the interaction would like to give the 4 involved partons.
* Assumes the states are not yet rotated.
*/
virtual InteractionRecoil
recoil(const DipoleInteraction &) const;
/**
* Does the transverse recoil on the four partons.
*/
void doTransverseRecoils(RealInteraction i, InteractionRecoil recs) const;
/**
* creates and initialises a RealInteraction with realpartons, ranges and max.
*/
virtual RealInteraction
initialiseInteraction(const pair<DipolePtr, DipolePtr> inter,
RealPartonStatePtr lrs, RealPartonStatePtr rrs,
pair<pair<bool, bool>, pair<bool, bool> > doesInt,
const ImpactParameters & b) const;
/**
* creates and initialises a RealInteraction with realpartons, ranges and max.
*/
virtual RealInteraction
initialiseInteraction(const DipoleInteraction & di,
RealPartonStatePtr lrs, RealPartonStatePtr rrs) const;
/**
* Decides which of the four partons actually interact.
*/
virtual pair<pair<bool, bool>, pair<bool, bool> >
doesInt(const pair<tPartonPtr, tPartonPtr> left,
const pair<tPartonPtr, tPartonPtr> right,
const ImpactParameters & b) const;
/**
* Selects the two partons (out of the four) that are used in theInteraction == 0.
*/
virtual pair<bool, bool> int0Partons(tcPartonPtr p11, tcPartonPtr p12,
tcPartonPtr p21, tcPartonPtr p22,
const ImpactParameters & b) const;
/**
* calculates and sets the effective plus and minus of the partons of the interaction
**/
virtual void updateMomenta(RealInteraction * interaction) const;
// /**
// * undoes the last recoil for each of the four partons.
// */
// virtual void undoLastRecoil(RealPartonPtr p11, RealPartonPtr p12,
// RealPartonPtr p21, RealPartonPtr p22) const;
/**
* Does the transverse recoil on the real partons with the pT supplied in recs.
* Does the p+ alt. p- transfer of neededPlus nad neededMinus to put the
* other state on shell. If not consistent, return false.
* Also updates p_mu for the 4 effective partons. Assumes all are rightmoving.
*/
virtual bool doInteraction(InteractionRecoil recs, const FList::const_iterator inter,
RealPartonStatePtr lrs, RealPartonStatePtr rrs,
pair<pair<bool, bool>, pair<bool, bool> > doesInt,
const ImpactParameters & b) const;
/**
* Does the transverse recoil on the real partons with the pT supplied in recs.
* Does the p+ alt. p- transfer of neededPlus nad neededMinus to put the
* other state on shell. If not consistent, return false.
* Also updates p_mu for the 4 effective partons. Assumes all are rightmoving.
*/
virtual bool doInteraction(InteractionRecoil recs,
const DipoleInteraction::List::const_iterator inter,
RealPartonStatePtr lrs, RealPartonStatePtr rrs,
pair<pair<bool, bool>, pair<bool, bool> > doesInt) const;
/**
* reconnects the colour flow in an interaction, taking care that rescatterings are
* treated properly.
*/
bool reconnect(const DipoleInteraction &di) const {
return reconnect(di.dips.first, di.dips.second);
}
/**
* reconnects the colour flow in an interaction, taking care that rescatterings are
* treated properly.
*/
virtual bool reconnect(tDipolePtr d1, tDipolePtr d2) const;
/**
* checks the interactions in the two real states, ans picks out the ones that
* did non-colour singlet exchanges.
*/
virtual vector<pair<DipolePtr, DipolePtr> >
getColourExchanges(tRealPartonStatePtr lrs, tRealPartonStatePtr rrs) const;
/**
* Mark two dipoles as having interacted with eachother.
*/
virtual void interact(Dipole & d1, Dipole & d2) const;
//@}
/**
* Return the value of the sine-interaction strength.
*/
double fSinFn(const Parton::Point & rho1, const Parton::Point & rho2,
const TransverseMomentum & pt) const;
/**
* Return the interaction
**/
inline int interaction() const {
return theInteraction;
}
/**
* Flag determining if only one parton in each dipole is considered
* interacting or both.
*/
inline bool partonicInteraction() const {
return usePartonicInteraction;
}
/**
* The confinement scale.
*/
InvEnergy rMax() const;
/**
* for debugging
*/
mutable int Nfij;
mutable int NScalVeto;
mutable int NBVeto;
mutable double scalVeto;
mutable double bVeto;
private:
/**
* return the fractions of p+ and p- to be transfered to put everything on shell.
*/
pair<double, double> findBoosts(Energy intPlus1, Energy intPlus2,
Energy intMinus1, Energy intMinus2,
Energy evoPlus2, Energy evoMinus1) const;
/**
* return the fractions of p+ and p- to be transfered to put everything on shell.
**/
pair<double, double> findBoosts(RealInteraction i) const;
/**
* boosts the interacting particles with the scales provided
**/
void doBoosts(RealInteraction i, pair<double, double> boosts) const;
/**
* boosts the affective partons with the scale x.
*/
void doBoost(tRealPartonPtr p1, InvEnergy range1,
tRealPartonPtr p2, InvEnergy range2, double x) const;
/**
* reduce the recoil until the interacting partons are ordered in rapidity.
**/
void reduceRecoil(RealInteraction RI, InteractionRecoil recs) const;
/**
* checks if the interaction is ordered or not
**/
bool ordered(RealInteraction i, InteractionRecoil recs, const ImpactParameters & b) const;
/**
* sets all particles in the backwards cone on shell
**/
void setOnShell(RealInteraction i) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* 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.
*/
inline virtual IBPtr clone() const {
return new_ptr(*this);
}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
inline virtual IBPtr fullclone() const {
return new_ptr(*this);
}
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/**
* Exception class to signal bad kinematics.
*/
struct InteractionKinematicException: public Exception {};
private:
/**
* The confinement scale.
*/
InvEnergy theRMax;
/**
* Flag determining which interaction to be used.
*/
int theInteraction;
/**
* Flag determining which approximation to the sine-functions in the
* interaction strength to use.
*/
int sinFunction;
/**
* Flag determining if only one parton in each dipole is considered
* interacting or both.
*/
bool usePartonicInteraction;
/**
* What kind of kinematical ordering is requeired in the interaction.
*/
int theIntOrdering;
/**
* What to do with large recoils, if anything.
*/
int theRecoilReduction;
/**
* Check for off-shell incoming particles.
*/
bool checkOffShell;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DipoleXSec & operator=(const DipoleXSec &);
};
}
#endif /* DIPSY_DipoleXSec_H */
diff --git a/DIPSY/EventFiller.cc b/DIPSY/EventFiller.cc
--- a/DIPSY/EventFiller.cc
+++ b/DIPSY/EventFiller.cc
@@ -1,1657 +1,1662 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the EventFiller class.
//
#include "DipoleEventHandler.h"
#include "EventFiller.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/Utilities/Current.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/EventRecord/Step.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "RealPartonState.h"
#include "ParticleInfo.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/DebugItem.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Handlers/LuminosityFunction.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "../Cascade/EmitterBase.h"
#include <iostream>
#include <fstream>
using namespace DIPSY;
EventFiller::EventFiller()
: currentWeight(0.0), theRecoilScheme(0), theSingleMother(0),
theDGLAPinPT(0), theEffectiveWeights(0), theFSSwingTime(0.0),
theFSSwingTimeStep(0.1),
theValenceChargeNormalisation(0), thePTCut(ZERO), theSoftRemove(1), onlyOnce(false) {}
EventFiller::~EventFiller() {}
IBPtr EventFiller::clone() const {
return new_ptr(*this);
}
IBPtr EventFiller::fullclone() const {
return new_ptr(*this);
}
double EventFiller::fill(Step & step, DipoleEventHandler & eh, tPPair inc,
DipoleState & dl, DipoleState & dr,
const ImpactParameters & b) const {
if ( dl.hasShadows() ) return fillWithShadows(step, eh, inc, b, dl, dr);
if ( mode() == 4 ) return fill(step, eh, inc, b, dl, dr);
// Get a list of possible dipole-dipole interactions
FList fl = eh.xSecFn().flist(dl, dr, b);
//Sum over the (ununitarised) interaction probabilities 2*f_{ij}.
//Note that the factor 2 is included in flist(), to give the correct
//ND interaction probability.
double sum = 0.0;
for ( FList::iterator it = fl.begin(); it != fl.end(); ++it )
sum += it->first.first;
//Just check that they are not all 0.
if ( sum == 0.0 ) {
return 0.0;
}
//Non-diffractive interaction probability is 1 - exp(-Sum(2*f_{ij}))
double weight = eh.xSecFn().unitarize(sum);
//Combine the weights from all the sources.
currentWeight =
weight*sqr(hbarc)*dr.weight()*dl.weight()*b.weight()/eh.maxXSec();
// Select the interactions which should be performed.
pair<RealPartonStatePtr, RealPartonStatePtr> realStates =
selectInteractions(fl, b, eh.xSecFn());
//If no interactions found, discard event.
if ( realStates.first->interactions.empty() ) {
return 0.0;
}
//Counts the number of participants in a HI event.
//TODO: interface, or trigger on AA.
countParticipants(dl, dr, b.bVec().pt());
//Figure out which partons to keep, and which to remove.
//Also sort out colour flow, momenta etc.
vector<String> strings = extractStrings(dl, dr, realStates, b);
DipoleState & finalState = dl;
//Discard event if no final state partons were found.
if ( strings.empty() || ! fillStep(step, inc, strings) ) {
weight = 0.0;
currentWeight = 0.0;
return weight;
}
// The last thing we do is to fix up valens configurations.
finalState.fixValence(step);
return weight;
}
double EventFiller::
fill(Step & step, DipoleEventHandler & eh, tPPair inc,
const ImpactParameters & b, DipoleState & dl, DipoleState & dr) const {
if ( dl.hasShadows() ) return fillWithShadows(step, eh, inc, b, dl, dr);
// Get a list of possible dipole-dipole interactions
InteractionList intl = eh.xSecFn().flist(b, dl, dr);
//Sum over the (ununitarised) interaction probabilities 2*f_{ij}.
//ND interaction probability.
double sum = 0.0;
for ( InteractionList::iterator it = intl.begin(); it != intl.end(); ++it )
sum += it->f2;
//Just check that they are not all 0.
if ( sum <= 0.0 ) {
return 0.0;
}
//Non-diffractive interaction probability is 1 - exp(-Sum(2*f_{ij}))
double weight = eh.xSecFn().unitarize(sum);
//Combine the weights from all the sources.
currentWeight =
weight*sqr(hbarc)*dr.weight()*dl.weight()*b.weight()/eh.maxXSec();
// Select the interactions which should be performed.
pair<RealPartonStatePtr, RealPartonStatePtr> realStates =
selectInteractions(intl, eh.xSecFn());
//If no interactions found, discard event.
if ( realStates.first->interactions.empty() ) {
return 0.0;
}
//Counts the number of participants in a HI event.
//TODO: interface, or trigger on AA.
countParticipants(dl, dr, b.bVec().pt());
//Figure out which partons to keep, and which to remove.
//Also sort out colour flow, momenta etc.
vector<String> strings = extractStrings(dl, dr, realStates, b);
DipoleState & finalState = dl;
//Discard event if no final state partons were found.
if ( strings.empty() || ! fillStep(step, inc, strings) ) {
weight = 0.0;
currentWeight = 0.0;
return weight;
}
// The last thing we do is to fix up valens configurations.
finalState.fixValence(step);
return weight;
}
double EventFiller::
fillWithShadows(Step & step, DipoleEventHandler & eh, tPPair inc,
const ImpactParameters & b, DipoleState & dl, DipoleState & dr) const {
// Get a list of possible dipole-dipole interactions
InteractionList intl = eh.xSecFn().flist(b, dl, dr);
//Sum over the (ununitarised) interaction probabilities 2*f_{ij}.
//ND interaction probability.
double sum = 0.0;
for ( InteractionList::iterator it = intl.begin(); it != intl.end(); ++it )
sum += it->f2;
//Just check that they are not all 0.
if ( intl.empty() || sum <= 0.0 ) return 0.0;
//Non-diffractive interaction probability is 1 - exp(-Sum(2*f_{ij}))
double weight = eh.xSecFn().unitarize(sum);
//Combine the weights from all the sources.
currentWeight =
weight*sqr(hbarc)*dr.weight()*dl.weight()*b.weight()/eh.maxXSec();
InteractionVector interactions = selectShadowInteractions(intl, eh.xSecFn());
//If no interactions found, discard event.
if ( interactions.empty() ) return 0.0;
//Figure out which partons to keep, and which to remove.
//Also sort out colour flow, momenta etc.
vector<String> strings = extractShadowStrings(dl, dr, interactions, b);
DipoleState & finalState = dl;
//Discard event if no final state partons were found.
if ( strings.empty() || ! fillStep(step, inc, strings) ) {
weight = 0.0;
currentWeight = 0.0;
return weight;
}
// The last thing we do is to fix up valens configurations.
finalState.fixValence(step);
return weight;
}
void EventFiller::countParticipants(const DipoleState & dl,
const DipoleState & dr, const InvEnergy b) const {
//The initial dipoles. Note that these will in general have emitted
//something, and thus no longer be active. They may have been
//reconnected by a new dipole in the absorption process though.
vector<DipolePtr> valenceDip = dl.initialDipoles();
valenceDip.insert(valenceDip.end(),
dr.initialDipoles().begin(), dr.initialDipoles().end());
//count the nonparticipating nucleons
int untouchedNucleons = 0;
for ( int i = 0; i < int(valenceDip.size()); i++ ) {
//First find the three valence partons in the nucleon that
//the dipole belongs to.
//Two are connected to the dipole, easy.
tPartonPtr p1 = valenceDip[i]->partons().first;
tPartonPtr p2 = valenceDip[i]->partons().second;
//The three dipoles in valenceDip are always after each other,
//so the third parton can be found in the dipole before or
//after (or both).
tPartonPtr p3;
if ( i != int(valenceDip.size())-1 &&
valenceDip[i+1]->partons().first == p2 )
p3 = valenceDip[i+1]->partons().second;
else if ( i != int(valenceDip.size())-1 &&
valenceDip[i+1]->partons().second == p1 )
p3 = valenceDip[i+1]->partons().first;
else if ( i != 0 && valenceDip[i-1]->partons().first == p2 )
p3 = valenceDip[i-1]->partons().second;
else if ( i != 0 && valenceDip[i-1]->partons().second == p1 )
p3 = valenceDip[i-1]->partons().first;
//If none of the valence partons have interacted
//ie (have children that interacted), the nucleon is considered
//to not have interacted. Note that with swings, a single interaction
//can make several nucleons participating.
if ( !(p1->interacted()) &&
!(p2->interacted()) &&
!( p3 && p3->interacted() ) ) {
valenceDip[i]->participating(false);
untouchedNucleons++;
}
else
valenceDip[i]->participating(true);
}
//We will triple count, as we loop through all 3 dipoles in
//each nucleon. Divide by 3 to make up for this.
untouchedNucleons /= 3;
}
/*
* This is the probability that a parton is interacting, given that
* the state interacts, but the partons is not selected as primary.
*
* Prob of p interacting: p
* Prob of p interacting, given that the state interacts: x = p/totalP
* Prob of p being selected as primary, given state interacts: y = p/sumP
* Prob of p being selected as non-primary, given state
* interacts: z
* x = y + z ==> z = p/totalP - p/sumP
* Prob of p not being primary, given state interacts: A = (sumP-p)/sumP
* And now what we are looking for:
* Prob of p being selected as non-primary, given state interacts
* and not selected as primary: B
* z = A*B ==> B = z/A = p(1/totalP - 1/sumP)/((sumP-p)/sumP) =
* = p(sumP/totalP - 1)/(sumP-p)
*/
double correctedProb(double totalP, double sumP, double p) {
return p*(sumP/totalP - 1.0)/(sumP - p);
}
pair<RealPartonStatePtr, RealPartonStatePtr>
EventFiller::selectInteractions(const FList & fl, const ImpactParameters & b,
const DipoleXSec & xSec) const {
double sumfij = 0.0; //the sum of the individual nonuntarised int probs
double sumUP = 0.0; //The sum of the individual unitarised int probs
//Set up a selector, and a map sorted on pt of the interaction.
Selector<FList::const_iterator, double> sel;
DipolePairMap ordered;
//Loop through all possible interactions and fill
//@sel and @ordered.
for ( FList::const_iterator it = fl.begin(); it != fl.end(); ++it ) {
//If no interaction probability, don't bother.
if ( it->first.second == 0.0 ) continue;
//insert into the selector.
sel.insert(it->first.second, it);
sumfij += it->first.first;
sumUP += it->first.second;
//Calculate interaction recoils by calling the DipoleXSec object.
DipoleXSec::InteractionRecoil rec =
xSec.recoil(it->second.first->partons(), it->second.second->partons(), b);
//Insert in @ordered, sorting on max pt recoil.
double pt = max(max(rec.first.first.pt(), rec.first.second.pt()),
max(rec.second.first.pt(), rec.second.second.pt()))/GeV;
ordered.insert(make_pair(pt, it));
// ordered.insert(make_pair(it->first.first, it));
}
//interaction probability (for at least one interaction)
//of the two states.
double totalP = Current<DipoleEventHandler>()->xSecFn().unitarize(sumfij);
//create the real states.
RealPartonStatePtr rrs = new_ptr(RealPartonState());
RealPartonStatePtr lrs = new_ptr(RealPartonState());
//Add the valence partons (as they are always real) and save.
lrs->addValence(fl.begin()->second.first->dipoleState());
rrs->addValence(fl.begin()->second.second->dipoleState());
lrs->saveState();
rrs->saveState();
DipolePairVector interactions;
DipolePairMap potential;
DipolePairMap failedPrims;
bool found = false;
int counter = 0;
double maxpt = 0.0;
double maxpte = 0.0;
while ( !found ) {
potential.clear();
interactions.clear();
maxpt = 0.0;
maxpte = 0.0;
//select a first interaction (since there has to be at least one).
FList::const_iterator prim = sel[UseRandom::rnd()];
//Go through the other interactions and check if they are also
//interacting. A modified probability is used, to make up for the
//bias introduced in selecting a primary interaction.
for ( FList::const_iterator it = fl.begin(); it != fl.end(); ++it )
if ( it == prim || correctedProb(totalP,sumUP,it->first.second)
> UseRandom::rnd() ) {
DipoleXSec::InteractionRecoil rec =
xSec.recoil(it->second.first->partons(), it->second.second->partons(), b);
double pt = max(max(rec.first.first.pt(), rec.first.second.pt()),
max(rec.second.first.pt(), rec.second.second.pt()))/GeV;
//If the interaction passed the amplitude test, add it to the list
//of potential interations.
potential.insert(make_pair(pt, it));
// potential.insert(make_pair(it->first.first, it));
}
//Keep track on how many times we have tried to find a consistent
//set of interactions. Give up after 10 tries.
counter++;
if ( counter > 10 ) {
overTenEvents += currentWeight;
return make_pair(lrs, rrs);
}
//test all potetential interactions in order, skipping already failed dips.
//failed first interactions inserted in failedprims,
//and from sel and ordered if first in ordered.
int i = 0;
set<tDipolePtr> ldips, rdips;
for ( DipolePairMap::iterator it = potential.begin(); it != potential.end(); ++it ) {
i++;
//Check first interaction, and already failed as first.
//Then autofail without checking again.
if ( failedPrims.find(it->first) != failedPrims.end() && !found ) {
continue;
}
if ( onlyOnce && ( ldips.find(it->second->second.first) != ldips.end() ||
rdips.find(it->second->second.second) != rdips.end() ) ) continue;
//Try to add the interaction.
if (addInteraction(it->second, lrs, rrs, interactions, b, xSec)) {
found = true;
ldips.insert(it->second->second.first);
rdips.insert(it->second->second.second);
maxpt = max(maxpt, it->first);
maxpte = max(maxpte, max(max(it->second->second.first->partons().first->pT().pt(),
it->second->second.first->partons().second->pT().pt()),
max(it->second->second.second->partons().first->pT().pt(),
it->second->second.second->partons().second->pT().pt()))/GeV);
}
else {
//remember if it was first interaction and failed, to not
//try it again later.
if ( !found ) {
failedPrims.insert(*it);
if ( it == ordered.begin() ) {
ordered.erase(it);
sel.erase(it->second);
}
}
}
}
//if sel (or ordered) empty, give up event. no int can be first.
if ( sel.empty() ) return make_pair(lrs, rrs);
}
//Make sure that the real state partons are on-shell.
//May not actually be needed, not sure... TODO: check!
for( DipolePairVector::iterator it = interactions.begin();
it != interactions.end(); it++ ) {
lrs->setOnShell((*it)->second.first);
rrs->setOnShell((*it)->second.second);
}
for ( RealParton::RealPartonSet::iterator it = lrs->valence.begin();
it != lrs->valence.end(); it++ )
(*it)->setOnShell();
for ( RealParton::RealPartonSet::iterator it = rrs->valence.begin();
it != rrs->valence.end(); it++ )
(*it)->setOnShell();
// *** TO REMOVE *** Temporary analysis
if ( histptmax ) {
histptmax->fill(maxpt, currentWeight);
histptmaxe->fill(maxpte, currentWeight);
sumwmaxpt += currentWeight;
}
//return the real states.
return make_pair(lrs, rrs);
}
pair<RealPartonStatePtr, RealPartonStatePtr>
EventFiller::selectInteractions(const InteractionList & intl,
const DipoleXSec & xSec) const {
InteractionVector interactions;
double sumfij = 0.0; //the sum of the individual nonuntarised int probs
double sumUP = 0.0; //The sum of the individual unitarised int probs
//Set up a selector, and a map sorted on pt of the interaction.
Selector<InteractionList::const_iterator, double> sel;
InteractionPTSet ordered;
//Loop through all possible interactions and fill
//@sel and @ordered.
for ( InteractionList::const_iterator it = intl.begin(); it != intl.end(); ++it ) {
//If no interaction probability, don't bother.
if ( it->f2 <= 0.0 ) continue;
//insert into the selector.
sel.insert(it->uf2, it);
sumfij += it->f2;
sumUP += it->uf2;
ordered.insert(it);
}
//interaction probability (for at least one interaction)
//of the two states.
double totalP = Current<DipoleEventHandler>()->xSecFn().unitarize(sumfij);
//create the real states.
RealPartonStatePtr rrs = new_ptr(RealPartonState());
RealPartonStatePtr lrs = new_ptr(RealPartonState());
//Add the valence partons (as they are always real) and save.
lrs->addValence(intl.begin()->dips.first->dipoleState());
rrs->addValence(intl.begin()->dips.second->dipoleState());
lrs->saveState();
rrs->saveState();
InteractionPTSet potential;
InteractionPTSet failedPrims;
bool found = false;
int counter = 0;
while ( !found ) {
potential.clear();
//select a first interaction (since there has to be at least one).
InteractionList::const_iterator prim = sel[UseRandom::rnd()];
//Go through the other interactions and check if they are also
//interacting. A modified probability is used, to make up for the
//bias introduced in selecting a primary interaction.
for ( InteractionList::const_iterator it = intl.begin(); it != intl.end(); ++it )
if ( it == prim || correctedProb(totalP,sumUP,it->uf2) > UseRandom::rnd() )
potential.insert(it);
//Keep track on how many times we have tried to find a consistent
//set of interactions. Give up after 10 tries.
counter++;
if ( counter > 10 ) {
overTenEvents += currentWeight;
return make_pair(lrs, rrs);
}
//test all potetential interactions in order, skipping already failed dips.
//failed first interactions inserted in failedprims,
//and from sel and ordered if first in ordered.
int i = 0;
set<tDipolePtr> ldips, rdips;
set< pair<tPartonPtr, tPartonPtr> > intpairs;
for ( InteractionPTSet::iterator iit = potential.begin();
iit != potential.end(); ++iit ) {
InteractionList::const_iterator it = *iit;
i++;
//Check first interaction, and already failed as first.
//Then autofail without checking again.
if ( failedPrims.find(it) != failedPrims.end() && !found ) {
continue;
}
if ( onlyOnce && ( ldips.find(it->dips.first) != ldips.end() ||
rdips.find(it->dips.second) != rdips.end() ) ) continue;
TransverseMomentum recoil = it->rec;
it->norec = false;
if ( intpairs.find(it->ints) != intpairs.end() ) {
recoil = TransverseMomentum();
it->norec = true;
}
//Try to add the interaction.
if ( addInteraction(it, lrs, rrs, recoil, interactions, xSec) ) {
found = true;
ldips.insert(it->dips.first);
rdips.insert(it->dips.second);
intpairs.insert(it->ints);
}
else {
//remember if it was first interaction and failed, to not
//try it again later.
if ( !found ) {
failedPrims.insert(it);
if ( iit == ordered.begin() ) {
ordered.erase(it);
sel.erase(it);
}
}
}
}
//if sel (or ordered) empty, give up event. no int can be first.
if ( sel.empty() ) return make_pair(lrs, rrs);
}
//Make sure that the real state partons are on-shell.
//May not actually be needed, not sure... TODO: check!
for( InteractionVector::iterator it = interactions.begin();
it != interactions.end(); it++ ) {
lrs->setOnShell((*it)->dips.first);
rrs->setOnShell((*it)->dips.second);
}
for ( RealParton::RealPartonSet::iterator it = lrs->valence.begin();
it != lrs->valence.end(); it++ )
(*it)->setOnShell();
for ( RealParton::RealPartonSet::iterator it = rrs->valence.begin();
it != rrs->valence.end(); it++ )
(*it)->setOnShell();
//return the real states.
return make_pair(lrs, rrs);
}
EventFiller::InteractionVector
EventFiller::selectShadowInteractions(const InteractionList & intl,
const DipoleXSec & xSec) const {
InteractionVector interactions;
tDipoleStatePtr dl = &intl.begin()->dips.first->dipoleState();
tDipoleStatePtr dr = &intl.begin()->dips.second->dipoleState();
double sumfij = 0.0; //the sum of the individual nonuntarised int probs
double sumUP = 0.0; //The sum of the individual unitarised int probs
//Set up a selector, and a map sorted on pt of the interaction.
Selector<InteractionList::const_iterator, double> sel;
InteractionPTSet ordered;
//Loop through all possible interactions and fill
//@sel and @ordered.
for ( InteractionList::const_iterator it = intl.begin(); it != intl.end(); ++it ) {
//If no interaction probability, don't bother.
if ( it->f2 <= 0.0 ) continue;
//insert into the selector.
sel.insert(it->uf2, it);
sumfij += it->f2;
sumUP += it->uf2;
ordered.insert(it);
}
//interaction probability (for at least one interaction)
//of the two states.
double totalP = Current<DipoleEventHandler>()->xSecFn().unitarize(sumfij);
InteractionPTSet potential;
InteractionPTSet failedPrims;
bool found = false;
int counter = 0;
while ( !found ) {
potential.clear();
interactions.clear();
dl->resetShadows();
dr->resetShadows();
//select a first interaction (since there has to be at least one).
InteractionList::const_iterator prim = sel[UseRandom::rnd()];
//Go through the other interactions and check if they are also
//interacting. A modified probability is used, to make up for the
//bias introduced in selecting a primary interaction.
for ( InteractionList::const_iterator it = intl.begin(); it != intl.end(); ++it )
if ( it == prim || correctedProb(totalP,sumUP,it->uf2) > UseRandom::rnd() )
potential.insert(it);
//Keep track on how many times we have tried to find a consistent
//set of interactions. Give up after 10 tries.
counter++;
if ( counter > 10 ) {
overTenEvents += currentWeight;
return interactions;
}
// Test all potetential interactions in order, skipping already
// failed dips. Failed first interactions are inserted in
// failedprims, and removed from sel and ordered if first in
// ordered. Keep track of each dipole that has interacted, as a
// dipole can only interact once. Also keep track of all pairs of
// partons that has itneracted, as if the same paie is
// rescattered, they should not be given double recoil.
int i = 0;
set<tDipolePtr> ldips, rdips;
set< pair<tSPartonPtr, tSPartonPtr> > intpairs;
for ( InteractionPTSet::iterator iit = potential.begin();
iit != potential.end(); ++iit ) {
InteractionList::const_iterator it = *iit;
i++;
// Check first interaction, and if already failed as first
// autofail it without checking again.
if ( failedPrims.find(it) != failedPrims.end() && !found ) continue;
// Never let the same dipole interact more than once.
if ( ldips.find(it->dips.first) != ldips.end() ||
rdips.find(it->dips.second) != rdips.end() ) continue;
// Try to add the interaction. First prepare by inserting the
// interaction at a suitable ShadowParton.
it->prepare();
// If the same two partons interact twice only give recoil once.
it->norec = ( intpairs.find(it->sints) != intpairs.end() );
interactions.push_back(it);
+ it->id = interactions.size();
// We need to check previously accepted interactions as well, so
// we have to recheck all.
- dl->resetInteractedShadows();
- dr->resetInteractedShadows();
- if ( recheckInteractions(interactions) ) {
+ if ( recheckInteractions(interactions, false) ) {
found = true;
ldips.insert(it->dips.first); // Prevent dipoles from interacting again.
rdips.insert(it->dips.second);
intpairs.clear(); // Re-build all pairs of interacting partons.
for ( int i = 0, N = interactions.size(); i < N; ++i )
intpairs.insert(interactions[i]->sints);
} else {
interactions.pop_back();
it->reject();
+ it->id = -it->id;
+ // *** TODO *** remove this - it's only for debugging.
+ if ( !recheckInteractions(interactions, false) )
+ Throw<ConsistencyException>()
+ << "In DIPSY::EventFiller: A previously accepted set of interactions was not "
+ << "accepted in subsequent debug check!" << Exception::abortnow;
//remember if it was first interaction and failed, to not
//try it again later.
if ( !found ) {
failedPrims.insert(it);
if ( iit == ordered.begin() ) {
ordered.erase(it);
sel.erase(it);
}
}
}
}
//if sel (or ordered) empty, give up event. no int can be first.
if ( sel.empty() ) return interactions;
}
return interactions;
}
bool EventFiller::
addInteraction(FList::const_iterator inter, RealPartonStatePtr lrs,
RealPartonStatePtr rrs, DipolePairVector & inters,
const ImpactParameters & b, const DipoleXSec & xSec) const {
rescatter += currentWeight;
//With some settings, only some of the partons actually interact.
//Check which here. Other tunes will just return 4x true.
pair<pair<bool, bool>, pair<bool, bool> > doesInt =
xSec.doesInt(inter->second.first->partons(),
inter->second.second->partons(), b);
//Calculate the recoils from the interaction.
DipoleXSec::InteractionRecoil recs =
xSec.recoil(inter->second.first->partons(),
inter->second.second->partons(), b, doesInt);
//Add the interacting partons and their parents to the real state
//and go through the evolution checking for ordering, local pt-max, etc.
if ( !lrs->fullControlEvolution
(inter->second.first, inter->second.second,
doesInt.first.first, doesInt.first.second,
recs.first.first.pt(), recs.first.second.pt()) ) {
return false;
}
//Add the interaction to the other state, and check that as well.
if ( !rrs->fullControlEvolution
(inter->second.second, inter->second.first,
doesInt.second.first, doesInt.second.second,
recs.second.first.pt(), recs.second.second.pt()) ) {
//If this second evolution failed, remove the interaction from
//the first
//state by reverting it to the previously saved state.
lrs->revertToPrevious(inter->second.first);
return false;
}
//If both evolutions could accomodate the interacting partons,
//add the interaction recoil as well,
//and check that things still are ok.
inters.push_back(inter);
if ( !controlRecoils(inters, lrs, rrs, b, xSec, doesInt) ) {
//If failed, remove the interaction from the evolutions.
lrs->revertToPrevious(inter->second.first);
rrs->revertToPrevious(inter->second.second);
inters.pop_back();
return false;
}
//If the interactions was ok, mark the dipoles as interacting
//and save the real states.
xSec.interact(*inter->second.first, *inter->second.second);
// inter->second.first->interact(*inter->second.second);
// inter->second.second->interact(*inter->second.first);
lrs->saveState();
rrs->saveState();
return true;
}
bool EventFiller::
addInteraction(InteractionList::const_iterator inter, RealPartonStatePtr lrs,
RealPartonStatePtr rrs, const TransverseMomentum & recoil,
InteractionVector & inters, const DipoleXSec & xSec) const {
rescatter += currentWeight;
//Add the interacting partons and their parents to the real state
//and go through the evolution checking for ordering, local pt-max, etc.
if ( !lrs->singleControlEvolution(inter->dips.first, inter->dips.second,
inter->ints.first, recoil) )
return false;
//Add the interaction to the other state, and check that as well.
if ( !rrs->singleControlEvolution(inter->dips.second, inter->dips.first,
inter->ints.second, inter->b->invRotatePT(-recoil)) )
return false;
//If both evolutions could accomodate the interacting partons,
//add the interaction recoil as well,
//and check that things still are ok.
inters.push_back(inter);
if ( !controlRecoils(inters, lrs, rrs, xSec) ) {
//If failed, remove the interaction from the evolutions.
lrs->revertToPrevious(inter->dips.first);
rrs->revertToPrevious(inter->dips.second);
inters.pop_back();
return false;
}
//If the interactions was ok, mark the dipoles as interacting
//and save the real states.
xSec.interact(*inter->dips.first, *inter->dips.second);
// inter->second.first->interact(*inter->second.second);
// inter->second.second->interact(*inter->second.first);
lrs->saveState();
rrs->saveState();
return true;
}
-bool EventFiller::recheckInteractions(const InteractionVector & interactions) const {
- for ( int i = 0, N = interactions.size(); i < N; ++i )
- if ( !interactions[i]->check() ) return false;
- else interactions[i]->accept();
+bool EventFiller::recheckInteractions(const InteractionVector & interactions, bool final) const {
+ for ( int i = 0, N = interactions.size(); i < N; ++i ) {
+ if ( i == 0 ) {
+ interactions[0]->dips.first->dipoleState().resetInteractedShadows();
+ interactions[0]->dips.second->dipoleState().resetInteractedShadows();
+ }
+ if ( !interactions[i]->check(final) ) return false;
+ }
return true;
}
vector<EventFiller::String>
EventFiller::extractStrings(DipoleState & dl, DipoleState & dr,
pair<RealPartonStatePtr, RealPartonStatePtr> realStates,
const ImpactParameters & b ) const {
//Just rename a bit for convenience.
DipoleStatePtr rightState = &dr;
DipoleStatePtr leftState = &dl;
RealPartonStatePtr lrs = realStates.first;
RealPartonStatePtr rrs = realStates.second;
//grow back the partons marked as virtuals.
//This will be the ones without interacting children, but not the
//ones removed during ordering/pt-max-fixing.
removeVirtuals(leftState);
removeVirtuals(rightState);
//The unordered and not-ok pt-max were removed by pairing them
//up with other partons that should stay. Now merge all the momentum
//into a single parton, and flag the others as virtual.
lrs->mergeVirtuals();
rrs->mergeVirtuals();
//Save to update the parton information from the realPartons.
lrs->saveState();
rrs->saveState();
//A high-pt parton (so very localised in x_T) can not recoil all
//of a low-pt (so smeared out in x_T) parent, but will rather shoot
//out a recoiling part of the parent. Fix this by replacing recoiled
//low-pt parents by a remnant, and a recoiler.
lrs->addRecoilers();
rrs->addRecoilers();
//Save states to update partons.
lrs->saveState();
rrs->saveState();
//balance p+/p- by moving the entire state. Was first done in the
//interaction (where the interacting partons had to provide the
//energy), but the modifications here changes kinematics so that
//some more balancing may be needed. Should very rarely be
//large changes.
fixBoost(lrs, rrs);
//Translate the right state to it's right place in x_T,
//and mirror it in rapidity so that it actually comes from the
//other side. Finally merge the two cascades into one state.
rightState->translate(b);
rightState->mirror(0.0);
DipoleStatePtr finalState = leftState->merge(rightState);
vector<pair<DipolePtr, DipolePtr> > swinging =
Current<DipoleEventHandler>()->xSecFn().getColourExchanges(lrs, rrs);
//swing the interacting dipoles.
for ( int i = 0; i < int(swinging.size()); i++ ) {
Current<DipoleEventHandler>()->
xSecFn().reconnect(swinging[i].first, swinging[i].second);
}
//now remove the partons that were made virtual by mergeVirtuals.
//This is so that it should be easier to reconnect the colour
//through the interacting dipoles with the state from the other side.
removeVirtuals(finalState);
//do final state swings, aka pythias colour reconnection.
double yrange = FSSwingTime();
double ystep = FSSwingTimeStep();
for ( double y = 0.0; y < yrange; y += ystep ) {
finalState->swingFS(y, y + ystep);
}
// finalState->lambdaMeasure(0.36*GeV2, histDipLength2, histDipMass2);
//compensate for the 6 valence charges in the proton
finalState->normaliseValenceCharge(theValenceChargeNormalisation);
//make sure the non-particpating nucleons (if AA) have the
//original colour flow.
finalState->restoreNonparticipants();
//try to find and fix any problems or inconsistencies that may have
//popped up. Error messages are sent if anything is found.
dodgeErrors(finalState);
//If you want to save the initial state event to file (mainly AA).
//TODO: interface!
//TODO: stop event after this.
//finalState->saveGluonsToFile(currentWeight);
//If you want to make a fancy movie in mathematica of your event. :)
// finalState->printForMovie(0, 1000);
//Sort the remaining partons in strings.
vector<String> ret = finalState->strings();
static DebugItem printfinal("DIPSY::PrintFinal", 6);
if ( histptfi ) {
for ( int i = 0, N = ret.size(); i < N; ++i )
for ( int j = 0, M = ret[i].size(); j < M; ++j ) {
if ( abs(ret[i][j]->y()) < 1.0 ) {
histptf->fill(ret[i][j]->pT().pt()/GeV, currentWeight);
if ( ret[i][j]->interacted() )
histptfi->fill(ret[i][j]->pT().pt()/GeV, currentWeight);
}
if ( printfinal ) {
cerr << ( ret[i][j]->interacted()? "*": "-")
<< ( ret[i][j]->valence()? "v": " ")
<< ( ret[i][j]->rightMoving()? ">": "<")
<< setw(9) << ret[i][j]->y()
<< setw(10) << ret[i][j]->pT().pt()/GeV
<< setw(10) << ret[i][j]->position().x()/InvGeV
<< setw(10) << ret[i][j]->position().y()/InvGeV
<< endl;
}
}
if ( printfinal ) cerr << endl;
}
return ret;
}
vector<EventFiller::String>
EventFiller::extractShadowStrings(DipoleState & dl, DipoleState & dr,
const InteractionVector & interactions,
const ImpactParameters & b ) const {
// First we check all interactions again, this time putting relevant
// partons on shell.
- dl.resetInteractedShadows();
- dr.resetInteractedShadows();
-
- for ( int i = 0, N = interactions.size(); i < N; i++ )
- if ( !interactions[i]->perform() )
- Throw<ConsistencyException>()
- << "In DIPSY::EventFiller: A previously accepted set of interactions was not "
- << "accepted in the final check!" << Exception::abortnow;
+ if ( !recheckInteractions(interactions, true) )
+ Throw<ConsistencyException>()
+ << "In DIPSY::EventFiller: A previously accepted set of interactions was not "
+ << "accepted in the final check!" << Exception::abortnow;
// Translate the right state to it's right place in x_T,
// and mirror it in rapidity so that it actually comes from the
// other side. Finally merge the two cascades into one state.
dr.translate(b);
dr.mirror(0.0);
DipoleStatePtr finalState = dl.merge(&dr);
//swing the interacting dipoles.
for ( int i = 0, N = interactions.size(); i < N; i++ )
Current<DipoleEventHandler>()->xSecFn().reconnect(*interactions[i]);
// Remove all virtual partons
removeVirtuals(finalState);
// Compensate for the 6 valence charges in the proton
finalState->normaliseValenceCharge(theValenceChargeNormalisation);
// Make sure the non-particpating nucleons (if AA) have the
// original colour flow.
finalState->restoreNonparticipants();
//try to find and fix any problems or inconsistencies that may have
//popped up. Error messages are sent if anything is found.
dodgeErrors(finalState);
//Sort the remaining partons in strings.
vector<String> ret = finalState->strings();
static DebugItem printfinal("DIPSY::PrintFinal", 6);
if ( histptfi ) {
for ( int i = 0, N = ret.size(); i < N; ++i )
for ( int j = 0, M = ret[i].size(); j < M; ++j ) {
if ( abs(ret[i][j]->y()) < 1.0 ) {
histptf->fill(ret[i][j]->pT().pt()/GeV, currentWeight);
if ( ret[i][j]->interacted() )
histptfi->fill(ret[i][j]->pT().pt()/GeV, currentWeight);
}
if ( printfinal ) {
cerr << ( ret[i][j]->interacted()? "*": "-")
<< ( ret[i][j]->valence()? "v": " ")
<< ( ret[i][j]->rightMoving()? ">": "<")
<< setw(9) << ret[i][j]->y()
<< setw(10) << ret[i][j]->pT().pt()/GeV
<< setw(10) << ret[i][j]->position().x()/InvGeV
<< setw(10) << ret[i][j]->position().y()/InvGeV
<< endl;
}
}
if ( printfinal ) cerr << endl;
}
return ret;
}
bool EventFiller::fillStep (Step & step, tPPair incoming,
const vector<String> & strings) const {
Direction<0> dir(true);
vector< vector<PPtr> > particles(strings.size());
SubProPtr sub = new_ptr(SubProcess(incoming));
for ( int is = 0, NS = strings.size(); is < NS; ++is ) {
int N = strings[is].size();
vector<bool> removed(N, false);
if ( theSoftRemove && pT2Cut() > ZERO && N > 2 ) {
bool changed = true;
while ( changed ) {
changed = false;
for ( int i = 0; i < N; ++i ) {
if ( removed[i] ) continue;
if ( strings[is][i]->valence() && theSoftRemove == 2 ) continue;
if ( strings[is][i]->flavour() != ParticleID::g ) continue;
int i1 = (N + i - 1)%N;
while ( removed[i1] ) i1 = (N + i1 - 1)%N;
int i3 = (i + 1)%N;
while ( removed[i3] ) i3 = (i3 + 1)%N;
if ( i1 == i3 ) continue;
if ( invPT2(strings[is], i1, i, i3) < pT2Cut() ) {
if ( removeGluon(strings[is], i1, i, i3) ) {
removed[i] = changed = true;
} else {
return false;
}
}
}
}
}
particles[is].resize(N);
for ( int i = 0; i < N; ++i ) {
if ( removed[i] ) continue;
particles[is][i] = strings[is][i]->produceParticle();
int ip = i - 1;
while ( ip >= 0 && removed[ip] ) --ip;
if ( ip >= 0 )
particles[is][i]->colourConnect(particles[is][ip]);
if ( i == N - 1 && strings[is][i]->flavour() == ParticleID::g ) {
int i0 = 0;
while ( i0 < i && removed[i0] ) ++i0;
particles[is][i0]->colourConnect(particles[is][i]);
}
sub->addOutgoing(particles[is][i]);
}
if ( strings[is][0]->flavour() == ParticleID::g ) {
int i0 = 0;
int in = N - 1;
while ( i0 < in && removed[i0] ) ++i0;
while ( in > i0 && removed[in] ) --in;
particles[is][i0]->colourConnect(particles[is][in]);
}
}
step.addSubProcess(sub);
return true;
}
void EventFiller::fixBoost(RealPartonStatePtr lrs, RealPartonStatePtr rrs) const {
//access the CoM energy, through a kindof roundabout path...
//TODO: must be a better way! Current?
PartonPtr dummy = lrs->partons.begin()->first;
Energy sqrtS = ((dummy->dipoles().first) ? dummy->dipoles().first:dummy->dipoles().second)->dipoleState().handler().lumiFn().maximumCMEnergy();
//sum left and right p+ and p-.
Energy leftPlus = ZERO;
Energy leftMinus = ZERO;
Energy rightPlus = ZERO;
Energy rightMinus = ZERO;
//loop and sum over the to-keep partons.
for ( map<tPartonPtr,RealPartonPtr>::const_iterator
it =lrs->partons.begin();
it != lrs->partons.end(); it++ ) {
if ( it->second->keep == RealParton::NO ) continue;
leftPlus += it->second->plus;
leftMinus += it->second->minus;
}
for ( map<tPartonPtr,RealPartonPtr>::const_iterator
it = rrs->partons.begin();
it != rrs->partons.end(); it++ ) {
if ( it->second->keep == RealParton::NO ) continue;
rightPlus += it->second->minus;
rightMinus += it->second->plus;
}
//Solve the 2:nd degree equation for how much energy has to be
//transfered to set both states on shell.
double A = (- rightPlus*rightMinus - sqr(sqrtS) + leftPlus*leftMinus)/(2.0*rightMinus*sqrtS);
double B = rightPlus/rightMinus;
double y1 = -1.0;
double x1 = -1.0;
if ( sqr(A) - B > 0.0 ) {
//the factor to change right p- with.
y1 = - A + sqrt(sqr(A) - B);
// double y2 = - A - sqrt(sqr(A) - B);
//The factor to change left p+ with.
x1 = (y1*sqrtS - rightPlus)/(y1*leftPlus);
// double x2 = (y2*sqrtS - rightPlus)/(y2*leftPlus);
}
//error handling
if ( x1 < 0 || y1 < 0 ) {
Throw<SpaceLikeGluons>()
<< "EventFiller::fixBoost gave negative or nan solution to boost equation, "
<< "will not balance momentum.\n"
<< "This is probably caused by a rouge gluon being far too virtual.\n"
<< Exception::warning;
return;
}
//loop again, scaling the p+ and p- according to the solution above.
for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = lrs->partons.begin();
it != lrs->partons.end(); it++) {
if ( it->second->keep == RealParton::NO ) continue;
it->second->plus *= x1;
it->second->minus /= x1;
}
for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = rrs->partons.begin();
it != rrs->partons.end(); it++) {
if ( it->second->keep == RealParton::NO ) continue;
it->second->minus /= y1;
it->second->plus *= y1;
}
//Save state to transfer the changes to the partons.
lrs->saveState();
rrs->saveState();
}
// void EventFiller::fixValence(Step & step, DipoleState & dl, DipoleState & dr) const {
// list<PartonPtr> alll = dl.getPartons();
// set<tcPartonPtr> vall;
// for ( list<PartonPtr>::iterator pit = alll.begin(); pit != alll.end(); ++pit )
// if ( (**pit).valence() ) vall.insert(*pit);
// list<PartonPtr> allr = dr.getPartons();
// set<tcPartonPtr> valr;
// for ( list<PartonPtr>::iterator pit = allr.begin(); pit != allr.end(); ++pit )
// if ( (**pit).valence() ) valr.insert(*pit);
// vector<PPtr> lval;
// vector<PPtr> rval;
// for ( ParticleSet::iterator pit = step.particles().begin();
// pit != step.particles().end(); ++pit ) {
// tcPartonPtr p = ParticleInfo::getParton(**pit);
// if ( !p ) continue;
// if ( member(vall, p) ) lval.push_back(*pit);
// else if ( member(valr, p) ) rval.push_back(*pit);
// }
// if ( UseRandom::rndbool() ) {
// dl.fixValence(step, lval);
// dr.fixValence(step, rval);
// } else {
// dr.fixValence(step, rval);
// dl.fixValence(step, lval);
// }
// }
void EventFiller::dodgeErrors(DipoleStatePtr finalState) const {
list<PartonPtr> partons = finalState->getPartons();
for ( list<PartonPtr>::iterator it = partons.begin(); it != partons.end(); it++ ) {
PartonPtr p = *it;
//check for empty pointers
if ( !p ) {
Throw<SpaceLikeGluons>()
<< "dodgeError found empty pointer from getPartons()! :o" << Exception::warning;
continue;
}
//check for colour flow to itself
//check for 0 monetum
if ( p->pT().pt() == ZERO ) {
Throw<SpaceLikeGluons>()
<< "dodgeError found 0 pt gluon." << Exception::warning;
p->pT(TransverseMomentum(UseRandom::rnd()*GeV,UseRandom::rnd()*GeV));
}
if ( p->plus() == ZERO || p->minus() == ZERO ) {
Throw<SpaceLikeGluons>()
<< "dodgeError found 0 lightcone momentum." << Exception::warning;
p->minus(UseRandom::rnd()*GeV);
p->plus(UseRandom::rnd()*GeV);
}
//check for nan momentum
if ( isnan(p->pT().pt()/GeV) ) {
Throw<SpaceLikeGluons>()
<< "dodgeError found NAN pt gluon, fix pt." << Exception::warning;
p->pT(TransverseMomentum(UseRandom::rnd()*GeV,UseRandom::rnd()*GeV));
}
if ( isnan(p->plus()/GeV) || isnan(p->minus()/GeV) ) {
Throw<SpaceLikeGluons>()
<< "dodgeError found NAN lightcone momentum, fix plus, minus." << Exception::warning;
p->minus(UseRandom::rnd()*GeV);
p->plus(UseRandom::rnd()*GeV);
}
//check for negative momentum
if ( p->pT().pt() < ZERO ) {
Throw<SpaceLikeGluons>()
<< "dodgeError found negative pt gluon.... >_>, fix pt." << Exception::warning;
p->pT(TransverseMomentum(UseRandom::rnd()*GeV,UseRandom::rnd()*GeV));
}
if ( p->plus() < ZERO || p->minus() < ZERO ) {
Throw<SpaceLikeGluons>()
<< "dodgeError found negative lightcone momentum, fix plus,minus." << Exception::warning;
p->minus(UseRandom::rnd()*GeV);
p->plus(UseRandom::rnd()*GeV);
//check for off-shell momentum
if ( sqr(p->plus()*p->minus() - p->pT().pt2() - sqr(p->mass())) > sqr(sqr(0.01*GeV)) ) {
Throw<SpaceLikeGluons>()
<< "dodgeError found off-shell parton, fix pt." << Exception::warning;
if ( p->plus()*p->minus() < sqr(p->mass()) ) {
Throw<SpaceLikeGluons>()
<< "dodgeError found insufficient energy for mass, fix plus,minus." << Exception::warning;
double ratio = 2.0*sqr(p->mass())/(p->plus()*p->minus()); //give a bit extra for pt
p->plus(p->plus()*sqrt(ratio));
p->minus(p->minus()*sqrt(ratio));
}
double mod = sqrt(p->plus()*p->minus() - sqr(p->mass()))/p->pT().pt();
p->pT(p->pT()*mod);
}
}
}
}
bool EventFiller::
controlRecoils(DipolePairVector & sel,
RealPartonStatePtr lrs, RealPartonStatePtr rrs,
const ImpactParameters & b, const DipoleXSec & xSec,
pair<pair<bool, bool>, pair<bool, bool> > doesInt) const {
//Keep track on which partons are interacting.
list<pair<bool, bool> >::const_iterator
leftDoesInt = lrs->doesInts.begin();
list<pair<bool, bool> >::const_iterator
rightDoesInt = rrs->doesInts.begin();
//Loop through the interaction, and
for ( DipolePairVector::const_iterator inter = sel.begin();
inter != sel.end(); inter++ ) {
//Which partons are interacting in this dip-dip interaction
pair<pair<bool, bool>, pair<bool, bool> >
trueDoesInt = make_pair(*leftDoesInt, *rightDoesInt);
//calculate recoils
DipoleXSec::InteractionRecoil recoil =
xSec.recoil((*inter)->second.first->partons(),
(*inter)->second.second->partons(), b, trueDoesInt);
//call the DipoleXSec object to check kinematics and vetos
//for the interaction. If things don't work, undo the changes.
if ( !xSec.doInteraction(recoil, *inter, lrs, rrs, trueDoesInt, b) ) {
//revert the total recoil between the states.
lrs->totalRecoil -= recoil.first.first + recoil.first.second;
rrs->totalRecoil -= recoil.second.first + recoil.second.second;
}
}
return true;
}
bool EventFiller::
controlRecoils(InteractionVector & sel, RealPartonStatePtr lrs, RealPartonStatePtr rrs,
const DipoleXSec & xSec) const {
//Keep track on which partons are interacting.
list<pair<bool, bool> >::const_iterator
leftDoesInt = lrs->doesInts.begin();
list<pair<bool, bool> >::const_iterator
rightDoesInt = rrs->doesInts.begin();
//Loop through the interaction, and
for ( InteractionVector::const_iterator inter = sel.begin();
inter != sel.end(); inter++ ) {
//Which partons are interacting in this dip-dip interaction
pair<pair<bool, bool>, pair<bool, bool> >
trueDoesInt = make_pair(*leftDoesInt, *rightDoesInt);
//calculate recoils
DipoleXSec::InteractionRecoil recoil = xSec.recoil(**inter);
//call the DipoleXSec object to check kinematics and vetos
//for the interaction. If things don't work, undo the changes.
if ( !xSec.doInteraction(recoil, *inter, lrs, rrs, trueDoesInt) ) {
//revert the total recoil between the states.
lrs->totalRecoil -= recoil.first.first + recoil.first.second;
rrs->totalRecoil -= recoil.second.first + recoil.second.second;
}
}
return true;
}
void EventFiller::removeVirtuals(DipoleStatePtr state) const {
list<PartonPtr> vP = state->getPartons();
//loop through the partons, and remove the off shell ones.
//this only does colour flow. RealPartonState does the kinematics
//in mergeVirtuals.
for ( list<PartonPtr>::iterator it = vP.begin(); it != vP.end(); it++) {
tPartonPtr p = *it;
//only handle off-shell parton.
if ( !(p->onShell()) ) {
//if there are only 2 or fewer on-shell partons in the colour chain,
//it has to be swinged at some point, and better early than late,
//as last-second forced swings can lead to silly colour connections.
if ( p->nOnShellInChain() < 2 ) {
//tell the absorber to find a swing anywhere in the colour chain
if ( p->dipoles().first ) absorber()->swingLoop(p->dipoles().first, *state);
else absorber()->swingLoop(p->dipoles().second, *state);
}
//once off-shell loops are handled, absorb the parton.
absorber()->removeParton(p);
}
}
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void EventFiller::doinit() throw(InitException) {
HandlerBase::doinit();
}
void EventFiller::dofinish() {
HandlerBase::dofinish();
if ( sumwmaxpt > 0.0 ) {
histptmax->scale(1.0/sumwmaxpt);
histptmaxe->scale(1.0/sumwmaxpt);
histptf->scale(1.0/sumwmaxpt);
histptfi->scale(1.0/sumwmaxpt);
}
}
void EventFiller::doinitrun() {
HandlerBase::doinitrun();
histptmax = FactoryBase::tH1DPtr();
histptmaxe = FactoryBase::tH1DPtr();
histptf = FactoryBase::tH1DPtr();
histptfi = FactoryBase::tH1DPtr();
sumwmaxpt = 0.0;
static DebugItem plotmaxpt("DIPSY::PlotMaxPT", 6);
if ( plotmaxpt && generator()->histogramFactory() ) {
generator()->histogramFactory()->initrun();
generator()->histogramFactory()->registerClient(this);
histptmax = generator()->histogramFactory()->createHistogram1D
("ptmi",100,0.0,100.0);
histptmaxe = generator()->histogramFactory()->createHistogram1D
("ptme",100,0.0,100.0);
histptf = generator()->histogramFactory()->createHistogram1D
("ptfa",100,0.0,100.0);
histptfi = generator()->histogramFactory()->createHistogram1D
("ptfi",100,0.0,100.0);
}
}
double EventFiller::pTScale(DipoleState & state) const {
return state.handler().emitter().pTScale();
}
Energy2 EventFiller::invPT2(const String & str, int i1, int i2, int i3) {
LorentzMomentum p1 = str[i1]->momentum();
LorentzMomentum p2 = str[i2]->momentum();
LorentzMomentum p3 = str[i3]->momentum();
Energy2 s = (p1 + p2 + p3).m2();
if ( s < sqr(str[i1]->mass() + str[i3]->mass() ) )
Throw<SpaceLikeGluons>()
<< "DIPSY produced space-like gluons. Three neighboring ones had "
<< "negative mass squared. This cannot be fixed at the moment. "
<< "Event discarded." << Exception::eventerror;
Energy2 s12 = (p1 + p2).m2();
Energy2 s23 = (p2 + p3).m2();
if ( s12 < ZERO || s23 < ZERO ) return -1.0*GeV2;
return s12*s23/s;
}
bool EventFiller::removeGluon(const String & str, int i1, int i2, int i3) {
LorentzMomentum p1 = str[i1]->momentum();
LorentzMomentum p2 = str[i2]->momentum();
LorentzMomentum p3 = str[i3]->momentum();
LorentzMomentum ptest = p1 + p2;
if ( ptest.m2() < Constants::epsilon*1000.0*sqr(ptest.e()) ) {
str[i1]->plus(ptest.plus());
str[i1]->pT(TransverseMomentum(ptest.x(), ptest.y()));
return true;
}
ptest = p3 + p2;
if ( ptest.m2() < Constants::epsilon*1000.0*sqr(ptest.e()) ) {
str[i3]->plus(ptest.plus());
str[i3]->pT(TransverseMomentum(ptest.x(), ptest.y()));
return true;
}
Energy2 S = (p1 + p2 + p3).m2();
if ( S <= ZERO ) return false;
LorentzRotation R = Utilities::boostToCM(makeTriplet(&p1, &p2, &p3));
double Psi = Constants::pi - p3.theta();
double beta = 0.0;
Energy W = sqrt(S);
double x1 = 2.0*p1.e()/W;
double x3 = 2.0*p3.e()/W;
bool g1 = str[i1]->flavour() == ParticleID::g;
bool g3 = str[i3]->flavour() == ParticleID::g;
if ( ( g1 && g3 ) || (!g1 && !g3 ) )
beta = Psi*sqr(x3)/(sqr(x1) + sqr(x3)); // minimize pt
else if ( g1 )
beta = Psi;
R.rotateY(-beta);
R.invert();
Lorentz5Momentum p1n(p1.m());
Lorentz5Momentum p3n(p3.m());
try {
SimplePhaseSpace::CMS(p1n, p3n, S, 1.0, 0.0);
} catch ( ImpossibleKinematics ) {
return false;
}
p1n.transform(R);
p3n.transform(R);
str[i1]->plus(p1n.plus());
str[i1]->pT(TransverseMomentum(p1n.x(), p1n.y()));
str[i3]->plus(p3n.plus());
str[i3]->pT(TransverseMomentum(p3n.x(), p3n.y()));
return true;
}
void EventFiller::persistentOutput(PersistentOStream & os) const {
os << theAbsorber << theRecoilScheme << theMode << theSingleMother
<< theDGLAPinPT << theEffectiveWeights << theFSSwingTime
<< theFSSwingTimeStep << theValenceChargeNormalisation
<< ounit(thePTCut, GeV) << theSoftRemove << onlyOnce;
}
void EventFiller::persistentInput(PersistentIStream & is, int) {
is >> theAbsorber >> theRecoilScheme >> theMode >> theSingleMother
>> theDGLAPinPT >> theEffectiveWeights >> theFSSwingTime
>> theFSSwingTimeStep >> theValenceChargeNormalisation
>> iunit(thePTCut, GeV) >> theSoftRemove >> onlyOnce;
}
DescribeClass<EventFiller,HandlerBase>
describeDIPSYEventFiller("DIPSY::EventFiller", "libAriadne5.so libDIPSY.so");
void EventFiller::Init() {
static ClassDocumentation<EventFiller> documentation
("The EventFiller class is able to produce an initial ThePEG::Step "
"from two colliding DipoleStates.");
static Reference<EventFiller,DipoleAbsorber> interfaceDipoleAbsorber
("DipoleAbsorber",
"The object used to absorb non-interacting dipoles.",
&EventFiller::theAbsorber, true, false, true, true, false);
static Switch<EventFiller,int> interfaceMode
("Mode",
"How the real state is found from the virtual cascade. Speed versus consistency.",
&EventFiller::theMode, 0, true, false);
static SwitchOption interfaceModeMode
(interfaceMode,
"Consistent",
"The fully consistent version. Not currently reccomended.",
0);
static SwitchOption interfaceModeFast
(interfaceMode,
"Fast",
"Checking only the new real partons introduced by each new interaction, rather than rechecking them all. Good for heavy ions. Not currently reccomended.",
1);
static SwitchOption interfaceModeSingle
(interfaceMode,
"SingleSweep",
"Checking all partons, but sweeping just once in each direction, making "
"it a lot faster. Good for heavy ions. May give an occasional unordered "
"chain, but hopefully not too often. The recommended option.",
2);
static SwitchOption interfaceModeNonRecursive
(interfaceMode,
"NonRecursive",
"Does all evo no matter what. Then removes biggest problem, and does it "
"all over again. Never turns partons back on once switched off.",
3);
static SwitchOption interfaceModeNewSingle
(interfaceMode,
"NewSingle",
"A new implementation of SingleSweep.",
4);
static Switch<EventFiller,int> interfaceEffectiveWeights
("EffectiveWeights",
"How the p+ and pT recoils are distributed among the single partons in an effective parton.",
&EventFiller::theEffectiveWeights, 0, true, false);
static SwitchOption interfaceEffectiveWeightsPlusWeighted
(interfaceEffectiveWeights,
"PlusWeighted",
"Weight pt and plus according to p+ of the individual partons.",
0);
static SwitchOption interfaceEffectiveWeightsPlusEvenWeighted
(interfaceEffectiveWeights,
"PlusEvenWeighted",
"The plus is distibuted according to the plus of the partons, but the pt is shared evenly among the partons.",
1);
static SwitchOption interfaceEffectiveWeightsPlusSingleWeighted
(interfaceEffectiveWeights,
"PlusSingleWeighted",
"The plus is distibuted according to the plus of the partons, but the pt is taken only by the colour connected parton",
2);
static Switch<EventFiller,int> interfaceRecoilScheme
("RecoilScheme",
"How to give tread the positive light-cone momentum of a recoiling parton.",
&EventFiller::theRecoilScheme, 0, true, false);
static SwitchOption interfaceRecoilSchemePreservePlus
(interfaceRecoilScheme,
"PreservePlus",
"blaha",
0);
static SwitchOption interfaceRecoilSchemeFixedY
(interfaceRecoilScheme,
"FixedY",
"blahaaa",
1);
static SwitchOption interfaceRecoilSchemeFrameFan
(interfaceRecoilScheme,
"FrameFan",
"dssklajhls",
2);
static Parameter<EventFiller,int> interfaceSingleMother
("SingleMother",
"If an emission is regarded to come from a single parton rather than both "
" partons in the emitting dipole.",
&EventFiller::theSingleMother, 1, 1, 0, 0,
true, false, Interface::lowerlim);
static Parameter<EventFiller,int> interfaceDGLAPinPT
("DGLAPinPT",
"If the DGLAP supression should be made in terms of pt rather than r. "
" partons in the emitting dipole.",
&EventFiller::theDGLAPinPT, 1, 1, 0, 0,
true, false, Interface::lowerlim);
static Parameter<EventFiller,double> interfaceFSSwingTimeStep
("FSSwingTimeStep",
"How long time steps is are to be used for FS colour reconnections.",
&EventFiller::theFSSwingTimeStep, 0.1, 0.0, 0,
true, false, Interface::lowerlim);
static Parameter<EventFiller,double> interfaceFSSwingTime
("FSSwingTime",
"How long time is allowed for FS colour reconnections. 0 turns it off.",
&EventFiller::theFSSwingTime, 0.0, 0.0, 0,
true, false, Interface::lowerlim);
static Switch<EventFiller,int> interfaceValenceChargeNormalisation
("ValenceChargeNormalisation",
"How to treat the (too large) colour charge of the valence partons.",
&EventFiller::theValenceChargeNormalisation, 0, true, false);
static SwitchOption interfaceValenceChargeNormalisationNone
(interfaceValenceChargeNormalisation,
"None",
"Dont do anything.",
0);
static SwitchOption interfaceValenceChargeNormalisationSwing
(interfaceValenceChargeNormalisation,
"Swing",
"Swing some of the dipoles going to the valence partons",
1);
static Parameter<EventFiller,Energy> interfacePTCut
("PTCut",
"The minimum invariant transverse momentum allowed for a gluon. "
"Gluons below the cut will be removed from the final state. "
"If zero, no gluons will be removed.",
&EventFiller::thePTCut, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
true, false, Interface::lowerlim);
static Switch<EventFiller,int> interfaceSoftRemove
("SoftRemove",
"Determines if gluons with invariant transverse momentum below "
"<interface>PTCut</interface> should be reabsorbed.",
&EventFiller::theSoftRemove, 1, true, false);
static SwitchOption interfaceSoftRemoveOff
(interfaceSoftRemove,
"Off",
"No gluons are absorbed",
0);
static SwitchOption interfaceSoftRemoveAll
(interfaceSoftRemove,
"All",
"All soft gluons below the cut are absorbed.",
1);
static SwitchOption interfaceSoftRemoveNoValence
(interfaceSoftRemove,
"NoValence",
"All except valence gluons are absorbed if below the cut.",
2);
static Switch<EventFiller,bool> interfaceOnlyOnce
("OnlyOnce",
"Do not allow a dipole to interact more than once.",
&EventFiller::onlyOnce, false, true, false);
static SwitchOption interfaceOnlyOnceManyInteractionsPerDipole
(interfaceOnlyOnce,
"ManyInteractionsPerDipole",
"Allow a dipole to interact several times.",
false);
static SwitchOption interfaceOnlyOnceOneInteractionPerDipole
(interfaceOnlyOnce,
"OneInteractionPerDipole",
"Only one interaction per dipole.",
true);
}
diff --git a/DIPSY/EventFiller.h b/DIPSY/EventFiller.h
--- a/DIPSY/EventFiller.h
+++ b/DIPSY/EventFiller.h
@@ -1,571 +1,572 @@
// -*- C++ -*-
#ifndef DIPSY_EventFiller_H
#define DIPSY_EventFiller_H
//
// This is the declaration of the EventFiller class.
//
#include "EventFiller.fh"
#include "ThePEG/Handlers/HandlerBase.h"
#include "DipoleEventHandler.fh"
#include "RealPartonState.fh"
#include "RealParton.fh"
#include "DipoleState.h"
#include "ImpactParameters.h"
#include "DipoleXSec.h"
#include "DipoleAbsorber.h"
#include "ThePEG/Utilities/Current.h"
#include "ThePEG/Analysis/FactoryBase.h"
#ifndef LWH_AIAnalysisFactory_H
#ifndef LWH
#define LWH ThePEGLWH
#endif
#include "ThePEG/Analysis/LWH/AnalysisFactory.h"
#endif
namespace DIPSY {
using namespace ThePEG;
/**
* The EventFiller class is able to produce an initial
* ThePEG::Collision from two colliding DipoleStates.
*
* @see \ref EventFillerInterfaces "The interfaces"
* defined for EventFiller.
*/
class EventFiller: public HandlerBase {
public:
/**
* Copy FList typedef from DipoleXSec.
*/
typedef DipoleXSec::FList FList;
/**
* Copy InteractionList typedef from DipoleXSec.
*/
typedef DipoleInteraction::List InteractionList;
/**
* Copy InteractionList typedef from DipoleXSec.
*/
typedef DipoleInteraction::PTSet InteractionPTSet;
/**
* A String is simply a vector of colour-connected partons.
*/
typedef DipoleState::String String;
/**
* A vector of dipole pairs represented as iterators into an FList.
*/
typedef DipoleXSec::DipolePairVector DipolePairVector;
/**
* An ordered map of dipole pairs represented as iterators into an
* FList.
*/
typedef DipoleXSec::DipolePairMap DipolePairMap;
/**
* Simple vector of interactions.
*/
typedef vector<DipoleInteraction::List::const_iterator> InteractionVector;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
EventFiller();
/**
* The destructor.
*/
virtual ~EventFiller();
//@}
public:
/** @name Virtual functions which can be overridden in subclasses. */
//@{
/**
* Fill the current collision object in the given DipoleEventHandler
* with the final state gluons produced when dipole states \a dl and
* \a dr collides with impact parameters \a b.
* @return the total interaction probability.
*/
virtual double fill(Step & step, DipoleEventHandler & eh, tPPair inc,
DipoleState & dl, DipoleState & dr,
const ImpactParameters & b) const;
/**
* Fill the current collision object in the given DipoleEventHandler
* with the final state gluons produced when dipole states \a dl and
* \a dr collides with impact parameters \a b.
* @return the total interaction probability.
*/
virtual double fill(Step & step, DipoleEventHandler & eh, tPPair inc,
const ImpactParameters & b,
DipoleState & dl, DipoleState & dr) const;
/**
* Fill the current collision object in the given DipoleEventHandler
* with the final state gluons produced when dipole states \a dl and
* \a dr collides with impact parameters \a b.
* @return the total interaction probability.
*/
virtual double fillWithShadows(Step & step, DipoleEventHandler & eh, tPPair inc,
const ImpactParameters & b,
DipoleState & dl, DipoleState & dr) const;
/**
* Select which dipole pairs should interact, and return them in the
* order in which they should be processed.
*/
virtual pair<RealPartonStatePtr, RealPartonStatePtr>
selectInteractions(const FList & fl, const ImpactParameters & b, const DipoleXSec & xSec) const;
/**
* Select which dipole pairs should interact, and return them in the
* order in which they should be processed.
*/
virtual pair<RealPartonStatePtr, RealPartonStatePtr>
selectInteractions(const InteractionList & intl, const DipoleXSec & xSec) const;
/**
* Select which dipole pairs should interact, and return them in the
* order in which they should be processed. (shadow parton version)
*/
virtual InteractionVector
selectShadowInteractions(const InteractionList & intl, const DipoleXSec & xSec) const;
/**
* Tries to do the interaction between the real states lrs and rrs with the two dipoles
* pointed to by inter, assuming that the interactions in inters already has been tested
* and accepted. States collide at impact paramter b and using the recoils in xSec.
*/
bool addInteraction(FList::const_iterator inter,
RealPartonStatePtr lrs, RealPartonStatePtr rrs,
DipolePairVector & inters,
const ImpactParameters & b, const DipoleXSec & xSec) const ;
/**
* Tries to do the interaction between the real states lrs and rrs with the two dipoles
* pointed to by inter, assuming that the interactions in inters already has been tested
* and accepted. States collide at impact paramter b and using the recoils in xSec.
*/
bool addInteraction(InteractionList::const_iterator inter,
RealPartonStatePtr lrs, RealPartonStatePtr rrs,
const TransverseMomentum & recoil,
InteractionVector & inters, const DipoleXSec & xSec) const ;
/**
* Go through all previously accepted interactions and check that
- * they can still be performed.
+ * they can still be performed. If \a final, this is the last check
+ * and partons are set on-shell.
*/
- bool recheckInteractions(const InteractionVector & interactions) const;
+ bool recheckInteractions(const InteractionVector & interactions, bool final) const;
/**
* Extract all strings after according to the given interactions.
*/
virtual vector<String>
extractStrings(DipoleState & dl, DipoleState & dr,
pair<RealPartonStatePtr, RealPartonStatePtr>,
const ImpactParameters & b) const;
/**
* Extract all strings after according to the given interactions.
*/
virtual vector<String>
extractShadowStrings(DipoleState & dl, DipoleState & dr,
const InteractionVector & interactions,
const ImpactParameters & b) const;
/**
* Fill the given Step with a SubProcess object using the given
* incoming particles and list of strings.
*/
virtual bool fillStep(Step & step, tPPair incoming,
const vector<String> & strings) const;
/**
* Fix up valens partons if they were not of the correct flavour or
* if they should collapse into a hadron.
*/
// virtual void fixValence(Step & step, DipoleState & dl, DipoleState & dr) const;
/**
* get the recoil scheme.
*/
inline int recoilScheme() const {
return theRecoilScheme;
}
/**
* set the recoil scheme.
*/
inline void recoilScheme(int x) {
theRecoilScheme = x;
}
/**
* get the mode.
*/
inline int mode() const {
return theMode;
}
/**
* set the mode.
*/
inline void mode(int x) {
theMode = x;
}
/**
* get the singleMother.
*/
inline int singleMother() const {
return theSingleMother;
}
/**
* set the singleMother.
*/
inline void singleMother(int x) {
theSingleMother = x;
}
/**
* get the DGLAPinPT.
*/
inline int DGLAPinPT() const {
return theDGLAPinPT;
}
/**
* set the DGLAPinPT.
*/
inline void DGLAPinPT(int x) {
theDGLAPinPT = x;
}
/**
* get the ValenceChargeNormalisation
**/
inline int valenceChargeNormalisation() const {
return theValenceChargeNormalisation;
}
/**
* The minimum squared invariant transverse momentum allowed for a
* gluon in a string.
*/
Energy2 pT2Cut() const {
return sqr(thePTCut);
}
/**
* Return the squared invariant transverse momentum of gluon \a i2
* in a string given its colour-neighbours \a i1 and \a i3.
*/
static Energy2 invPT2(const String & str, int i1, int i2, int i3);
/**
* Remove a gluon \a i2 in a string shuffling its momenum to its
* colour-neighbours \a i1 and \a i3.
*/
static bool removeGluon(const String & str, int i1, int i2, int i3);
/**
* get the effectiveWeights.
*/
inline int effectiveWeights() const {
return theEffectiveWeights;
}
/**
* get the FS swing time.
*/
inline double FSSwingTime() const {
return theFSSwingTime;
}
/**
* get the step size for the FS swing time.
*/
inline double FSSwingTimeStep() const {
return theFSSwingTimeStep;
}
//@}
public:
/**
* Get the object used to absorb non-interacting dipoles.
*/
inline DipoleAbsorberPtr absorber() const {
return theAbsorber;
}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* 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;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/**
* The weight which will be used for the event being generated.
*/
mutable double currentWeight;
mutable bool fail;
/**
* Debug/testing
*/
protected:
mutable double nInt;
mutable double nInt1;
mutable double nTried;
mutable double nEvents;
mutable double nTestedInts;
mutable double foundInt;
mutable double failedEvo;
mutable double failedRec;
mutable double rescatter;
mutable double overTenEvents;
mutable double rejectedEvents;
mutable double nLoops;
mutable double nYes;
mutable double nNo;
mutable double avYInEvo;
mutable double avYInInt;
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit() throw(InitException);
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Is this where this should be declared?? :o
*/
virtual void dofinish();
private:
protected:
/**
* balances p+/p- by moving the entire state in rapidity. Assumed to be
* UNMIRRORED states, ie, plus and minus are flipped for the right state.
**/
void fixBoost(RealPartonStatePtr lrs, RealPartonStatePtr rrs) const;
/**
* Checks through the final state and tries to fix the most blatant errors.
* Specially tries to catch bad momenta that can cause crashes.
**/
void dodgeErrors(DipoleStatePtr finalState) const;
/**
* Checks if the interacting partons in the current states have enough
* p+- to set the other state on shell.
*/
bool controlRecoils(DipolePairVector & sel,
RealPartonStatePtr rrs, RealPartonStatePtr lrs,
const ImpactParameters & b, const DipoleXSec & xSec,
pair<pair<bool, bool>, pair<bool, bool> > doesInt) const;
/**
* Checks if the interacting partons in the current states have enough
* p+- to set the other state on shell.
*/
bool controlRecoils(InteractionVector & sel,
RealPartonStatePtr rrs, RealPartonStatePtr lrs,
const DipoleXSec & xSec) const;
/**
* Removes the off shell partons and recouples colour flow. p_mu unchanged.
*/
void removeVirtuals(DipoleStatePtr state) const;
/**
* Removes the parton and recouples colour flow. p_mu unchanegd.
*/
void removeParton(tPartonPtr p) const;
/**
* find the pT scale through the dipolestates eventhandlers emitter object.
*/
double pTScale(DipoleState &) const;
/**
* Takes statistics on the number of participants
* in a heavy ion collision.
*/
void countParticipants(const DipoleState & dl, const DipoleState & dr, const InvEnergy b) const;
/**
* The object used to absorb non-interacting dipoles.
*/
DipoleAbsorberPtr theAbsorber;
/**
* What scheme to use when doing recoils.
*/
int theRecoilScheme;
/**
* In what mode to create the real state. Fast vs consistent.
*/
int theMode;
/**
* If partons have one or two mothers.
*/
int theSingleMother;
/**
* If DGLAP supression is made in terms of pt rather than r.
*/
int theDGLAPinPT;
/**
* How to distribute recoil among the members of an effective parton.
*/
int theEffectiveWeights;
/**
* How long time (in GeV-1) the FS get's to do colour reconnections.
*/
double theFSSwingTime;
/**
* How large time steps (in GeV-1) used in the FS colour
* reconnections.
*/
double theFSSwingTimeStep;
/**
* How the valence charge is handled.
**/
int theValenceChargeNormalisation;
/**
* The minimum invariant transverse momentum allowed for a gluon in a string.
*/
Energy thePTCut;
/**
* Options for removing gluons with too small invariant transverse
* momentum.
*/
int theSoftRemove;
/**
* Do not allow a dipole to interact more than once.
*/
bool onlyOnce;
// *** TO REMOVE *** Temporary analysis
mutable FactoryBase::tH1DPtr histptmax, histptmaxe, histptfi, histptf;
mutable double sumwmaxpt;
public:
/**
* Exception class for space-like gluon momenta.
*/
struct SpaceLikeGluons: public Exception {};
/**
* Exception class for space-like gluon momenta.
*/
struct ConsistencyException: public Exception {};
/**
* Exception class for failed gluon removal.
*/
struct RemoveGluonException: public Exception {};
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
EventFiller & operator=(const EventFiller &);
};
}
#endif /* DIPSY_EventFiller_H */
diff --git a/DIPSY/ShadowParton.cc b/DIPSY/ShadowParton.cc
--- a/DIPSY/ShadowParton.cc
+++ b/DIPSY/ShadowParton.cc
@@ -1,460 +1,501 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ShadowParton class.
//
#include "ShadowParton.h"
#include "Parton.h"
#include "DipoleState.h"
#include "DipoleEventHandler.h"
#include "Dipole.h"
#include "ImpactParameters.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "ThePEG/Utilities/Current.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/DescribeClass.h"
using namespace DIPSY;
ShadowParton::ShadowParton(Parton & p)
: theOriginal(&p), thePosition(p.position()), thePlus(p.plus()),
thePT(p.pT()), thePT0(p.pT()), theMass(p.mass()), theMinus(p.minus()),
theY(p.y()), theY0(p.y()), theFlavour(p.flavour()),
theDifference(first), theDipoles(p.dipoles()), theEmissionFactor(ZERO),
theRes(ZERO), hasInteracted(false), isOnShell(false), theInteraction(0),
- isInteracting(false), isValence(p.valence()) {}
+ isInteracting(0), isValence(p.valence()) {}
ShadowParton::~ShadowParton() {}
void ShadowParton::setupParton() {
if ( tSPartonPtr sp = original()->shadow() ) {
sp->theNext = this;
thePrevious = sp;
}
original()->shadow(this);
}
SPartonPtr ShadowParton::createValence(Parton & p) {
SPartonPtr sp = new_ptr(ShadowParton(p));
sp->setupParton();
sp->isValence = true;
return sp;
}
void ShadowParton::setupSwing(Parton & pa1, Parton & pa2, Parton & pb1, Parton & pb2) {
SPartonPtr sa1 = new_ptr(ShadowParton(pa1));
sa1->setupParton();
SPartonPtr sa2 = new_ptr(ShadowParton(pa2));
sa2->setupParton();
SPartonPtr sb1 = new_ptr(ShadowParton(pb1));
sb1->setupParton();
SPartonPtr sb2 = new_ptr(ShadowParton(pb2));
sb2->setupParton();
sa1->theDifference = sa2->theDifference = sb1->theDifference = sb2->theDifference = swing;
sa1->setupNeighbors();
sa2->setupNeighbors();
sb1->setupNeighbors();
sb2->setupNeighbors();
}
void ShadowParton::setupEmission(Parton & emitter, Parton & produced, Parton & recoiler) {
SPartonPtr se = new_ptr(ShadowParton(emitter));
se->setupParton();
SPartonPtr sp = new_ptr(ShadowParton(produced));
sp->setupParton();
SPartonPtr sr = new_ptr(ShadowParton(recoiler));
sr->setupParton();
se->theDifference = emission;
se->theSibling = sp;
se->setupNeighbors();
sp->theDifference = emission;
sp->theParent = this;
sp->theSibling = se;
theChild = sp;
sp->setupNeighbors();
sr->theDifference = recoil;
sr->setupNeighbors();
InvEnergy2 d2 = se->theRes = sp->theRes = se->dist2(*sp);
se->theEmissionFactor = sp->theEmissionFactor = alphaS(d2)*d2/UseRandom::rnd();
se->pT0(sp->pT0());
se->y0(sp->y0());
}
void ShadowParton::setupNeighbors() {
if ( dipoles().first )
theNeighbors.first = dipoles().first->partons().first->shadow();
if ( dipoles().second )
theNeighbors.second = dipoles().second->partons().second->shadow();
}
tSPartonPtr ShadowParton::last() {
if ( next() ) return next();
return this;
}
tSPartonPtr ShadowParton::initial() {
if ( previous() ) return previous();
return this;
}
tcSPartonPtr ShadowParton::last() const {
if ( next() ) return next();
return this;
}
double ShadowParton::alphaS(InvEnergy2 r2) {
return Current<DipoleEventHandler>()->alphaS(sqrt(r2));
}
tcSPartonPtr ShadowParton::resolve(InvEnergy2 r2, tPartonPtr stopp) const {
if ( resolved(r2, stopp) || !mother() ) return this;
return mother()->resolve(r2, stopp);
}
tSPartonPtr ShadowParton::resolve(InvEnergy2 r2, tPartonPtr stopp) {
if ( resolved(r2, stopp) || !mother() ) return this;
return mother()->resolve(r2, stopp);
}
void ShadowParton::flagOnShell(InvEnergy2 r2, tPartonPtr stopp) {
tSPartonPtr resolved = resolve(r2, stopp);
if ( resolved->onShell() ) return;
resolved->onShell(true);
if ( resolved->sibling() ) {
r2 = resolved->res();
resolved->sibling()->onShell(true);
}
if ( resolved->mother() ) resolved->mother()->flagOnShell(r2, stopp);
// *** TODO *** probably the whole function should be deleted
// else dipoleState().flagShadowsOnShell();
}
void ShadowParton::pTplus(const TransverseMomentum & qt, Energy qp) {
pT(qt);
plus(qp);
minus(mt2()/plus());
y(log(mt()/plus()));
}
void ShadowParton::setOnShell(bool always) {
original()->plus(plus());
- original()->minus(minus());
original()->pT(pT());
+ original()->minus(original()->mt2()/plus());
original()->y(log(original()->mt()/original()->plus()));
- if ( !interacted() || always ) original()->onShell(true);
+ if ( interaction() ) return;
+ if ( !interacted() || always ) {
+ original()->onShell(true);
+ onShell(true);
+ }
}
void ShadowParton::interact() {
hasInteracted = true;
if ( !interaction() && mother() ) mother()->interact();
}
void ShadowParton::resetInteracted() {
hasInteracted = false;
- interacting(false);
+ interacting(0);
if ( next() ) next()->resetInteracted();
if ( child() ) child()->resetInteracted();
}
void ShadowParton::reset() {
onShell(false);
interaction(0);
hasInteracted = false;
- interacting(false);
+ interacting(0);
if ( next() ) next()->reset();
if ( child() ) child()->reset();
}
void ShadowParton::insertInteraction(const DipoleInteraction * i) {
//if this shadow has already interacted we do nothing.
if ( interaction() ) return;
// We insert here if this is a valens or was emitted from an interaction chain.
if ( interacted() || !mother() || mother()->interacted() ) interaction(i);
// Otherwise we insert further up the chain.
else mother()->insertInteraction(i);
}
void ShadowParton::rejectInteraction(const DipoleInteraction * i) {
- if ( interaction() ) {
- if ( interaction() == i ) interaction(0);
- } else {
- if ( mother() ) mother()->rejectInteraction(i);
- }
+ if ( interaction() && interaction() == i ) interaction(0);
+ else if ( mother() ) mother()->rejectInteraction(i);
}
void ShadowParton::acceptInteraction(const DipoleInteraction * i) {
hasInteracted = true;
if ( !mother() ) dipoleState().acceptShadowInteraction(this);
else if ( interaction() != i ) mother()->acceptInteraction(i);
}
void ShadowParton::makeIncoming() {
if ( interacted() && onShell() ) onShell(false);
else if ( mother() ) mother()->makeIncoming();
}
bool ShadowParton::setEmissionMomenta(const LorentzMomentum & p, bool valencefix) {
if ( parent() ) {
// If this parton was emitted pretend it was actually he emitter
// and continues as a propagator. It will retain its original
// rapidity, but will get the transverse momentum of the incoming
// propagator. The exception is if we are putting a valence on
// shell in an otherwise unresolved emission - in that case the
// transverse momentum is given by the initial valence.
TransverseMomentum qT = pT0();
- if ( valencefix ) qT = -sibling()->initial()->pT();
+ if ( valencefix ) qT = -sibling()->initial()->pT()/2.0; // *** TODO *** Think this through!
pT(TransverseMomentum(p.x(), p.y()) + qT);
plus(mt0()*exp(y0()));
y(y0());
// The sibling will get minus the original transverse momentum and
// what ever is left of the positive light-cone momenta. (return
// false if not enough to put on shell).
sibling()->pT(-qT);
sibling()->plus(p.plus() - mt()*exp(-y()));
if ( sibling()->plus() < ZERO ) return false;
sibling()->minus(mt2()/sibling()->plus());
sibling()->y(log(sibling()->mt()/sibling()->plus()));
// The propagators negative light-cone momentum is trivial,
minus(p.minus() - sibling()->minus());
} else {
// If this parton was assumed in the evolution, life is simpler.
sibling()->pT(pT0());
sibling()->plus(sibling()->mt()*exp(-y0()));
sibling()->minus(sibling()->mt2()/sibling()->plus());
sibling()->y(y0());
pT(TransverseMomentum(p.x(), p.y()) - pT0());
plus(p.plus() - sibling()->plus());
if ( plus() < ZERO ) return false;
minus(mt2()/plus());
y(log(mt()/plus()));
}
return true;
}
ShadowParton::Propagator ShadowParton::
propagator(InvEnergy2 r2, tPartonPtr stopp, const DipoleInteraction * i, bool final) {
// If we reach a parton that has already been directly involved with
- // a previous interaction, simply return its momentum.
- if ( interacting() ) return momentum();
+ // a previous interaction, simply return its momentum. If it was
+ // flagged inteacting reset the on-shell flag as it is now
+ // rescattering.
+ if ( interacting() ) {
+ if ( final ) original()->onShell(false);
+ return momentum();
+ }
// If this was a valence, ask the DipoleState about its momentum.
if ( !mother() ) return dipoleState().incomingMomentum(this, final);
// Now, if this emission was not resolved then simply return the
// propagator of the mother. However, if the sibling is a valence
// that would not be set on shell, make the emission anyway, but
// with a special procedure.
- bool valencefix = false;
- if ( !resolved(r2, stopp) ) {
- if ( sibling() && sibling()->valence() ) valencefix = true;
- else return mother()->propagator(r2, stopp, i, final);
- }
+ bool emitvalence = ( sibling() && sibling()->valence() );
+ bool unresolved = !resolved(r2, stopp);
+
+ if ( unresolved && !emitvalence )
+ return mother()->propagator(r2, stopp, i, final);
- // Get the propagator of the mother with the new resolution scale.
- Propagator prop = mother()->propagator(res(), stopp, i, final);
+ // Get the propagator of the mother with the new resolution
+ // scale. But don't put anything on-shell yet in case the emission is rejected later.
+ Propagator prop = mother()->propagator(res(), stopp, i, false);
// If it is not possible to set the emitted parton on shell, flag
// the propagator as failed.
- if ( !setEmissionMomenta(prop.p, valencefix) ) {
+ if ( !setEmissionMomenta(prop.p, false) ) {
prop.fail = true;
// If the emitted parton is necessary for a subsequent
// interaction, return the failed propagator.
if ( sibling()->interaction() ) return prop;
// Otherwise just ignore the emission.
return mother()->propagator(r2, stopp, i, final);
}
// Now depending on the colour line of the emitted parton, it must
// be ordered in light-cone momenta with previous partons on the
- // same side except if it patrakes in a subsequent interaction.
+ // same side except if it partakes in a subsequent interaction.
if ( colourSibling() ) {
if ( sibling()->interaction() ) {
prop.colpos = Constants::MaxEnergy;
prop.colneg = ZERO;
} else {
+ // Check if a resolved emission was not ordered, in that case
+ // flag it unresolved.
if ( sibling()->plus() > prop.colpos ||
- sibling()->minus() < prop.colneg )
- return mother()->propagator(r2, stopp, i, final);
- if ( !valencefix ) { // Artificially emitted valence does not affect ordering.
+ sibling()->minus() < prop.colneg )
+ unresolved = true;
+ else {
prop.colpos = sibling()->plus();
prop.colneg = sibling()->minus();
}
}
} else {
- if ( sibling()->interaction() ) {
+ if ( sibling()->interaction() ) {
prop.acopos = Constants::MaxEnergy;
prop.aconeg = ZERO;
} else {
if ( sibling()->plus() > prop.acopos ||
sibling()->minus() < prop.aconeg )
- return mother()->propagator(r2, stopp, i, final);
- if ( !valencefix ) { // Artificially emitted valence does not affect ordering.
+ unresolved = true;
+ else {
prop.acopos = sibling()->plus();
prop.aconeg = sibling()->minus();
}
}
}
+ // If the emission was found to be unresolved we need to recalculate the incoming momentum.
+ if ( unresolved ) {
+ prop = mother()->propagator(r2, stopp, i, final);
+ // If the unresolved emission was a valence we need to emit it
+ // anyway with a special treatment.
+ if ( emitvalence ) {
+ if ( !setEmissionMomenta(prop.p, true) ) {
+ prop.fail = true;
+ return prop;
+ }
+ } else
+ return prop;
+ }
+ else if ( final )
+ // In any case we need to redo the incoming momentum to mark it final.
+ prop = mother()->propagator(res(), stopp, i, final);
+
// Now we are done. if this is the final check, then put emitted
// partons on shell (if not interacted) and return the propagator.
if ( final ) sibling()->setOnShell();
prop.p -= sibling()->momentum();
return prop;
}
tSPartonPtr ShadowParton::findFirstOnShell() {
if ( original()->onShell() ) return this;
if ( !next() ) return tSPartonPtr();
if ( next()->colourSibling() ) {
tSPartonPtr ch = child()->findFirstOnShell();
if ( ch ) return ch;
}
return next()->findFirstOnShell();
}
tSPartonPtr ShadowParton::findSecondOnShell() {
if ( original()->onShell() ) return this;
if ( !next() ) return tSPartonPtr();
if ( !next()->colourSibling() ) {
tSPartonPtr ch = child()->findSecondOnShell();
if ( ch ) return ch;
}
return next()->findSecondOnShell();
}
Ariadne5::ClonePtr ShadowParton::clone() const {
return new_ptr(*this);
}
void ShadowParton::mirror(double y0) {
y(2.0*y0 - y());
swap(thePlus, theMinus);
if ( next() ) next()->mirror(y0);
if ( child() ) child()->mirror(y0);
}
void ShadowParton::translate(const ImpactParameters & b) {
thePosition = b.translate(thePosition);
thePT = b.rotatePT(thePT);
if ( next() ) next()->translate(b);
if ( child() ) child()->translate(b);
}
void ShadowParton::propagateShadowMomenta() {
if ( onShell() ) {
original()->plus(plus());
original()->pT(pT());
original()->y(log(original()->mt()/original()->plus()));
original()->onShell(true);
}
if ( next() ) next()->propagateShadowMomenta();
if ( child() ) child()->propagateShadowMomenta();
}
void ShadowParton::rebind(const TranslationMap & trans) {
thePrevious = trans.translate(thePrevious);
theParent = trans.translate(theParent);
theSibling = trans.translate(theSibling);
theNext = trans.translate(theNext);
theChild = trans.translate(theChild);
theDipoles.first = trans.translate(theDipoles.first);
theDipoles.second = trans.translate(theDipoles.second);
theNeighbors.first = trans.translate(theNeighbors.first);
theNeighbors.second = trans.translate(theNeighbors.second);
}
DipoleState& ShadowParton::dipoleState() const {
if ( dipoles().first ) return dipoles().first->dipoleState();
else return dipoles().second->dipoleState();
}
double ShadowParton::pTScale() const {
return dipoleState().handler().emitter().pTScale();
}
Energy ShadowParton::mass() const {
if ( theMass < ZERO )
theMass = CurrentGenerator::current().getParticleData(flavour())->mass();
return theMass;
}
bool ShadowParton::partonOnShell() const {
return onShell() || ( previous() && previous()->partonOnShell() );
}
void ShadowParton::debugme() {
cout << "data for ShadowParton " << this << endl;
cout << "thePosition1: " << thePosition.first*GeV
<< ", thePosition2: " << thePosition.second*GeV
<< ", thePlus: " << thePlus/GeV
<< ", thePT: " << thePT.pt()/GeV
<< ", theMinus: " << theMinus/GeV
<< ", theY: " << theY
<< ", theFlavour: " << theFlavour
<< ", theDipoles1: " << theDipoles.first
<< ", theDipoles2: " << theDipoles.second
<< ", hasInteracted: " << hasInteracted
<< ", isOnShell: " << isOnShell
<< ", isValence: " << isValence
<< ", theMass: " << theMass/GeV << endl;
}
void ShadowParton::debugTree(string indent) {
if ( mother() && !child() && next() && !onShell() ) {
next()->debugTree(indent);
return;
}
- cerr << indent << "-"
- << (valence()? "v": "-")
- << (onShell()? "*": "-")
- << (original()->onShell()? "o": "-")
- << (interacted()? "x": "-") << "->";
- if ( onShell() )
- cerr << "[" << y() << ", " << pt()/GeV << "]";
+ cerr << indent << "--"
+ << (valence()? "v": "-");
+ if ( interaction() ) cerr << interaction()->id;
+ cerr << (interacted()? "+": "-");
+ if ( interacting() ) cerr << interacting();
+ cerr << (onShell()? "*": "-") << "->";
+ if ( original()->onShell() )
+ cerr << "[" << original()->pT().x()/GeV << ", "
+ << original()->pT().y()/GeV << ", "
+ << original()->y() << "]";
cerr << endl;
- if ( next() ) next()->debugTree(indent + " ");
- if ( child() ) child()->debugTree(indent + " ");
+ if ( indent.size() > 1 && indent[indent.size() - 1] == '\\' )
+ indent[indent.size() - 1] = ' ';
+ if ( next() && child() ) next()->debugTree(indent + " |");
+ if ( next() && !child() ) next()->debugTree(indent + " ");
+ if ( child() ) child()->debugTree(indent + " \\");
}
+/*
+----->
+ |---->
+ | |---->
+ | \---->
+ | |---->
+ | \---->
+ \---->
+ */
void ShadowParton::persistentOutput(PersistentOStream & os) const {
os << theOriginal << ounit(thePosition, InvGeV) << ounit(thePlus, GeV)
<< ounit(thePT, GeV) << ounit(theMass, GeV) << ounit(theMinus, GeV)
<< theY << theFlavour << thePrevious << theParent << theSibling
<< theNext << theChild << theDipoles << theNeighbors
<< ounit(theEmissionFactor, 1.0/GeV2) << ounit(theRes, 1.0/GeV2)
<< hasInteracted << isOnShell << isValence << oenum(theDifference);
}
void ShadowParton::persistentInput(PersistentIStream & is, int) {
is >> theOriginal >> iunit(thePosition, InvGeV) >> iunit(thePlus, GeV)
>> iunit(thePT, GeV) >> iunit(theMass, GeV)>> iunit(theMinus, GeV)
>> theY >> theFlavour >> thePrevious >> theParent >> theSibling
>> theNext >> theChild >> theDipoles >> theNeighbors
>> iunit(theEmissionFactor, 1.0/GeV2) >> iunit(theRes, 1.0/GeV2)
>> hasInteracted >> isOnShell >> isValence >> ienum(theDifference);
}
DescribeClass<ShadowParton,Ariadne5::CloneBase>
describeDIPSYShadowParton("DIPSY::ShadowParton", "libAriadne5.so libDIPSY.so");
// Definition of the static class description member.
void ShadowParton::Init() {}
diff --git a/DIPSY/ShadowParton.h b/DIPSY/ShadowParton.h
--- a/DIPSY/ShadowParton.h
+++ b/DIPSY/ShadowParton.h
@@ -1,895 +1,895 @@
// -*- C++ -*-
#ifndef DIPSY_ShadowParton_H
#define DIPSY_ShadowParton_H
//
// This is the declaration of the ShadowParton class.
//
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Vectors/Transverse.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Repository/UseRandom.h"
#include "Ariadne/Config/CloneBase.h"
#include "ShadowParton.fh"
#include "Parton.fh"
#include "Dipole.fh"
#include "DipoleXSec.fh"
namespace DIPSY {
using namespace ThePEG;
class ImpactParameters;
/**
* Here is the documentation of the ShadowParton class.
*/
class ShadowParton: public Ariadne5::CloneBase {
public:
/**
* Typedef for position in transverse coordinate space.
*/
typedef Transverse<InvEnergy> Point;
/**
* A pair of partons.
*/
typedef pair<tSPartonPtr,tSPartonPtr> tSPartonPair;
/**
* A pair of dipoles.
*/
typedef pair<tDipolePtr,tDipolePtr> tDipolePair;
/**
* The Dipole is a friend.
*/
friend class Dipole;
/**
* Enum difference from the previous version of the same parton.
*/
enum Difference {
first = 0, /**< No difference - this is the first instance. */
emission = 1, /**< The previous instance emitted. */
recoil = 2, /**< The previous received a recoil from a neighboring emission */
swing = 3, /**< The previous underwent a swing. */
interactn = 4 /**< The previous underwent an interaction */
};
/**
* Helper class for traversing along a propagator to an interaction
*/
struct Propagator {
/**
* Constructor.
*/
Propagator(): colpos(Constants::MaxEnergy), colneg(ZERO),
acopos(Constants::MaxEnergy), aconeg(ZERO), fail(false) {}
/**
* Construct from momentum.
*/
Propagator(const LorentzMomentum & mom)
: p(mom), colpos(Constants::MaxEnergy), colneg(ZERO),
acopos(Constants::MaxEnergy), aconeg(ZERO), fail(false) {}
/**
* The Momentum of the propagator.
*/
LorentzMomentum p;
/**
* The positive light-cone momentum of the previous emission on the colour line.
*/
Energy colpos;
/**
* Thenegative light-cone momentum of the previous emission on the colour line.
*/
Energy colneg;
/**
* The positive light-cone momentum of the previous emission on
* the anti-colour line.
*/
Energy acopos;
/**
* Thenegative light-cone momentum of the previous emission on the
* anti-colour line.
*/
Energy aconeg;
/**
* Indicate that the propagator has failed.
*/
bool fail;
};
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
inline ShadowParton()
: thePosition(Point()), thePlus(ZERO), theMass(-1*GeV), theMinus(0*GeV),
theY(0), theY0(0), theFlavour(ParticleID::g), theDifference(first),
theEmissionFactor(ZERO), theRes(ZERO), hasInteracted(false), isOnShell(false),
- theInteraction(0), isInteracting(false), isValence(false) {}
+ theInteraction(0), isInteracting(0), isValence(false) {}
/**
* Construct from ordinary parton.
*/
ShadowParton(Parton & p);
/**
* The destructor.
*/
virtual ~ShadowParton();
/**
* Create a valence shadow parton from a given ordinary parton.
*/
static SPartonPtr createValence(Parton & p);
/**
* Easy access to alphaS.
*/
static double alphaS(InvEnergy2 r2);
/**
* Setup shadow partons in an emission.
*/
void setupEmission(Parton & emitter, Parton & produced, Parton & recoiler);
/**
* Setup shadow partons in a swing.
*/
static void setupSwing(Parton & pa1, Parton & pa2, Parton & pb1, Parton & pb2);
/**
* Setup pointers to and from neighboring partons.
*/
void setupParton();
/**
* Setup neighboring shadow partons.
*/
void setupNeighbors();
//@}
protected:
/** @name The virtual functions to be overridden in sub-classes. */
//@{
/**
* Return a simple clone of this object. Should be implemented as
* <code>return new_ptr(*this);</code> by a derived class.
*/
virtual Ariadne5::ClonePtr clone() const;
/**
* Rebind pointers to other CloneBase objects. Called after a number
* of interconnected CloneBase objects have been cloned, so that
* the cloned objects will refer to the cloned copies afterwards.
*
* @param trans a TranslationMap relating the original objects to
* their respective clones.
*/
virtual void rebind(const TranslationMap & trans);
//@}
public:
/**
* Return the DipoleState to which this parton belongs. If it was
* originally belonged to one state which was subsequently merged
* into another, the parton will belong to the latter state.
*/
DipoleState & dipoleState() const;
/**
* finds the pTscale through the dipolestates emitter.
*/
double pTScale() const;
/**
* Calculate the squared transverse distance to the given parton.
*/
inline InvEnergy2 dist2(const ShadowParton & p) const {
return (position() - p.position()).pt2();
}
/**
* Calculate the transverse distance to the given parton.
*/
inline InvEnergy dist(const ShadowParton & p) const {
return sqrt(dist2(p));
}
/**
* Produce a ThePEG::Particle corresponding to this parton.
*/
PPtr produceParticle() const;
/**
* The mass of this parton.
*/
Energy mass() const;
/**
* The transverse momentum squared of this parton.
*/
Energy2 pt2() const {
return pT().pt2();
}
/**
* The transverse mass squared of this parton.
*/
Energy2 mt2() const {
return pt2() + sqr(mass());
}
/**
* The transverse mass of this parton in the emission.
*/
Energy2 mt02() const {
return pT0().pt2() + sqr(mass());
}
/**
* The transverse mass of this parton in the emission.
*/
Energy mt0() const {
return sqrt(mt02());
}
/**
* The transverse momentum of this parton.
*/
Energy pt() const {
return sqrt(pt2());
}
/**
* The transverse momentum of this parton.
*/
Energy mt() const {
return sqrt(mt2());
}
/**
* The final-state momentum of this particle.
*/
LorentzMomentum momentum() const {
- return lightCone(plus(), minus(), pT());
+ return lightCone(plus(), mt2()/plus(), pT());
}
/** @name Simple access functions. */
//@{
/**
* Get the position in impact parameter space.
*/
inline const Point & position() const {
return thePosition;
}
/**
* Get the positive light-cone momentum.
*/
inline Energy plus() const {
return thePlus;
}
/**
* Get the transverse momentum.
*/
inline TransverseMomentum pT() const {
return thePT;
}
/**
* Get the transverse momentum.
*/
inline TransverseMomentum pT0() const {
return thePT0;
}
/**
* Get the accumulated negative light-cone momentum deficit.
*/
inline Energy minus() const {
return theMinus;
}
/**
* Get the rapidity.
*/
inline double y() const {
return theY;
}
/**
* Get the rapidity generated in the emission.
*/
inline double y0() const {
return theY0;
}
/**
* Get the flavour of this parton.
*/
inline long flavour() const {
return theFlavour;
}
/**
* Indicate what happened to the previous instance of this parton.
*/
Difference diff() const {
return theDifference;
}
/**
* Get the original DIPSY parton.
*/
tPartonPtr original() const {
return theOriginal;
}
/**
* Get the previous instance of this parton befor it was emitted or
* swinged.
*/
inline tSPartonPtr previous() const {
return thePrevious;
}
/**
* Get the parent parton if this was emitted.
*/
inline tSPartonPtr parent() const {
return theParent;
}
/**
* Get the sibling of this parton if any.
*/
inline tSPartonPtr sibling() const {
return theSibling;
}
/**
* Get the next version of this parton if emitted or swinged.
*/
inline tSPartonPtr next() const {
return theNext;
}
/**
* Get the first version of this parton if emitted or swinged.
*/
tSPartonPtr initial();
/**
* Get the last version of this parton if emitted or swinged.
*/
tSPartonPtr last();
/**
* Get the last version of this parton if emitted or swinged.
*/
tcSPartonPtr last() const;
/**
* Get the parton this has emitted if any.
*/
inline tSPartonPtr child() const {
return theChild;
}
/**
* Get the mother parton, either previous() or parent().
*/
inline tSPartonPtr mother() const {
return parent()? parent(): previous();
}
/**
* Get the connecting dipoles.
*/
inline tDipolePair dipoles() const {
return theDipoles;
}
/**
* Indicate if this parton has interacted.
*/
inline bool interacted() const {
return hasInteracted;
}
/**
* Return true if the sibling is on the colour side of this parton.
*/
inline bool colourSibling() const {
return sibling() == theNeighbors.first;
}
/**
* Return if this shadow was the root of an interacting propagator.
*/
inline const DipoleInteraction * interaction() const {
return theInteraction;
}
/**
* Return if this shadow was directly involved in the interaction.
*/
- inline bool interacting() const {
+ inline int interacting() const {
return isInteracting;
}
/**
* Return if this parton is on shell.
*/
inline bool onShell() const {
return isOnShell;
}
/**
* Return true if the actual parton should be put on-shell;
*/
bool partonOnShell() const;
/**
* Set if this parton is on shell.
*/
inline void onShell(bool b) {
isOnShell = b;
}
/**
* Set if this parton is on shell.
*/
inline void interaction(const DipoleInteraction * i) {
theInteraction = i;
}
/**
* Set if this shadow is directly involved in an interaction.
*/
- inline void interacting(bool b) {
- isInteracting = b;
+ inline void interacting(int i) {
+ isInteracting = i;
}
/**
* Return if this parton is a valence parton.
*/
inline bool valence() const {
return isValence;
}
/**
* Set if this parton is a valence parton.
*/
inline void valence(bool b) {
isValence = b;
}
/**
* Set the position in impact parameter space.
*/
inline void position(const Point & x) {
thePosition = x;
}
/**
* Set the positive light-cone momentum.
*/
inline void plus(Energy x) {
thePlus = x;
}
/**
* Set the transverse momentum.
*/
inline void pT(const TransverseMomentum & x) {
thePT = x;
}
/**
* Set the momentum using transverse momentum and positive
* light-cone momentum. Set also other components assuming real.
*/
void pTplus(const TransverseMomentum & qt, Energy qp);
/**
* Set the transverse momentum.
*/
inline void pT0(const TransverseMomentum & x) {
thePT0 = x;
}
/**
* Set the accumulated negative light-cone momentum deficit.
*/
inline void minus(Energy x) {
theMinus = x;
}
/**
* Set the rapidity.
*/
inline void y(double x) {
theY = x;
}
/**
* Set the rapidity.
*/
inline void y0(double x) {
theY0 = x;
}
/**
* Set the flavour of this parton.
*/
inline void flavour(long x) {
theFlavour = (x == 0? ParticleID::g: x);
}
/**
* Set the parent parton.
*/
inline void parent(tSPartonPtr p) {
theParent = p;
}
/**
* Set the connecting dipoles.
*/
inline void dipoles(tDipolePair x) {
theDipoles = x;
}
/**
* Indicate that this parton has interacted and mark all parent and
* original partons if needed. Also mark all resolved emissions on-shell.
*/
void interact();
/**
* The emission factor for this parton if emitted or emitter.
*/
InvEnergy2 emissionFactor() const {
return theEmissionFactor;
}
/**
* The resolusion scale factor for this parton if emitted or emitter.
*/
InvEnergy2 res() const {
return theRes;
}
/**
* Return true if this parton and its sibling are resolved. An
* emission involving \a stopp is always resolved.
*/
bool resolved(InvEnergy2 r2, tPartonPtr stopp) const {
return ( ( sibling() && ( sibling()->original() == stopp
|| sibling()->interaction()
|| emissionFactor() > alphaS(r2)*r2 ) )
|| interacting() );
}
/**
* Go back in the history of this shadow parton and find the
* original parton which would be resolved at a given size. Always
* resolve emissions involving \a stopp.
*/
tSPartonPtr resolve(InvEnergy2 r2, tPartonPtr stopp);
/**
* Go back in the history of this shadow parton and find the
* original parton which would be resolved at a given size. Always
* resolve emissions involving \a stopp. (const version).
*/
tcSPartonPtr resolve(InvEnergy2 r2, tPartonPtr stopp) const;
/**
* Go back in the history and flag all partons which are resoved by
* the given size to be set on-shell. Always
* resolve emissions involving \a stopp.
*/
void flagOnShell(InvEnergy2 r2, tPartonPtr stopp);
/**
* If this parton is not on shell, find an emitted parton
* with the same colour lines (modulo swing) that is.
*/
tSPartonPtr findFirstOnShell();
/**
* If this parton is not on shell, find an emitted parton
* with the same colour lines (modulo swing) that is.
*/
tSPartonPtr findSecondOnShell();
/**
* Copy the momentum to the original parton and flag on-shell if it
* has not interacted.
*/
void setOnShell(bool always = false);
/**
* Go forward in the evolution and reset all interaction flags.
*/
void reset();
/**
* Go forward in the evolution and reset all interacted flags.
*/
void resetInteracted();
/**
* Set the momenta in an emission, assuming the incoming momentum \a p.
*/
bool setEmissionMomenta(const LorentzMomentum & p, bool valencefix);
/**
* Go back in the history and find the momentum of this incoming
* parton. Always resolve emissions involving \a stopp. If \a final
* is true, propagate the momentum of outgoing partons to the
* original partons and flag them on-shell if non-interacting.
*/
Propagator propagator(InvEnergy2 r2, tPartonPtr stopp,
const DipoleInteraction * i, bool final = false);
/**
* Recursively tell all this and any offspring to propagate their
* momentum to the original partons if set on shell.
*/
void propagateShadowMomenta();
/**
* Indicate that this parton should be tested for an
* interaction. Insert the corresponding interaction object at the
* correct place in the chain. The primary interaction is always at
* one of the incoming shadows. Secondary interactions are inserted
* at a parton which has been or should be made real by a previous
* emission. If this parton already had an interaction, simply leave
* it as it is.
*/
void insertInteraction(const DipoleInteraction * i);
/**
* Indicate that a previously tested interaction is accepted. This
* is done by marking all propagators interacted() in the chain down
* to the one where the interaction was found.
*/
void acceptInteraction(const DipoleInteraction * i);
/**
* Make this and all its anceseters incoming.
*/
void makeIncoming();
/**
* Indicate that the interaction that has been tested was rejected
* by simply removing the interaction inserted by acceptInteraction().
*/
void rejectInteraction(const DipoleInteraction * i);
/**
* Change direction before producing collision.
*/
void mirror(double yframe);
/**
* Translate according to the impagt parameter.
*/
void translate(const ImpactParameters & b);
/**
* Prints all the member variables to cerr.
*/
void debugme();
/**
* Debug the structure of the shadow tree.
*/
void debugTree(string indent);
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* 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();
private:
/**
* The original DIPSY parton
*/
tPartonPtr theOriginal;
/**
* The position in impact parameter space.
*/
Point thePosition;
/**
* The positive light-cone momentum.
*/
Energy thePlus;
/**
* The transverse momentum.
*/
TransverseMomentum thePT;
/**
* The transverse momentum generated in the emission.
*/
TransverseMomentum thePT0;
/**
* The mass of this parton.
*/
mutable Energy theMass;
/**
* The accumulated negative light-cone momentum deficit.
*/
Energy theMinus;
/**
* The rapidity.
*/
double theY;
/**
* The rapidity generated in the emission.
*/
double theY0;
/**
* The flavour of this parton.
*/
long theFlavour;
/**
* Indicate what happened to the previous instance of this parton.
*/
Difference theDifference;
/**
* The state of this parton, if any, before the swing or the emission of this
* parton's sibling. Is aways null if theParent exists.
*/
tSPartonPtr thePrevious;
/**
* The parent parton if this has emitted or swinged. Is aways null if
* thePrevious exists.
*/
tSPartonPtr theParent;
/**
* The sibling of this parton if any.
*/
tSPartonPtr theSibling;
/**
* The state of this parton after the emission of the child. Is null
* if this parton has not emitted.
*/
SPartonPtr theNext;
/**
* The child this parton has emitted, if any.
*/
SPartonPtr theChild;
/**
* The dipoles connected to this parton.
*/
tDipolePair theDipoles;
/**
* The partons connected to this.
*/
tSPartonPair theNeighbors;
/**
* The emission factor for this parton if emitted or emitter.
*/
InvEnergy2 theEmissionFactor;
/**
* The resolution factor for this parton if emitted or emitter.
*/
InvEnergy2 theRes;
/**
* Indicate if this parton has interacted.
*/
bool hasInteracted;
/**
* Indicate if this parton is on shell, that is if it has
* been supplied with the neccesary p+ or p- for left or rightmoving
* particles respectively.
*/
bool isOnShell;
/**
* This is root of an interaction chain.
*/
const DipoleInteraction * theInteraction;
/**
* Indicate that this shadow was directly involved with an interaction.
*/
- bool isInteracting;
+ int isInteracting;
/**
* Indicate if this parton is a valence parton.
*/
bool isValence;
protected:
/**
* Exception class for bad kinematics.
*/
struct ShadowPartonKinematicsException: public Exception {};
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ShadowParton & operator=(const ShadowParton &);
};
}
#endif /* DIPSY_ShadowParton_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:49 AM (1 d, 8 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983146
Default Alt Text
(254 KB)

Event Timeline