Page MenuHomeHEPForge
Files F11221236

No OneTemporary
Actions

View Options

diff --git a/DIPSY/DipoleState.cc b/DIPSY/DipoleState.cc
--- a/DIPSY/DipoleState.cc
+++ b/DIPSY/DipoleState.cc
@@ -1,639 +1,639 @@
// -*- 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/Utilities/Current.h"
#include "ThePEG/Repository/UseRandom.h"
#include <iostream>
#include <fstream>
using namespace DIPSY;
DipoleState::DipoleState(const DipoleState & x)
: theHandler(x.theHandler), theInitialDipoles(x.theInitialDipoles),
theWFInfo(x.theWFInfo), theWeight(x.theWeight), allDipoles(x.allDipoles) {
for ( set<DipolePtr>::iterator it = allDipoles.begin();
it != allDipoles.end(); ++it )
(**it).theDipoleState = this;
}
DipoleState::~DipoleState() {}
void DipoleState::printDipoles(ostream & output){
theSwingCandidates.clear();
theSwingCandidates.resize(Current<DipoleEventHandler>()->nColours());
for ( int i = 0, N = initialDipoles().size(); i < N; ++i ) {
sortDipole(*(initialDipoles()[i]));
}
for(int i=0, N=theSwingCandidates.size(); i < N;i++){
for(int j=0, N=theSwingCandidates[i].size(); j < N;j++){
output << i << '\t' << theSwingCandidates[i][j]->partons().first->position().first*GeV << '\t' <<
theSwingCandidates[i][j]->partons().first->position().second*GeV << '\t' <<0<< endl;
output << i << '\t' <<theSwingCandidates[i][j]->partons().second->position().first*GeV << '\t' <<
theSwingCandidates[i][j]->partons().second->position().second*GeV << '\t' <<0<<endl<<endl;
output << i << '\t' << theSwingCandidates[i][j]->partons().first->position().first*GeV << '\t' <<
theSwingCandidates[i][j]->partons().first->position().second*GeV << '\t' << 1 << endl;
output << i << '\t' << theSwingCandidates[i][j]->partons().first->position().first*GeV << '\t' <<
theSwingCandidates[i][j]->partons().first->position().second*GeV << '\t' << 2 << '\t'<<
theSwingCandidates[i][j]->partons().first->pT().first/GeV << '\t' <<
theSwingCandidates[i][j]->partons().first->pT().second/GeV << '\t'<<endl<< endl;
output << i << '\t' << theSwingCandidates[i][j]->partons().second->position().first*GeV <<'\t'<<
theSwingCandidates[i][j]->partons().second->position().second*GeV << '\t' << 1 << endl;
output << i << '\t' << theSwingCandidates[i][j]->partons().second->position().first*GeV <<'\t'<<
theSwingCandidates[i][j]->partons().second->position().second*GeV << '\t' << 2 << '\t'<<
theSwingCandidates[i][j]->partons().second->pT().first/GeV << '\t' <<
theSwingCandidates[i][j]->partons().second->pT().second/GeV << '\t'<<endl<< endl;
}
}
}
void DipoleState::sortDipoles() {
theSwingCandidates.clear();
theSwingCandidates.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
theSwingCandidates[d.colour()].push_back(& d);
}
void DipoleState::evolve(double ymin, double ymax) {
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();
continue;
}
break;
}
}
void DipoleState::generateColourIndex(tDipolePtr d) {
do {
d->colour(UseRandom::irnd(handler().nColours()));
} 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()) )
partons.push_back((*it)->partons().first);
if ( !((*it)->partons().second->interacted()) )
partons.push_back((*it)->partons().second);
}
return partons;
}
bool DipoleState::diagnosis() {
bool ok = true;
int activeDipoles = 0;
int partons = 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;
for(set< DipolePtr >::iterator it = allDipoles.begin();
it != allDipoles.end(); it++) {
if( ((*it)->children().first || (*it)->children().second) )
continue;
activeDipoles++;
PartonPtr p = (*it)->partons().first;
partons++;
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++;
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( abs((plus + plusNeeded - plusMissing) - (minus + minusNeeded - minusMissing)) > 1*GeV )
ok = false;
if ( !ok )
cout <<
" 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 <<
"allDipoles has " << allDipoles.size() << " dipoles." << endl <<
"number of active dipoles: " << activeDipoles << endl <<
"number of partons " << partons << endl <<
"hasQuark: " << hasQuark << endl;
return ok;
}
void DipoleState::reabsorb(PartonPtr p) {
DipolePtr d1 = p->dipoles().first;
DipolePtr d2 = p->dipoles().second;
PartonPtr p1 = d1->partons().first;
PartonPtr p2 = d2->partons().second;
if ( p1 == p2 ) {
// cout << "trying to reabsorb a loop..." << endl;
sortDipoles();
if( !forceSwing(p) ) return; //return if no swing was found. do not reabsorb.
d1 = p->dipoles().first;
d2 = p->dipoles().second;
p1 = d1->partons().first;
p2 = d2->partons().second;
}
DipolePtr d3 = createDipole();
d1->children().second = d3;
d2->children().first = d3;
p1->dipoles(make_pair(p1->dipoles().first,d3));
p2->dipoles(make_pair(d3,p2->dipoles().second));
d3->partons(make_pair(p1,p2));
d3->firstNeighbor(d1->neighbors().first);
d3->secondNeighbor(d2->neighbors().second);
d3->neighbors().first->secondNeighbor(d3);
d3->neighbors().second->firstNeighbor(d3);
bool right = p->rightMoving();
bool right1 = p1->rightMoving();
bool right2 = p2->rightMoving();
//remember how much p+/p- is supplied:
Energy suppliedMinus1 = p1->pT().pt()*exp( p1->y() ) - p1->minus();
Energy suppliedPlus1 = p1->pT().pt()*exp( -p1->y() ) - p1->plus();
Energy suppliedMinus2 = p2->pT().pt()*exp( p2->y() ) - p2->minus();
Energy suppliedPlus2 = p2->pT().pt()*exp( -p2->y() ) - p2->plus();
//Conserve pT and p+ as if reverse emission. y and p- adapt.
Parton::Point r1 = p->position() - p1->position();
Parton::Point r2 = p->position() - p2->position();
double P1 = r2.pt()/( r1.pt() + r2.pt() );
double P2 = r1.pt()/( r1.pt() + r2.pt() );
bool projection = false; // project out pT along the dipoles.
// conserves angular momentum.
bool justDistance = true; // distribute pT according to distance only.
// Bounds rapidity shifts.
if ( projection ) {
TransverseMomentum e1 = r1/sqr(r1.pt());
TransverseMomentum e2 = r2/sqr(r2.pt());
p2->pT(p2->pT() + e2*(p->pT().first*e1.second - p->pT().second*e1.first)/
(e2.first*e1.second-e2.second*e1.first));
p1->pT(p1->pT() + e1*(p->pT().first*e2.second - p->pT().second*e2.first)/
(e1.first*e2.second - e1.second*e2.first));
}
else if ( justDistance ) {
p1->pT( p1->pT() + P1*p->pT() );
p2->pT( p2->pT() + P2*p->pT() );
}
//for rightmoving partons, plus is conserved and minus is just measuring missing
//momentum. other way around for leftmoving.
if( right ) {
Energy supPm = p->pT().pt()*exp( p->y() ) - p->minus();
if(right1) {
p1->plus( p1->plus() + P1*p->plus() );
p1->y( log( p1->pT().pt()/p1->plus() ) );
p1->minus( p1->pT().pt()*exp( p1->y() ) - suppliedMinus1 - P1*supPm);
}
else {
p1->minus( p1->minus() + P1*supPm );
p1->y( log( p1->minus()/p1->pT().pt() ) );
p1->plus( p1->pT().pt()*exp( -p1->y() ) - suppliedPlus1 - P1*p->plus() );
}
if(right2) {
p2->plus( p2->plus() + P2*p->plus() );
p2->y( log( p2->pT().pt()/p2->plus() ) );
p2->minus( p2->pT().pt()*exp( p2->y() ) - suppliedMinus2 - P2*supPm );
}
else {
p2->minus( p2->minus() + P2*supPm );
p2->y( log( p2->minus()/p2->pT().pt() ) );
p2->plus( p2->pT().pt()*exp( -p2->y() ) - suppliedPlus2 - P2*p->plus() );
}
}
else {
Energy supPp = p->pT().pt()*exp( -p->y() ) - p->plus();
if(!right1) {
p1->minus( p1->minus() + P1*p->minus() );
p1->y( -log( p1->pT().pt()/p1->minus() ) );
p1->plus( p1->pT().pt()*exp( -p1->y() ) - suppliedPlus1 - P1*supPp );
}
else {
p1->plus( p1->plus() + P1*supPp );
p1->y( -log( p1->plus()/p1->pT().pt() ) );
p1->minus( p1->pT().pt()*exp( p1->y() ) - suppliedMinus1 - P1*p->minus() );
}
if(!right2) {
p2->minus( p2->minus() + P2*p->minus() );
p2->y( -log( p2->pT().pt()/p2->minus() ) );
p2->plus( p2->pT().pt()*exp( -p2->y() ) - suppliedPlus2 - P2*supPp );
}
else {
p2->plus( p2->plus() + P2*supPp );
p2->y( -log( p2->plus()/p2->pT().pt() ) );
p2->minus( p2->pT().pt()*exp( p2->y() ) - suppliedMinus2 - P2*p->minus() );
}
}
if(d1->size()/(d1->size()+d2->size()) > UseRandom::rnd() )
d3->colour(d1->colour());
else
d3->colour(d2->colour());
}
bool DipoleState::forceSwing(PartonPtr p) {
p->dipoles().first->reset();
bool foundFirstDipole = p->dipoles().first->forceGenerateRec();
p->dipoles().second->reset();
bool foundSecondDipole = p->dipoles().second->forceGenerateRec();
if( foundFirstDipole || foundSecondDipole ) {
if ( !foundFirstDipole )
p->dipoles().second->recombine();
else if ( !foundSecondDipole )
p->dipoles().first->recombine();
else {
double y1 = p->dipoles().first->generatedY();
double y2 = p->dipoles().second->generatedY();
double ymax = max(y1,y2);
if ( p->dipoles().first->swingDipole()->colour() != p->dipoles().first->colour() ) y1 += ymax;
if ( p->dipoles().second->swingDipole()->colour() != p->dipoles().second->colour() ) y2 += ymax;
if ( y1 < y2 )
p->dipoles().first->recombine();
else
p->dipoles().second->recombine();
}
return true;
}
else {
- cout << "couldn't find dipole to swing with! forced to leave it!!" << endl;
+// cout << "couldn't find dipole to swing with! forced to leave it!!" << endl;
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());
allDipoles.insert(otherState->allDipoles.begin(),otherState->allDipoles.end());
return this;
}
DipoleStatePtr DipoleState::collide( DipoleStatePtr otherState,
const vector<FList::const_iterator> & sel,
const ImpactParameters & b ) {
//mirror the other state in y
otherState->mirror();
//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() ) { //P_el/P_tot = f_ij/2 correct?
(*it)->second.first->swingDipole( (*it)->second.second );
(*it)->second.first->recombine();
//transfer pT!!
}
}
return this;
}
void DipoleState::mirror() {
for(set< DipolePtr >::iterator it = allDipoles.begin();
it != allDipoles.end(); it++) {
if( !((*it)->children().first || (*it)->children().second) ) {
PartonPtr p = (*it)->partons().first;
p->y( -p->y() );
p->plus( p->pT().pt()*exp(-p->y()) );
p->minus( p->pT().pt()*exp(p->y()) );
p->rightMoving( !(p->rightMoving()) );
if( !((*it)->neighbors().second) ) {
PartonPtr p2 = (*it)->partons().second;
p2->y( -p2->y() );
p2->plus( p2->pT().pt()*exp(-p2->y()) );
p2->minus( p2->pT().pt()*exp(p2->y()) );
p->rightMoving( !(p->rightMoving()) );
}
}
}
}
void DipoleState::translate(const ImpactParameters & b) {
for(set< DipolePtr >::iterator it = allDipoles.begin();
it != allDipoles.end(); it++) {
if( !((*it)->children().first || (*it)->children().second) ) {
(*it)->partons().first->position( Parton::Point
((*it)->partons().first->position().x()*cos(b.phi()) +
(*it)->partons().first->position().y()*sin(b.phi()) + b.bVec().x() ,
(*it)->partons().first->position().y()*cos(b.phi()) -
(*it)->partons().first->position().x()*sin(b.phi()) + b.bVec().y() ) );
if( !((*it)->neighbors().second) )
(*it)->partons().second->position( Parton::Point
((*it)->partons().second->position().x()*cos(b.phi()) +
(*it)->partons().second->position().y()*sin(b.phi()) + b.bVec().x() ,
(*it)->partons().second->position().y()*cos(b.phi()) -
(*it)->partons().second->position().x()*sin(b.phi()) + b.bVec().y() ) );
}
}
}
list< PartonPtr > DipoleState::getPartons() {
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;
}
//This is moved to EventFiller, and should not be called.
void DipoleState::absorbSmallDipoles() {
for(set< DipolePtr >::iterator it = allDipoles.begin(); //which order to check?
it != allDipoles.end(); it++) { //doesn't work for quarks
if ( (*it)->interacted() || (*it)->neighbors().first->interacted() )
continue;
if( !((*it)->children().first || (*it)->children().second) &&
( (*it)->neighbors().first ) &&
( (*it)->neighbors().first->partons().first!=(*it)->partons().second ) ) {
PartonPtr p3 = (*it)->partons().first;
PartonPtr p1 = p3->dipoles().first->partons().first;
PartonPtr p2 = p3->dipoles().second->partons().second;
Parton::Point r1 = p3->position() - p1->position();
Parton::Point r2 = p3->position() - p2->position();
TransverseMomentum e1 = r1/sqr(r1.pt());
TransverseMomentum e2 = r2/sqr(r2.pt());
TransverseMomentum pT2prime = p2->pT() +
e2*(p3->pT().first*e1.second - p3->pT().second*e1.first)/
(e2.first*e1.second-e2.second*e1.first);
TransverseMomentum pT1prime = p1->pT() +
e1*(p3->pT().first*e2.second - p3->pT().second*e2.first)/
(e1.first*e2.second - e1.second*e2.first);
//the following determines the pT spectrum. Think through well!!
if( ( pT1prime.pt() + pT2prime.pt() )/
( p1->pT().pt() + p2->pT().pt() + p3->pT().pt()) < UseRandom::rnd() ) {
reabsorb(p3);
}
}
}
}
vector<DipoleState::String> DipoleState::strings() {
// cout << "entering strings" << endl;
vector<DipoleState::String> ret;
set<DipolePtr>::iterator it = allDipoles.begin();
while( !(allDipoles.empty()) ) {
it = allDipoles.begin();
//if it has children, remove it and continue.
if( (*it)->children().first || (*it)->children().second ) {
allDipoles.erase(it);
continue;
}
//start extracting partons in the forward direction
DipolePtr dip = (*it);
DipolePtr startingDip = dip;
DipoleState::String s;
while( dip->neighbors().second &&
dip->neighbors().second != startingDip ) {
s.push_back( dip->partons().second );
allDipoles.erase( dip );
dip = dip->neighbors().second;
}
//if we looped back, remove last connecting dipole
if( dip->neighbors().second ) {
s.push_back( dip->partons().second );
allDipoles.erase( dip );
}
//if we did not loop back, continue from start in other direction.
if( !(dip->neighbors().second) ) {
cout << dip << "didn't loop back!! O_o"<< endl;
while( dip->neighbors().first ) {
s.insert( s.begin(), dip->partons().first );
allDipoles.erase( dip );
dip = dip->neighbors().first;
}
}
// cout << "storing a string of length " << 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;
Energy plus = Plp + PpShuffled;
Energy minus = Prm + PmShuffled;
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;
diagnosis();
}
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;
}
// cout << "q1 = " << q1 << ", yshift1 = " << abs(log(q1)) << endl
// << "q2 = " << q2 << ", yshift2 = " << abs(log(q2)) << endl;
// cout << "new total p+: " << (q1*Plp + PrpNeeded/q2)/GeV << endl
// << "new total p-: " << (q2*Prm + PlmNeeded/q1)/GeV << 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;
// cout << "before: p+ = " << p->plus()/GeV <<
// ", y = " << p->y() <<
// ", p- = " << p->minus()/GeV << endl;
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() ) );
}
// cout << "after: p+ = " << p->plus()/GeV <<
// ", y = " << p->y() <<
// ", p- = " << p->minus()/GeV << endl;
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::persistentOutput(PersistentOStream & os) const {
os << theHandler << theInitialDipoles << theWFInfo << theWeight << allDipoles;
}
void DipoleState::persistentInput(PersistentIStream & is, int) {
is >> theHandler >> theInitialDipoles >> theWFInfo >> theWeight >> allDipoles;
}
ClassDescription<DipoleState> DipoleState::initDipoleState;
// Definition of the static class description member.
void DipoleState::Init() {}
diff --git a/DIPSY/EventFiller.cc b/DIPSY/EventFiller.cc
--- a/DIPSY/EventFiller.cc
+++ b/DIPSY/EventFiller.cc
@@ -1,317 +1,331 @@
// -*- 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/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/Current.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/EventRecord/Step.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace DIPSY;
EventFiller::EventFiller() {}
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,
const DipoleState & dl, const DipoleState & dr,
const ImpactParameters & b) const {
// Get a list of possible dipole-dipole interactions
FList fl = eh.xSecFn().flist(dl, dr, b);
// Select the interactions which should be performed.
vector<FList::const_iterator> selected = selectInteractions(fl);
vector<String> strings = extractStrings(dl, dr, selected, b);
fillStep(step, inc, strings);
double sum = 0.0;
for ( FList::iterator it = fl.begin(); it != fl.end(); ++it )
sum += it->first.first;
return eh.xSecFn().unitarize(sum);
}
vector<EventFiller::FList::const_iterator>
EventFiller::selectInteractions(const FList & fl) const {
vector<FList::const_iterator> ret;
// First we choose one interaction according to 1-exp(-f). This is
// to avoid the case of no interactions, which is taken care of by
// the event weight.
Selector<FList::const_iterator, double> sel;
for ( FList::const_iterator it = fl.begin(); it != fl.end(); ++it )
sel.insert(it->first.second, it);
FList::const_iterator prim = sel[UseRandom::rnd()];
// Now get all dipole pairs which wants to interact.
vector<FList::const_iterator> potential;
for ( FList::const_iterator it = fl.begin(); it != fl.end(); ++it )
if ( it == prim || it->first.second > UseRandom::rnd() )
potential.push_back(it);
// Mark the dipoles interacted.
for ( int i = 0, N = potential.size(); i < N; ++i )
if ( canInteract(*potential[i]->second.first, *potential[i]->second.second) )
ret.push_back(potential[i]);
return ret;
}
bool EventFiller::canInteract(Dipole & di, Dipole & dj) const {
if ( di.interacted() || dj.interacted() ) return false;
di.interact(dj);
dj.interact(di);
return true;
}
bool rapSort(PartonPtr part1,PartonPtr part2) {
return (abs(part1->y()) < abs(part2->y()));
}
vector<EventFiller::String>
EventFiller::extractStrings(const DipoleState & dl, const DipoleState & dr,
const vector<FList::const_iterator> & sel,
const ImpactParameters & b ) const {
DipoleStatePtr rightState = const_ptr_cast<DipoleStatePtr>(&dr);
DipoleStatePtr leftState = const_ptr_cast<DipoleStatePtr>(&dl);
+
+ double weight = 1.0 - exp(-Current<DipoleEventHandler>()->xSecFn().sumf(*leftState, *rightState, b));
// cout << "entering extractString__________________________________________________________" << endl;
// ofstream outfile ("left.dat");
// leftState->printDipoles(outfile);
// outfile.close();
// ofstream outfile4 ("right.dat");
// rightState->printDipoles(outfile4);
// outfile4.close();
//List and sort partons in rapidity.
list<PartonPtr> virtualPartonsL = leftState->virtualPartons();
virtualPartonsL.sort(rapSort);
list<PartonPtr> virtualPartonsR = rightState->virtualPartons();
virtualPartonsR.sort(rapSort);
//Go through them and reabsorb those that has not interacted.
for ( list<PartonPtr>::const_iterator it = virtualPartonsL.begin() ; it != virtualPartonsL.end(); it++ ) {
leftState->reabsorb((*it));
}
for ( list<PartonPtr>::const_iterator it = virtualPartonsR.begin() ; it != virtualPartonsR.end(); it++ ) {
rightState->reabsorb((*it));
}
list<PartonPtr> partonsL = leftState->getPartons();
list<PartonPtr> partonsR = rightState->getPartons();
partonsL.sort(rapSort); //how should they be sorted? what order to check?
partonsR.sort(rapSort);
bool left = true;
//absorb too small dipoles, depending on some kind of virtuality definition.
for ( list<PartonPtr>::const_iterator it = partonsL.begin() ; it != partonsR.end(); it++ ) {
if( it == partonsL.end() ) {
it = partonsR.begin();
left = false;
}
if ( !((*it)->dipoles().first) || !((*it)->dipoles().second) )
continue;
if ( (*it)->dipoles().first->interacted() || (*it)->dipoles().second->interacted() )
continue;
PartonPtr p1 = (*it)->dipoles().first->partons().first;
PartonPtr p2 = (*it)->dipoles().second->partons().second;
double P1 = sqrt((*it)->dist2(*p2))/( sqrt((*it)->dist2(*p1)) + sqrt((*it)->dist2(*p2)) );
double P2 = sqrt((*it)->dist2(*p1))/( sqrt((*it)->dist2(*p1)) + sqrt((*it)->dist2(*p2)) );
TransverseMomentum pT2prime = p2->pT() + P2*(*it)->pT();
TransverseMomentum pT1prime = p1->pT() + P1*(*it)->pT();
//Here is the physics. what critetion should we use??
if( ( pT1prime.pt() + pT2prime.pt() )/
- ( p1->pT().pt() + p2->pT().pt() + (*it)->pT().pt()) < 0.0*UseRandom::rnd() ) {
+ ( p1->pT().pt() + p2->pT().pt() + (*it)->pT().pt()) < UseRandom::rnd() ) {
if( left )
leftState->reabsorb(*it);
else
rightState->reabsorb(*it);
}
}
//Remove too small dipoles, even if they have intercting children.
//Much physics to be sorted out in here.
// leftState->absorbSmallDipoles();
// rightState->absorbSmallDipoles();
-// //merge the two states to one
-// rightState->mirror();
-// rightState->translate( b );
-// DipoleStatePtr finalState = leftState->merge( rightState );
+ //merge the two states to one
+ rightState->mirror();
+ rightState->translate( b );
+ DipoleStatePtr finalState = leftState->merge( rightState );
-// //do each individual dipole interaction. momentum transfer and colour flow.
-// for( vector< FList::const_iterator>::const_iterator it = sel.begin();
-// it != sel.end(); it++ ) {
+ //do each individual dipole interaction. momentum transfer and colour flow.
+ for( vector< FList::const_iterator>::const_iterator it = sel.begin();
+ it != sel.end(); it++ ) {
-// //tranfser pT in collision as 1/r.
-// PartonPtr p11 = (*it)->second.first->partons().first;
-// PartonPtr p12 = (*it)->second.first->partons().second;
-// PartonPtr p21 = (*it)->second.second->partons().first;
-// PartonPtr p22 = (*it)->second.second->partons().second;
-// TransverseMomentum recoil1 = ( p11->position() - p21->position() )/( p11->dist2( * p21 ) );
-// TransverseMomentum recoil2 = ( p12->position() - p22->position() )/( p12->dist2( * p22 ) );
+ //tranfser pT in collision as 1/r.
+ PartonPtr p11 = (*it)->second.first->partons().first;
+ PartonPtr p12 = (*it)->second.first->partons().second;
+ PartonPtr p21 = (*it)->second.second->partons().first;
+ PartonPtr p22 = (*it)->second.second->partons().second;
+ TransverseMomentum recoil11 = ( p11->position() - p21->position() )/( p11->dist2( * p21 ) );
+ TransverseMomentum recoil12 = ( p11->position() - p22->position() )/( p11->dist2( * p22 ) );
+ TransverseMomentum recoil21 = ( p12->position() - p21->position() )/( p12->dist2( * p21 ) );
+ TransverseMomentum recoil22 = ( p12->position() - p22->position() )/( p12->dist2( * p22 ) );
-// // cout << "recoils: " << recoil1.pt()/GeV<< ", "<<recoil2.pt()/GeV<<endl;
+ if ( p11->plus() < (p11->pT() + recoil11 + recoil12).pt() ||
+ p12->plus() < (p12->pT() + recoil21 + recoil22).pt() ||
+ p21->minus()< (p21->pT() - recoil11 - recoil21).pt() ||
+ p22->minus()< (p22->pT() - recoil21 - recoil22).pt() )
+ continue; //not enough momentum to do the recoil.
-// p11->pT( p11->pT() + recoil1 );
-// p12->pT( p12->pT() + recoil2 );
-// p21->pT( p21->pT() - recoil1 );
-// p22->pT( p22->pT() - recoil2 );
+ if ( recoil11.pt() > 200.0*GeV ) cout << "recoil of " << recoil11.pt()/GeV << endl;
+ if ( recoil12.pt() > 200.0*GeV ) cout << "recoil of " << recoil12.pt()/GeV << endl;
+ if ( recoil21.pt() > 200.0*GeV ) cout << "recoil of " << recoil21.pt()/GeV << endl;
+ if ( recoil22.pt() > 200.0*GeV ) cout << "recoil of " << recoil22.pt()/GeV << endl;
-// // cout<<"old ys: "<<p11->y()<<", "<<p12->y()<<", "<<p21->y()<<", "<<p22->y()<<endl;
-// //disregard p+, p- transfer for now.
-// p11->y( log( p11->pT().pt()/p11->plus() ) );
-// p11->minus( p11->pT().pt()*exp( p11->y() ) );
+ p11->pT( p11->pT() + recoil11 + recoil12 );
+ p12->pT( p12->pT() + recoil21 + recoil22 );
+ p21->pT( p21->pT() - recoil11 - recoil21 );
+ p22->pT( p22->pT() - recoil21 - recoil22 );
-// p12->y( log( p12->pT().pt()/p12->plus() ) );
-// p12->minus( p12->pT().pt()*exp( p12->y() ) );
+// cout<<"old ys: "<<p11->y()<<", "<<p12->y()<<", "<<p21->y()<<", "<<p22->y()<<endl;
+ //disregard p+, p- transfer for now.
+ p11->y( log( p11->pT().pt()/p11->plus() ) );
+ p11->minus( p11->pT().pt()*exp( p11->y() ) );
-// p21->y( log( p21->minus()/p21->pT().pt() ) );
-// p21->plus( p21->pT().pt()*exp( -p21->y() ) );
+ p12->y( log( p12->pT().pt()/p12->plus() ) );
+ p12->minus( p12->pT().pt()*exp( p12->y() ) );
-// p22->y( log( p22->minus()/p22->pT().pt() ) );
-// p22->plus( p22->pT().pt()*exp( -p22->y() ) );
-// // cout<<"new ys: "<<p11->y()<<", "<<p12->y()<<", "<<p21->y()<<", "<<p22->y()<<endl;
+ p21->y( log( p21->minus()/p21->pT().pt() ) );
+ p21->plus( p21->pT().pt()*exp( -p21->y() ) );
-// if ( p11->y()>0.0 || p12->y() > 0.0 || p21->y() < 0.0 || p22->y() < 0.0 )
-// cout<<"partons pushed over to wrong rapidity. probably due to multiple interactions"<<endl;
+ p22->y( log( p22->minus()/p22->pT().pt() ) );
+ p22->plus( p22->pT().pt()*exp( -p22->y() ) );
+// cout<<"new ys: "<<p11->y()<<", "<<p12->y()<<", "<<p21->y()<<", "<<p22->y()<<endl;
-// if( (*it)->first.second < 2.0*UseRandom::rnd() ) { //P_el/P_tot = f_ij/2 correct?
-// (*it)->second.first->swingDipole( (*it)->second.second );
-// (*it)->second.first->recombine();
-// }
-// }
+ if ( p11->y()>0.0 || p12->y() > 0.0 || p21->y() < 0.0 || p22->y() < 0.0 )
+ cout<<"partons pushed over to wrong rapidity. probably due to multiple interactions"<<endl;
+
+ if( (*it)->first.second < 2.0*UseRandom::rnd() ) { //P_el/P_tot = f_ij/2 correct?
+ (*it)->second.first->swingDipole( (*it)->second.second );
+ (*it)->second.first->recombine();
+ }
+ }
//Merge states, and recouple the interacting non-elastic dipole pairs.
- DipoleStatePtr finalState = leftState->collide(rightState,sel,b);
+// DipoleStatePtr finalState = leftState->collide(rightState,sel,b);
// ofstream outfile2 ("mergedStates.dat");
// finalState->printDipoles(outfile2);
// outfile2.close();
list<PartonPtr> partons = finalState->getPartons();
//now shuffle p+ and p- so that everything ends up on shell.
//here it is done assuming no P has been shuffled between states up to now.
Energy p = 0.0*GeV; //availible plus
Energy m = 0.0*GeV; //availible minus
Energy mP = 0.0*GeV; //missing plus
Energy mM = 0.0*GeV; //missing minus
for ( list<PartonPtr>::const_iterator it = partons.begin() ; it != partons.end(); it++ ) {
if( (*it)->rightMoving() ) {
p += (*it)->plus();
mM += (*it)->minus();
}
else {
m += (*it)->minus();
mP += (*it)->plus();
}
}
double q2 = 1.0/2.0 + ( mP / p - mM / m )/2.0 + sqrt( sqr( 1.0 + mP / p - mM / m )/4.0 - mP / p );
double q1 = (q2 - mP/p)/q2;
for ( list<PartonPtr>::const_iterator it = partons.begin() ; it != partons.end(); it++ ) {
if( (*it)->rightMoving() ) {
(*it)->plus( (*it)->plus()*q1 );
(*it)->y( log((*it)->pT().pt()/(*it)->plus() ) );
(*it)->minus( (*it)->pT().pt()*exp( (*it)->y() ) );
}
else {
(*it)->minus( (*it)->minus()*q2 );
(*it)->y( log( (*it)->minus()/(*it)->pT().pt() ) );
(*it)->plus( (*it)->pT().pt()*exp( -(*it)->y() ) );
}
}
histYShift->fill( - log(q1) );
histYShift->fill( - log(q2) );
Energy scalarpT = 0.0*GeV;
for ( list<PartonPtr>::const_iterator it = partons.begin() ; it != partons.end(); it++ ) {
scalarpT += (*it)->pT().pt();
}
- histpT->fill( scalarpT/GeV ); //weight with xsec
+ histpT->fill( scalarpT/GeV, weight ); //weight with xsec
// ofstream outfile3 ("after.dat");
// finalState->printDipoles(outfile3);
// outfile3.close();
//Sort the remaining partons in strings.
vector<String> ret = finalState->strings();
// cout << "b = " << b.bVec().pt()*GeV << endl;
return ret;
}
void 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 ) {
for ( int i = 0, N = strings[is].size(); i < N; ++i ) {
particles[is].push_back(strings[is][i]->produceParticle());
if ( i > 0 ) particles[is][i]->colourConnect(particles[is][i - 1]);
if ( i == N - 1 && strings[is][i]->flavour() == ParticleID::g )
particles[is][0]->colourConnect(particles[is][i]);
sub->addOutgoing(particles[is][i]);
}
}
step.addSubProcess(sub);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void EventFiller::doinitrun() {
generator()->histogramFactory()->initrun();
histYShift = generator()->histogramFactory()->createHistogram1D
("YShift",100,0.0,0.01);
histpT = generator()->histogramFactory()->createHistogram1D
("pT",50,0.0,500.0);
}
void EventFiller::persistentOutput(PersistentOStream & os) const {
// *** ATTENTION *** os << ; // Add all member variable which should be written persistently here.
}
void EventFiller::persistentInput(PersistentIStream & is, int) {
// *** ATTENTION *** is >> ; // Add all member variable which should be read persistently here.
}
ClassDescription<EventFiller> EventFiller::initEventFiller;
// Definition of the static class description member.
void EventFiller::Init() {
static ClassDocumentation<EventFiller> documentation
("The EventFiller class is able to produce an initial ThePEG::Step "
"from two colliding DipoleStates.");
}
diff --git a/DIPSY/TestDIPSY.in b/DIPSY/TestDIPSY.in
--- a/DIPSY/TestDIPSY.in
+++ b/DIPSY/TestDIPSY.in
@@ -1,60 +1,60 @@
mkdir /DIPSY
cd /DIPSY
create DIPSY::DipoleEventHandler LHCEventHandler libDIPSY.so
create DIPSY::ImpactParameterGenerator stdBGenerator
set LHCEventHandler:BGen stdBGenerator
create DIPSY::SimpleProton stdProton
set stdProton:Particle /Defaults/Particles/p+
set LHCEventHandler:WFL stdProton
set LHCEventHandler:WFR stdProton
set LHCEventHandler:PreSamples 1
create DIPSY::TotalXSecAnalysis TotXSec
insert LHCEventHandler:AnalysisHandlers[0] TotXSec
# create DIPSY::ElasticXSecAnalysis ElXSec ElasticXSecAnalysis.so
# set ElXSec:NBins 10000
# set ElXSec:DeltaB 0.005
# insert LHCEventHandler:AnalysisHandlers[0] ElXSec
#cp ElXSec ElXSec1
#set ElXSec1:DeltaB 0.01
#insert LHCEventHandler:AnalysisHandlers[0] ElXSec1
#cp ElXSec ElXSec2
#set ElXSec2:DeltaB 0.02
#insert LHCEventHandler:AnalysisHandlers[0] ElXSec2
#cp ElXSec ElXSec4
#set ElXSec4:DeltaB 0.04
#insert LHCEventHandler:AnalysisHandlers[0] ElXSec4
create DIPSY::PT1DEmitter pt1DEmitter PT1DEmitter.so
create DIPSY::OldStyleEmitter oldEmitter OldStyleEmitter.so
create DIPSY::Emitter stdEmitter
set stdEmitter:RScale 1.0
set LHCEventHandler:Emitter stdEmitter
create DIPSY::Swinger stdSwinger
set LHCEventHandler:Swinger stdSwinger
create DIPSY::EventFiller stdFiller
set LHCEventHandler:EventFiller stdFiller
create DIPSY::DipoleXSec stdXSec
set LHCEventHandler:XSecFn stdXSec
create ThePEG::FixedCMSLuminosity LHCLumi
set LHCLumi:Energy 2000.0
set LHCEventHandler:LuminosityFunction LHCLumi
set LHCEventHandler:LambdaQCD 0.2
set LHCEventHandler:RMax 2.93
set LHCEventHandler:NF 3
create ThePEG::LWHFactory HFac LWHFactory.so
set HFac:StoreType flat
set HFac:Suffix dat
create ThePEG::MultiEventGenerator LHCGenerator MultiEventGenerator.so
set LHCGenerator:RandomNumberGenerator /Defaults/Random
set LHCGenerator:StandardModelParameters /Defaults/StandardModel
set LHCGenerator:EventHandler LHCEventHandler
-set LHCGenerator:NumberOfEvents 1500
+set LHCGenerator:NumberOfEvents 1000
set LHCGenerator:PrintEvent 10
set LHCGenerator:DebugLevel 1
set LHCGenerator:HistogramFactory HFac
# insert LHCGenerator:DefaultObjects[0] stdEmitter
# insert LHCGenerator:DefaultObjects[0] oldStyleEmitter
# insert LHCGenerator:DefaultObjects[0] pt1DEmitter
# do LHCGenerator:AddInterface /DIPSY/LHCEventHandler:Emitter /DIPSY/stdEmitter, /DIPSY/oldStyleEmitter, /DIPSY/pt1DEmitter
# do LHCGenerator:AddInterface /DIPSY/LHCEventHandler:PreSamples 250, 1000, 4000, 16000
saverun TestDIPSY LHCGenerator

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 10:06 AM (1 d, 13 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111047
Default Alt Text
(40 KB)

Event Timeline