Page MenuHomeHEPForge
Files F7876925

No OneTemporary
Actions

View Options

This file is larger than 256 KB, so syntax highlighting was skipped.
This document is not UTF8. It was detected as ISO-8859-1 (Latin 1) and converted to UTF8 for display.
diff --git a/DIPSY/DIPSYDefaults.in b/DIPSY/DIPSYDefaults.in
--- a/DIPSY/DIPSYDefaults.in
+++ b/DIPSY/DIPSYDefaults.in
@@ -1,199 +1,208 @@
rrmdir /DIPSY
mkdir /DIPSY
cd /DIPSY
create ThePEG::Cuts NoCuts
create DIPSY::ImpactParameterGenerator stdBGenerator libDIPSY.so
set stdBGenerator:Width 4.0
create DIPSY::SimpleProton stdProton
set stdProton:Particle /Defaults/Particles/p+
create DIPSY::SimpleProton stdAntiProton
set stdAntiProton:Particle /Defaults/Particles/pbar-
set /Defaults/Particles/p+:PDF /Defaults/Partons/NoPDF
set /Defaults/Particles/pbar-:PDF /Defaults/Partons/NoPDF
create DIPSY::VirtualPhoton virtualPhoton VirtualPhoton.so
set virtualPhoton:Q2 14.0
set virtualPhoton:Polarisation 0
set virtualPhoton:Particle /Defaults/Particles/gamma
create DIPSY::SimpleNucleus Helium SimpleNucleus.so
create DIPSY::SimpleNucleus AntiHelium SimpleNucleus.so
do Helium:SetNucleus He
do AntiHelium:SetNucleus He
set Helium:WF stdProton
set AntiHelium:WF stdProton
set Helium:Particle /Defaults/Particles/p+
set AntiHelium:Particle /Defaults/Particles/pbar-
create DIPSY::SimpleNucleus Lithium SimpleNucleus.so
create DIPSY::SimpleNucleus AntiLithium SimpleNucleus.so
do Lithium:SetNucleus Li
do AntiLithium:SetNucleus Li
set Lithium:WF stdProton
set AntiLithium:WF stdProton
set Lithium:Particle /Defaults/Particles/p+
set AntiLithium:Particle /Defaults/Particles/pbar-
create DIPSY::SimpleNucleus Carbon SimpleNucleus.so
create DIPSY::SimpleNucleus AntiCarbon SimpleNucleus.so
do Carbon:SetNucleus C
do AntiCarbon:SetNucleus C
set Carbon:WF stdProton
set AntiCarbon:WF stdProton
set Carbon:Particle /Defaults/Particles/p+
set AntiCarbon:Particle /Defaults/Particles/pbar-
create DIPSY::SimpleNucleus Copper SimpleNucleus.so
create DIPSY::SimpleNucleus AntiCopper SimpleNucleus.so
do Copper:SetNucleus Cu
do AntiCopper:SetNucleus Cu
set Copper:WF stdProton
set AntiCopper:WF stdProton
set Copper:Particle /Defaults/Particles/p+
set AntiCopper:Particle /Defaults/Particles/pbar-
create DIPSY::SimpleNucleus Oxygen SimpleNucleus.so
create DIPSY::SimpleNucleus AntiOxygen SimpleNucleus.so
do Oxygen:SetNucleus O
do AntiOxygen:SetNucleus O
set Oxygen:WF stdProton
set AntiOxygen:WF stdProton
set Oxygen:Particle /Defaults/Particles/p+
set AntiOxygen:Particle /Defaults/Particles/pbar-
create DIPSY::SimpleNucleus Gold SimpleNucleus.so
create DIPSY::SimpleNucleus AntiGold SimpleNucleus.so
do Gold:SetNucleus Au
do AntiGold:SetNucleus Au
set Gold:WF stdProton
set AntiGold:WF stdProton
set Gold:Particle /Defaults/Particles/p+
set AntiGold:Particle /Defaults/Particles/pbar-
create DIPSY::SimpleNucleus Lead SimpleNucleus.so
create DIPSY::SimpleNucleus AntiLead SimpleNucleus.so
do Lead:SetNucleus Pb
do AntiLead:SetNucleus Pb
set Lead:WF stdProton
set AntiLead:WF stdProton
set Lead:Particle /Defaults/Particles/p+
set AntiLead:Particle /Defaults/Particles/pbar-
create DIPSY::TotalXSecAnalysis TotXSec
create DIPSY::ElasticXSecAnalysis ElXSec ElasticXSecAnalysis.so
set ElXSec:NBins 100
set ElXSec:DeltaB 0.5
set ElXSec:NqBins 200
set ElXSec:Deltaq 0.05
#create DIPSY::PT1DEmitter pt1DEmitter PT1DEmitter.so
#create DIPSY::OldStyleEmitter oldEmitter OldStyleEmitter.so
create DIPSY::Emitter stdEmitter
create DIPSY::Swinger stdSwinger
#create DIPSY::RecoilSwinger recSwinger RecoilSwinger.so
set stdSwinger:Lambda 1.0
create DIPSY::SmallDipoleAbsorber stdAbsorber
create DIPSY::EventFiller stdFiller
set stdFiller:DipoleAbsorber stdAbsorber
create DIPSY::DipoleXSec stdXSec
create ThePEG::FixedCMSLuminosity Lumi
set Lumi:Energy 7000.0
create ThePEG::LWHFactory HFac LWHFactory.so
set HFac:StoreType flat
set HFac:Suffix dat
cp /Defaults/Particles/pi0 pi0
set pi0:Stable Stable
create TheP8I::StringFragmentation Frag8 libTheP8I.so
create DIPSY::FSDipole5Ordering FSOrdering FSDipole5Ordering.so
cp /Ariadne5/Cascade/AriadneCascadePlain AriadneCascade
insert AriadneCascade:Reweighters[0] FSOrdering
create DIPSY::FSAnalysis FSAnalysis FSAnalysis.so
create DIPSY::DipoleEventHandler EventHandler libDIPSY.so
set EventHandler:WFL stdProton
set EventHandler:WFR stdProton
set EventHandler:PartonExtractor /Defaults/Handlers/StandardExtractor
set EventHandler:Cuts NoCuts
set EventHandler:BGen stdBGenerator
set EventHandler:PreSamples 100
# set EventHandler:CollisionType NonDiffractive
insert EventHandler:AnalysisHandlers[0] TotXSec
insert EventHandler:AnalysisHandlers[0] ElXSec
set EventHandler:EventFiller stdFiller
set EventHandler:XSecFn stdXSec
set EventHandler:Emitter stdEmitter
set EventHandler:Swinger stdSwinger
create ThePEG::MultiEventGenerator Generator MultiEventGenerator.so
set Generator:RandomNumberGenerator /Defaults/Random
set Generator:StandardModelParameters /Defaults/StandardModel
set Generator:EventHandler EventHandler
set Generator:NumberOfEvents 10000
set Generator:PrintEvent 0
set Generator:DebugLevel 0
set Generator:HistogramFactory HFac
set Generator:DumpPeriod 1000
insert Generator:LocalParticles[0] pi0
set EventHandler:DecayHandler /Defaults/Handlers/StandardDecayHandler
set EventHandler:HadronizationHandler Frag8
set EventHandler:CascadeHandler AriadneCascade
insert Generator:AnalysisHandlers[0] FSAnalysis
cp EventHandler LHCEventHandler
set LHCEventHandler:WFL stdProton
set LHCEventHandler:WFR stdProton
cp Lumi LHCLumi
set LHCEventHandler:LuminosityFunction LHCLumi
cp Generator LHCGenerator
set LHCGenerator:EventHandler LHCEventHandler
set LHCLumi:Energy 7000.0
cp EventHandler TevatronEventHandler
set TevatronEventHandler:WFL stdAntiProton
set TevatronEventHandler:WFR stdProton
cp Lumi TevatronLumi
set TevatronEventHandler:LuminosityFunction TevatronLumi
cp Generator TevatronGenerator
set TevatronGenerator:EventHandler TevatronEventHandler
set TevatronLumi:Energy 1800.0
cp EventHandler HeraEventHandler
set HeraEventHandler:WFL virtualPhoton
set HeraEventHandler:WFR stdProton
cp Lumi HeraLumi
set HeraEventHandler:LuminosityFunction HeraLumi
cp Generator HeraGenerator
set HeraGenerator:EventHandler HeraEventHandler
set HeraLumi:Energy 220.0
set virtualPhoton:Q2 14.0
cp EventHandler RHICEventHandler
set RHICEventHandler:WFL AntiGold
set RHICEventHandler:WFR Gold
cp Lumi RHICLumi
set RHICEventHandler:LuminosityFunction RHICLumi
cp Generator RHICGenerator
set RHICGenerator:EventHandler RHICEventHandler
set RHICLumi:Energy 39400.0
cp EventHandler LHCPbPbEventHandler
set LHCPbPbEventHandler:WFL AntiLead
set LHCPbPbEventHandler:WFR Lead
cp Lumi LHCPbPbLumi
set LHCPbPbEventHandler:LuminosityFunction LHCPbPbLumi
cp Generator LHCPbPbGenerator
set LHCPbPbGenerator:EventHandler LHCPbPbEventHandler
set LHCPbPbLumi:Energy 574000.0
cp EventHandler eRHICEventHandler
set eRHICEventHandler:WFL virtualPhoton
set eRHICEventHandler:WFR Gold
cp Lumi eRHICLumi
set eRHICEventHandler:LuminosityFunction eRHICLumi
cp Generator eRHICGenerator
set eRHICGenerator:EventHandler eRHICEventHandler
set eRHICLumi:Energy 1404.0
+cp EventHandler LHCpPbEventHandler
+set LHCpPbEventHandler:WFL stdProton
+set LHCpPbEventHandler:WFR Lead
+cp Lumi LHCpPbLumi
+set LHCpPbEventHandler:LuminosityFunction LHCpPbLumi
+cp Generator LHCpPbGenerator
+set LHCpPbGenerator:EventHandler LHCpPbEventHandler
+set LHCpPbLumi:Energy 40000.0
+
cd /
diff --git a/DIPSY/DipoleState.cc b/DIPSY/DipoleState.cc
--- a/DIPSY/DipoleState.cc
+++ b/DIPSY/DipoleState.cc
@@ -1,1941 +1,1900 @@
// -*- 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/Utilities/Throw.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/MaxCmp.h"
#include "ThePEG/Utilities/ObjectIndexer.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), theWFInfo(x.theWFInfo),
theWeight(x.theWeight), doTakeHistory(x.doTakeHistory), theYmax(x.theYmax),
theCollidingEnergy(x.theCollidingEnergy), theHistory(x.theHistory),
allDipoles(x.allDipoles) {
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<Ariadne::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());
}
return copy;
}
void DipoleState::printDipoles(ostream & output) const {
for(set< DipolePtr >::iterator it = allDipoles.begin();
it != allDipoles.end(); it++) {
if( ((*it)->children().first || (*it)->children().second) )
continue;
int i = (*it)->colour();
output << i << '\t' << (*it)->partons().first->position().first*GeV << '\t' <<
(*it)->partons().first->position().second*GeV << '\t' <<0<< '\t' <<
(*it)->partons().first->rightMoving() << '\t' <<
(*it)->partons().first->interacted() << '\t' <<
bool((*it)->interacted()) << '\t' <<
(*it)->DGLAPsafe() << '\t' <<
(*it)->partons().first->y() << '\t' <<
(*it)->partons().first->pT().pt()/GeV << endl;
output << i << '\t' <<(*it)->partons().second->position().first*GeV << '\t' <<
(*it)->partons().second->position().second*GeV << '\t' <<0<< '\t' <<
(*it)->partons().second->rightMoving() << '\t' <<
(*it)->partons().second->interacted() << '\t' <<
bool((*it)->interacted()) << '\t' <<
(*it)->DGLAPsafe() << '\t' <<
(*it)->partons().second->y() << '\t' <<
(*it)->partons().second->pT().pt()/GeV <<endl<<endl;
output << i << '\t' << (*it)->partons().first->position().first*GeV << '\t' <<
(*it)->partons().first->position().second*GeV << '\t' << 1 <<
'\t' << (*it)->partons().first->y()<< endl;
output << i << '\t' << (*it)->partons().first->position().first*GeV << '\t' <<
(*it)->partons().first->position().second*GeV << '\t' << 2 << '\t'<<
(*it)->partons().first->pT().first/GeV << '\t' <<
(*it)->partons().first->pT().second/GeV << '\t'<<
(*it)->partons().first->y()<< '\t' <<
(*it)->partons().first->rightMoving() << '\t' <<
(*it)->partons().first->plus()/GeV << '\t' <<
(*it)->partons().first->minus()/GeV << '\t' <<
bool((*it)->partons().first->valence()) << '\t' <<
(*it)->partons().first->pT().pt()/GeV << endl << endl;
output << i << '\t' << (*it)->partons().second->position().first*GeV <<'\t'<<
(*it)->partons().second->position().second*GeV << '\t' << 1 <<
'\t' << (*it)->partons().second->y()<< endl;
output << i << '\t' << (*it)->partons().second->position().first*GeV <<'\t'<<
(*it)->partons().second->position().second*GeV << '\t' << 2 << '\t'<<
(*it)->partons().second->pT().first/GeV << '\t' <<
(*it)->partons().second->pT().second/GeV << '\t'<<
(*it)->partons().second->y()<< '\t' <<
(*it)->partons().second->rightMoving() << '\t' <<
(*it)->partons().second->plus()/GeV << '\t' <<
(*it)->partons().second->minus()/GeV << '\t' <<
bool((*it)->partons().second->valence()) << '\t' <<
(*it)->partons().second->pT().pt()/GeV << endl<< endl;
}
for ( int i = 0; i < int(theInitialDipoles.size()); i++ ) {
if ( !(theInitialDipoles[i]->participating()) ) {
output << i << '\t' << theInitialDipoles[i]->partons().first->position().first*GeV << '\t' <<
theInitialDipoles[i]->partons().first->position().second*GeV << '\t' <<3<< '\t' <<
theInitialDipoles[i]->partons().first->rightMoving() << '\t' <<
theInitialDipoles[i]->partons().first->interacted() << '\t' <<
bool(theInitialDipoles[i]->interacted()) << '\t' <<
theInitialDipoles[i]->DGLAPsafe() << '\t' <<
theInitialDipoles[i]->partons().first->y() << endl;
output << i << '\t' <<theInitialDipoles[i]->partons().second->position().first*GeV << '\t' <<
theInitialDipoles[i]->partons().second->position().second*GeV << '\t' <<3<< '\t' <<
theInitialDipoles[i]->partons().second->rightMoving() << '\t' <<
theInitialDipoles[i]->partons().second->interacted() << '\t' <<
bool(theInitialDipoles[i]->interacted()) << '\t' <<
theInitialDipoles[i]->DGLAPsafe() << '\t' << theInitialDipoles[i]->partons().second->y() <<endl<<endl;
}
}
for ( int i = 0; i < int(theInitialDipoles.size()); i++ ) {
if ( (theInitialDipoles[i]->participating()) ) {
output << i << '\t' << theInitialDipoles[i]->partons().first->position().first*GeV << '\t' <<
theInitialDipoles[i]->partons().first->position().second*GeV << '\t' <<4<< '\t' <<
theInitialDipoles[i]->partons().first->rightMoving() << '\t' <<
theInitialDipoles[i]->partons().first->interacted() << '\t' <<
bool(theInitialDipoles[i]->interacted()) << '\t' <<
theInitialDipoles[i]->DGLAPsafe() << '\t' <<
theInitialDipoles[i]->partons().first->y() << '\t' <<
theInitialDipoles[i]->partons().first->pT().pt()/GeV << endl;
output << i << '\t' <<theInitialDipoles[i]->partons().second->position().first*GeV << '\t' <<
theInitialDipoles[i]->partons().second->position().second*GeV << '\t' <<4<< '\t' <<
theInitialDipoles[i]->partons().second->rightMoving() << '\t' <<
theInitialDipoles[i]->partons().second->interacted() << '\t' <<
(theInitialDipoles[i]->interacted()) << '\t' <<
theInitialDipoles[i]->DGLAPsafe() << '\t' <<
theInitialDipoles[i]->partons().second->y() << '\t' <<
theInitialDipoles[i]->partons().second->pT().pt()/GeV <<endl<<endl;
}
}
}
void DipoleState::sortDipolesFS() {
theSwingCandidates.clear();
theSwingCandidates.resize(Current<DipoleEventHandler>()->nColours());
list<DipolePtr> dips = getDipoles();
for ( list<DipolePtr>::const_iterator it = dips.begin(); it != dips.end(); it++ ) {
sortDipoleFS(**it);
}
}
void DipoleState::makeBinaryCollision(DipoleStatePtr dr, ImpactParameters b) {
//set range of binary interaction
InvEnergy range = 2.0/GeV;
list<PartonPtr> gluons = getPartons();
list<PartonPtr> otherGluons = dr->getPartons();
//loop through all pairs of gluons
for (list<PartonPtr>::const_iterator it = gluons.begin(); it != gluons.end(); it++) {
for (list<PartonPtr>::const_iterator jt = otherGluons.begin(); jt != otherGluons.end(); jt++) {
//if the gluons are within range, place a single gluon in mid rapidity between them.
PartonPtr p = *it;
PartonPtr op = *jt;
if ( (p->position() - (op->position() + b.bVec())).pt() < range ) {
DipolePtr dip = p->dipoles().first;
PartonPtr g = new_ptr(Parton());
g->position((p->position() + (op->position() + b.bVec()))/2.0);
dip->effectivePartons(EffectiveParton::create(*(dip->partons().first),
dip->size()/2.0),
EffectiveParton::create(*(dip->partons().second),
dip->size()/2.0));
dip->generatedGluon(g);
dip->generatedY(0.0);
dip->emit();
}
}
}
}
void DipoleState::makeGlauberCollision(DipoleStatePtr dr, ImpactParameters b) {
//set range of interaction
InvEnergy range = 5.0/GeV;
list<PartonPtr> gluons = getPartons();
list<PartonPtr> otherGluons = dr->getPartons();
vector<pair<Parton::Point, DipolePtr> > ln, rn;
for (list<PartonPtr>::const_iterator it = gluons.begin(); it != gluons.end(); it++,it++,it++) {
//identify nucleons and fill ln
PartonPtr p1 = *it;
PartonPtr p2 = p1->dipoles().first->partons().first;
PartonPtr p3 = p2->dipoles().first->partons().first;
ln.push_back(make_pair((p1->position() + p2->position() + p3->position())/3.0,
p1->dipoles().first));
}
for (list<PartonPtr>::const_iterator it = otherGluons.begin(); it != otherGluons.end(); it++,it++,it++) {
//identify nucleons and fill rn
PartonPtr p1 = *it;
PartonPtr p2 = p1->dipoles().first->partons().first;
PartonPtr p3 = p2->dipoles().first->partons().first;
rn.push_back(make_pair((p1->position() + p2->position() + p3->position())/3.0 + b.bVec(),
p1->dipoles().first));
}
for (int i = 0; i < int(ln.size());i++) {
bool participant = false;
for (int j = 0; j < int(rn.size());j++) {
//if there is a rn in range of this ln, set participant true and leave this for loop
if ( (ln[i].first - rn[j].first).pt() < range ) {
participant = true;
break;
}
}
if ( participant ) {
//create gluon at y = 0;
DipolePtr dip = ln[i].second;
PartonPtr g = new_ptr(Parton());
g->position(ln[i].first);
dip->effectivePartons(EffectiveParton::create(*(dip->partons().first),
dip->size()/2.0),
EffectiveParton::create(*(dip->partons().second),
dip->size()/2.0));
dip->generatedGluon(g);
dip->generatedY(0.0);
dip->emit();
}
}
for (int i = 0; i < int(rn.size());i++) {
bool participant = false;
DipolePtr intDip;
for (int j = 0; j < int(ln.size());j++) {
//if there is a rn in range of this ln, set participant true and leave this for loop
if ( (rn[i].first - ln[j].first).pt() < range ) {
participant = true;
intDip = ln[j].second;
break;
}
}
if ( participant ) {
//create gluon at y = 0;
DipolePtr dip = intDip;
PartonPtr g = new_ptr(Parton());
g->position(rn[i].first);
dip->effectivePartons(EffectiveParton::create(*(dip->partons().first),
dip->size()/2.0),
EffectiveParton::create(*(dip->partons().second),
dip->size()/2.0));
dip->generatedGluon(g);
dip->generatedY(0.0);
dip->emit();
}
}
}
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();
}
}
}
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);
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 theSwingCandidates[d.colour()].push_back(& d);
}
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 if ( !(d.participating()) );
else if ( !(d.isOn) );
else
theSwingCandidates[d.colour()].push_back(& d);
}
void DipoleState::evolve(double ymin, double ymax) {
save();
int count = 0;
double ycount = ymin;
// cout << "starting at y = " << ycount << endl;
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();
if ( ymin - ycount > 0.1 ) {
// cout << "now at y = " << ymin << endl;
count++;
// printForMovie(ymin, count);
ycount = ymin;
}
sel->emit();
save();
continue;
}
break;
}
theYmax = ymax;
}
void DipoleState::printForMovie(double y, int count) {
ostringstream os;
os << "movieState" << this << count << ".dat";
ofstream outfile (os.str().c_str());
for(set< DipolePtr >::iterator it = allDipoles.begin();
it != allDipoles.end(); it++) {
if( ((*it)->children().first || (*it)->children().second) )
continue;
tDipolePtr dip = *it;
tPartonPtr p1 = dip->partons().first;
tPartonPtr p2 = dip->partons().second;
outfile << p1->position().x()*GeV << '\t' << p1->position().y()*GeV << '\t'
<< y << '\t' << (p1->plus()+p1->minus())/GeV/2.0 << '\t'
<< p2->position().x()*GeV << '\t' << p2->position().y()*GeV << '\t'
<< y << '\t' << (p2->plus()+p2->minus())/GeV/2.0 << endl;
}
outfile.close();
}
void DipoleState::swingFS(double ymin, double ymax) {
// static CPUClock cpuclock("DIPSY::DipoleState::swingFS(loop)");
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 ) {
// CPUTimer timer(cpuclock);
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());
}
theYmax = ymax;
}
-void DipoleState::FSSwing(double ymin, double ymax) {
- save();
-// double ycount = ymin;
-// cout << "starting at y = " << ycount << endl;
- for ( set<DipolePtr>::const_iterator it = allDipoles.begin(); it != allDipoles.end(); it++ ) {
- (*it)->reset();
- (*it)->touch();
- }
- while ( true ) {
- sortFSDipoles();
- tDipolePtr sel = tDipolePtr();
- for ( int i = 0, N = initialDipoles().size(); i < N; ++i ) {
- tDipolePtr si = initialDipoles()[i]->getFSSwinger(ymin, ymax);
- if ( si && si->generatedY() < ymax && si->hasGen() &&
- ( !sel || si->generatedY() < sel->generatedY() )) {
- sel = si;
- }
- }
- if ( sel ) {
- ymin = sel->generatedY();
-// if ( ymin - ycount > 0.1 ||
-// (ymin > -1 && ymin - ycount > 0.04) ||
-// (ymin > -0.2 && ymin - ycount > 0.01) ) {
-// cout << "now at y = " << ymin << endl;
-// double dummy = UseRandom::rnd();
-// ycount = ymin;
-// }
- sel->emit();
- save();
- continue;
- }
- break;
- }
- 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 ) {
if ( Debug::level > 5 ) cout << "this doesn't look like protons. Will not do valence charge normalisation." << endl;
return;
}
if ( mode == 1 ) {
if ( Debug::level > 5 ) {
cout << "now do colour normalisation" << endl;
}
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();
}
}
if ( Debug::level > 5 ) {
cout << "done swinging" << endl;
}
}
}
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()) )
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;
if ( Debug::level > 5 ) cout << "start loop at head " << head->oY() << endl;
//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);
if ( Debug::level > 5 ) cout << " adding " << p->y() << " at front " << endl;
}
//remove from partons, so we know its added.
//could potentially be better to do this without destroy the dipoleState...
partons.erase(p);
//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;
if ( Debug::level > 5 ) cout << " found forward quark, start over in other direction from "
<< p->y() << endl;
}
else {
//want to turn around, but "head" is at the other end, so we are done.
p = head;
if ( Debug::level > 5 ) cout << " found forward quark, but no more partons on other side."
<< endl;
}
forward = false;
}
}
else {
if ( p->dipoles().first ) {
p = p->dipoles().first->partons().first;
if ( Debug::level > 5 ) cout << " continue backwards to " << p->y() << endl;
}
else {
if ( Debug::level > 5 ) cout << " found antiquark, done" << endl;
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;
}
void DipoleState::showHistory() {
if ( !doTakeHistory )
return;
if ( theHistory.size() == 0 ) {
cout << "empty history." << endl;
return;
}
cout << "*******************New History*******************\nUse\n n + enter (next)\n p + enter (previous)\n l + enter (last)\n f + enter (first)\n q + enter (next event)\n e + enter (evolved state)\nto scroll between intermediate states.\nTo plot the state in transverse space, run load plotstate.gp in gnuplot.\nTo plot the state in rapidity (x-axis) and transverse y, run load plotlogstate.gp in gnuplot.\n" << endl;
int i = 0;
// FILE *pipe = popen("gnuplot -persist","w");
char c = 0;
while (true) {
if( int(c) != 10 )
cout << "state " << i+1 << " out of " << theHistory.size() << endl;
ofstream outfile ("state.dat");
theHistory[i]->printDipoles(outfile);
outfile.close();
// fprintf(pipe, "load 'plotstate.gp'\n");
// system("gnuplot 'plotstate.gp'");
c = getchar();
if ( c == 'n' ) {
i++;
if ( i == int(theHistory.size()) ) {
cout << "that was last state already." << endl;
i--;
}
}
else if (c == 'p') {
i--;
if ( i == -1 ) {
cout << "that was first state already." << endl;
i++;
}
}
else if (c == 'l') {
i = theHistory.size()-1;
}
else if (c == 'f') {
i = 0;
}
else if (c == 'e') {
i = 0;
double startY = theHistory[i]->ymax();
while ( theHistory[i]->ymax() == startY )
i++;
}
else if ( c == 'q' ) {
cout << "shutting down this history." << endl;
break;
}
}
// fprintf(pipe, "exit\n");
// fclose(pipe);
}
double DipoleState::avYInInt() const {
double number = 0.0;
double length = 0.0;
list<DipolePtr> dips = getDipoles();
for(list< DipolePtr >::iterator it = dips.begin(); it != dips.end(); it++) {
if ( (*it)->interacted() ) {
number++;
if ( abs((*it)->partons().first->y()) > 100.0 )
cout << "parton at oY " << (*it)->partons().first->oY() << " has y "
<< (*it)->partons().first->y() << endl;
if ( abs((*it)->partons().second->y()) > 100.0 )
cout << "parton at oY " << (*it)->partons().second->oY() << " has y "
<< (*it)->partons().second->y() << endl;
length += abs((*it)->partons().first->y()-(*it)->partons().second->y());
}
}
if (number == 0.0 ) return 0.0;
return length/number;
}
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 ( !(abs((plus + plusNeeded)/GeV - 2000.0) < 0.01 || abs(plus/GeV) < 1.0) ||
// !(abs((minus + minusNeeded)/GeV - 2000.0) < 0.01 || abs(minus/GeV) < 1.0) ) {
// cout << "total p+, p- not consistent with CoME = 2000GeV!!!!" << endl;
// ok = false;
// }
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;
// 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->rightMoving() )
// cout << "| -- " << p->oY() << " has minus " << p->minus()/GeV << ", " << p->plus()/GeV << ", (" << p->pT().x()/GeV << ", " << p->pT().y()/GeV << ")" << endl;
// else
// cout << "| ++ " << p->oY() << " has plus " << p->plus()/GeV << ", " << p->minus()/GeV << ", (" << p->pT().x()/GeV << ", " << p->pT().y()/GeV << ")" << endl;
// }
cout << "----------------------------------------------------------------------\n";
if (ok ) cout << "| Found no errors! :)\n";
else cout << "| STATE NOT OK!! :o\n";
cout << "----------------------------------------------------------------------\n";
plotState(print);
}
if ( Debug::level > 1 ) cout << "dipole state checked, ok: " << ok << endl;
return ok;
}
void DipoleState::plotState(bool pause) const {
if ( Debug::level <= 5 ) return;
ofstream outfile ("state.dat");
printDipoles(outfile);
outfile.close();
cout << "State printed to state.dat." << endl;
if ( pause ) {
cout << "press any key to continue." << endl;
getchar();
}
}
void DipoleState::recursiveReabsorb( PartonPtr p ) {
DipolePtr d1 = p->dipoles().first;
DipolePtr d2 = p->dipoles().second;
reabsorb(p);
if ( d1->size() < d2->size() && d1->neighbors().first->size() < d1->size() )
recursiveReabsorb( d1->partons().first );
if ( d2->size() < d1->size() && d2->neighbors().second->size() < d2->size() )
recursiveReabsorb( d2->partons().second );
}
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, 0.0, 0.0) ) 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, exactly reverses emissions.
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, 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();
for ( int c=0;c < Current<DipoleEventHandler>()->nColours();c++ ) {
d1->colour(c);
d1->forceGenerateRec(ymax2);
}
d1->colour(trueColour);
foundFirstDipole = d1->swingDipole();
trueColour = d2->colour();
for ( int c=0;c < Current<DipoleEventHandler>()->nColours();c++ ) {
d2->colour(c);
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());
allDipoles.insert(otherState->allDipoles.begin(),otherState->allDipoles.end());
//merge the histories
if ( doTakeHistory )
for ( int i = 0; i < int(otherState->theHistory.size()); i++ )
theHistory.push_back(otherState->theHistory[i]);
//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() ) { //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(double y0) {
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( 2.0*y0 - p->y() );
Energy plus = p->plus();
p->plus( p->minus() );
p->minus( plus );
p->rightMoving( !(p->rightMoving()) );
if( !((*it)->neighbors().second) ) {
PartonPtr p2 = (*it)->partons().second;
p2->y( 2.0*y0 - p2->y() );
p2->plus( p2->pT().pt()*exp(-p2->y()) );
p2->minus( p2->pT().pt()*exp(p2->y()) );
p->rightMoving( !(p->rightMoving()) );
}
}
}
Energy x = thePlus;
thePlus = theMinus;
theMinus = x;
}
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() ) );
(*it)->partons().first->pT( b.rotatePT((*it)->partons().first->pT()));
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() ) );
(*it)->partons().second->pT( b.rotatePT((*it)->partons().second->pT()));
}
}
}
}
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;
}
//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);
}
}
}
}
void DipoleState::makeOriginal() {
// cout << "before make original: " << endl;
// coutData();
//clear the original dipoles
theInitialDipoles.clear(); //owns the dipoles
theSwingCandidates.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++;
}
// cout << "after makeoriginal: " << endl;
// coutData();
// for( list<PartonPtr>::iterator it = partons.begin();
// it != partons.end(); it++ )
// (**it).coutData();
// list<DipolePtr> dipoles = getDipoles();
// for( list<DipolePtr>::iterator it = dipoles.begin();
// it != dipoles.end(); it++ )
// (**it).coutData();
}
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;
// *** ATTENTION *** not used: Energy plus = Plp + PpShuffled;
// *** ATTENTION *** not used: 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;
}
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() ) );
}
}
}
}
bool DipoleState::rapiditySorted(const list<PartonPtr> & loop1, const list<PartonPtr> & loop2) {
return ( rapidityRange(loop1).first <= rapidityRange(loop2).first );
}
int DipoleState::nRapidityGaps() const {
set< list<PartonPtr> > temp = loops();
multiset< pair<double,double> > ranges;
for ( set< list<PartonPtr> >::iterator it = temp.begin(); it != temp.end(); ++it )
ranges.insert(rapidityRange(*it));
int nGaps = 0;
double highest = ranges.begin()->second;
for ( multiset< pair<double,double> >::iterator it = ranges.begin();
it != ranges.end(); it++ ) {
if ( it->first > highest ) nGaps++;
if ( it->second > highest ) highest = it->second;
}
return nGaps;
}
pair<double, double> DipoleState::rapidityRange(const list<PartonPtr> & loop ) {
double lower = (*(loop.begin()))->y();
double higher = lower;
for ( list<PartonPtr>::const_iterator it = loop.begin(); it != loop.end(); it++ ) {
PartonPtr p = *it;
if ( p->y() < lower ) lower = p->y();
else if ( p->y() > higher ) higher = p->y();
}
return make_pair(lower, higher);
}
void DipoleState::turnIntoDipole() {
allDipoles.clear();
theInitialDipoles.clear();
thePlus = 200*GeV;
Energy2 Q2 = 100.0*GeV2;
theMinus = Q2/(thePlus);
PartonPtr q1 = new_ptr(Parton());
PartonPtr q2 = new_ptr(Parton());
DipolePtr d = createDipole();
d->partons(make_pair(q1, q2));
generateColourIndex(d);
q1->plus(thePlus*0.55);
q2->plus(thePlus*0.45);
double phi = 2.0*M_PI*UseRandom::rnd();
q1->position(Parton::Point(cos(phi)*1.0/sqrt(Q2), sin(phi)*1.0/sqrt(Q2)));
q2->position(Parton::Point(-cos(phi)*1.0/sqrt(Q2), -sin(phi)*1.0/sqrt(Q2)));
q1->dipoles(make_pair(DipolePtr(), d));
q2->dipoles(make_pair(d, DipolePtr()));
q1->pT(q1->recoil(q2));
q2->pT(q2->recoil(q1));
q1->updateYMinus();
q2->updateYMinus();
q1->oY(q1->y());
q2->oY(q2->y());
q1->valence(true);
q2->valence(true);
q1->valencePT(q1->pT());
q2->valencePT(q2->pT());
q1->valencePlus(q1->plus());
q2->valencePlus(q2->plus());
theInitialDipoles.push_back(d);
}
void DipoleState::saveGluonsToFile(double weight) {
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;
- // os << "/scratch/parton/christof/DIPSY/events/" << filename << "/event" << event->number() << ".dat";
+ //TODO: interface!
+ os << "/home/christoffer/DIPSYevents/" << filename << "/event" << event->number() << ".dat";
- os << "/scratch/parton/christof/DIPSY/events/" << filename << "/event" << weight << ".dat";
+ // os << "/scratch/parton/christof/DIPSY/events/" << filename << "/event" << weight << ".dat";
weight = 2.0*int(sqrt(sqr(bx)+sqr(by))) + 1.0;
- // if ( weight == 1.0 )
- // os << "/scratch/parton/christof/DIPSY/events/" << filename << "/event" << 400 + event->number() << ".dat";
- // if ( weight == 2.0 )
- // os << "/scratch/parton/christof/DIPSY/events/" << filename << "/event" << 900 + event->number() << ".dat";
- // weight = 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() {
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::coutData() {
cout << "printing data for dipolestate " << this << endl;
cout << "thePlus: " << thePlus/GeV
<< ", theMinus: " << theMinus/GeV
<< ", theMinusDeficit: " << theMinusDeficit/GeV
<< ", theHandler: " << theHandler
<< ", theInitialDipoles: ";
for( int i = 0; i < int(theInitialDipoles.size()); i++ )
cout << theInitialDipoles[i] << ", ";
cout << "theSwingCandidates: ";
for( int i = 0; i < int(theSwingCandidates.size()); i++ ) {
cout << "colour " << i << ": ";
for( int j = 0; j < int(theSwingCandidates[i].size()); j++ ) {
cout << theSwingCandidates[i][j] << ", ";
}
}
cout << "theWFInfo: " << theWFInfo
<< ", theWeight: " << theWeight
<< ", doTakeHistory: " << doTakeHistory
<< ", theYmax: " << theYmax
<< ", theCollidingEnergy: " << theCollidingEnergy/GeV
<< ", theHistory: ";
for( int i = 0; i < int(theHistory.size()); i++ )
cout << theHistory[i] << ", ";
cout << "allDipoles: ";
for( set<DipolePtr>::iterator it = allDipoles.begin();
it != allDipoles.end(); it++ )
cout << *it << ", ";
cout << endl;
}
void DipoleState::persistentOutput(PersistentOStream & os) const {
os << ounit(thePlus, GeV) << ounit(theMinus, GeV)
<< ounit(theMinusDeficit, GeV) << theHandler << theInitialDipoles
<< theSwingCandidates << theWFInfo << theWeight << doTakeHistory
<< theYmax << ounit(theCollidingEnergy, GeV) << allDipoles;
}
void DipoleState::persistentInput(PersistentIStream & is, int) {
is >> iunit(thePlus, GeV)>> iunit(theMinus, GeV)
>> iunit(theMinusDeficit, GeV) >> theHandler >> theInitialDipoles
>> theSwingCandidates >> theWFInfo >> theWeight >> doTakeHistory
>> theYmax >> iunit(theCollidingEnergy, GeV) >> allDipoles;
}
// Static variable needed for the type description system in ThePEG.
#include "ThePEG/Utilities/DescribeClass.h"
DescribeClass<DipoleState,PersistentBase>
describeDIPSYDipoleState("DIPSY::DipoleState", "libDIPSY.so");
void DipoleState::Init() {}
diff --git a/DIPSY/DipoleState.h b/DIPSY/DipoleState.h
--- a/DIPSY/DipoleState.h
+++ b/DIPSY/DipoleState.h
@@ -1,572 +1,567 @@
// -*- C++ -*-
#ifndef DIPSY_DipoleState_H
#define DIPSY_DipoleState_H
//
// This is the declaration of the DipoleState class.
//
#include "ThePEG/Config/ThePEG.h"
#include "DipoleState.fh"
#include "Dipole.h"
#include "DipoleEventHandler.fh"
#include "WFInfo.h"
#include "DipoleXSec.h"
namespace DIPSY {
using namespace ThePEG;
/**
* Here is the documentation of the DipoleState class.
*/
class DipoleState: public PersistentBase {
public:
/**
* Copy FList typedef from DipoleXSec.
*/
typedef DipoleXSec::FList FList;
/**
* A String is simply a vector of colour-connected partons.
*/
typedef vector<PartonPtr> String;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The standard constructor taking the controlling
* DipoleEventHandler as argument.
*/
inline DipoleState(const DipoleEventHandler &, WFInfoPtr wf = WFInfoPtr());
/**
* The default constructor.
*/
inline DipoleState();
/**
* The copy constructor.
*/
DipoleState(const DipoleState &);
/**
* The destructor.
*/
virtual ~DipoleState();
//@}
public:
/**
* Runs through a final state evolution checking for swings only.
*/
- void FSSwing(double ymin, double ymax);
-
- /**
- * Runs through a final state evolution checking for swings only.
- */
void swingFS(double ymin, double ymax);
/**
* Does some colour recconection for the valence partons to account for
* original 3 colour colour structure of a proton.
*/
void normaliseValenceCharge(int mode);
/**
* Makes sure the non-participating nucleons of an AA collision have the
* colour flow of the original triangles.
*/
void restoreNonparticipants();
/**
* counts and returns the number of spectating nucleons.
*/
int numberOfSpectators();
/**
* Makes sure the the two partons are conected by a dipole with p1 as first parton.
* returns false if it fails, and they are left not connected.
*/
bool restoreDipole(PartonPtr p1, PartonPtr p2);
/**
* Evolve this state from the given minimum to the given maximum
* rapidity, assuming there will be p- coming from the other state.
*/
void evolve(double ymin, double ymax);
/**
* Makes the state merge with another colliding state.
*/
DipoleStatePtr collide(DipoleStatePtr otherState,
const vector<FList::const_iterator> & sel,
const ImpactParameters & b);
/**
* Makes the state merge with another colliding state.
*/
DipoleStatePtr merge(DipoleStatePtr otherState);
/**
* mirrors the state in rapidity around y0.
*/
void mirror(double y0);
/**
* moves all partons and their pT according to the imparctparameter b.
*/
void translate(const ImpactParameters & b );
/**
* returns the positions of all partons in the state,
* together with the size of the largest neighbouring dipole.
*/
vector<pair<Parton::Point, InvEnergy> > points();
/**
* Checks through all partons and reabsorbs the ones which are
* considered enough gain in some kind of virtuality.
*/
void absorbSmallDipoles();
/**
* Removes all the parents of the dipoles, and sets the currently
* active dipoles as the initial ones.
**/
void makeOriginal();
/**
* returns teh average distance in rapidity of the interacting dipoles.
*/
double avYInInt() const;
/**
* Controls and couts various info about the state. Returns false
* if something is too wrong.
*/
bool diagnosis(bool print) const;
/**
* Lets the user go through the states in the vector using gnuplot.
* Requires the file plotstate.gp in the directory to plot the states.
*/
void showHistory();
/**
* Returns the partons sorted in Strings.
*/
vector<DipoleState::String> strings();
/**
* Returns the active partons.
*/
list<PartonPtr> getPartons() const;
/**
* Returns the active dipoles.
*/
list<DipolePtr> getDipoles() const;
/**
* Returns all dipoles.
*/
set<DipolePtr> getAllDipoles() const {
return allDipoles;
}
/**
* Makes the state reabsorb a parton.
*/
void reabsorb(PartonPtr p);
/**
* Makes the state reabsorb a parton, and continues to absorb as
* long as the dipoles gets smaller along the colour chain.
*/
void recursiveReabsorb(PartonPtr p);
/**
* Forces one of the dipoles connected to the parton
* to immediately swing. But not with higher rapdity
* step than /a ymax1 of same colour, or /a ymax2 of
* another colour if same colour cant be found.
* /a ymax = 0 mean no limit.
*/
bool forceSwing(PartonPtr d, double ymax1, double ymax2);
/**
* Shuffles p+ and p- between the left and right part of the state
* so that the particles are really on shell. Assumes m=0 atm.
*/
void balanceMomenta();
/**
* Counts and returns the number of rapidity gaps in the state.
*/
int nRapidityGaps() const;
/**
* prints all the member variables to cout.
*/
void coutData();
/**
* Finds the maximum and minimum rapidity in a list of partons.
*/
static pair<double, double> rapidityRange(const list<PartonPtr> &);
/**
* checks if the first list has the lowest rapidity parton.
*/
static bool rapiditySorted(const list<PartonPtr> &, const list<PartonPtr> &);
/**
* Get the weight associated with the generation of this dipole state.
*/
inline double weight() const;
/**
* Set the weight associated with the generation of this dipole state.
*/
inline void weight(double);
/**
* Get the p+ that the original particle brought.
*/
inline Energy plus() const {
return thePlus;
}
/**
* Set the p+ that the original particle brought.
*/
inline void plus(Energy E) {
thePlus = E;
}
/**
* Get the p- that the original particle brought.
*/
inline Energy minus() const {
return theMinus;
}
/**
* Set the p- that the original particle brought.
*/
inline void minus(Energy E) {
theMinus = E;
}
/**
* Get the p- that the original particle was missing.
*/
inline Energy minusDeficit() const {
return theMinusDeficit;
}
/**
* Set if the state should save its history or not.
*/
inline void takeHistory(bool);
/**
* saves the current state into the states history, but only if
* doTakeHistory is true.
*/
inline void save();
/**
* The list of initial dipoles.
*/
inline const vector<DipolePtr> & initialDipoles() const;
/**
* Adds a dipole to the initial Dipoles.
* Mainly for testing purposes.
*/
inline void addDipole(Dipole & dip);
/**
* Get additional info about the wavefunction used to create this state.
*/
inline WFInfoPtr WFInfo() const;
/**
* Get additional info about the wavefunction used to create this state.
*/
inline const WaveFunction & wf() const;
/**
* Set additional info about the wavefunction used to create this state.
*/
inline void WFInfo(WFInfoPtr);
/**
* The controlling DipoleEventHandler.
*/
inline const DipoleEventHandler & handler() const;
/**
* Set the controlling DipoleEventHandler.
* added by CF to access from emitter.
*/
inline void handler(tcDipoleEventHandlerPtr);
/**
* Create a Dipole belonging to this state.
*/
inline tDipolePtr createDipole();
/**
* Generate a consistent colour index for the given Dipole.
*/
void generateColourIndex(tDipolePtr);
/**
* Returns all active dipoles of a certain colour.
*/
inline const vector<tDipolePtr> & swingCandidates(int) const;
/**
* adds the dipole, and all its children, to theSwingCandidates list.
* includes interacting dipoles.
*/
void sortDipoleFS(Dipole &);
/**
* Sort in colour all the dipoles originating from initialDipoles.
* includes interacting dipoles.
*/
void sortDipolesFS();
/**
* Sort in colour all the dipoles originating from initialDipoles.
* includes interacting dipoles.
*/
void sortFSDipoles();
/**
* adds the dipole, and all its children, to theSwingCandidates list.
*/
void sortDipole(Dipole &);
/**
* Sort in colour all the dipoles originating from initialDipoles.
*/
void sortDipoles();
/**
* Fakes a binary collision model with no fluctuation or saturation.
*/
void makeBinaryCollision(DipoleStatePtr dr, ImpactParameters b);
/**
* Fakes a glauber model.
*/
void makeGlauberCollision(DipoleStatePtr dr, ImpactParameters b);
/**
* Returns a list of all the active nonvalence partons that has not interacted.
*/
const list<PartonPtr> virtualPartons() const;
/**
* Calculates and returns the set of the loops in the state.
*/
const set< list<PartonPtr> > loops() const;
/**
* Returns the highest rapidity in any parton in the state.
*/
const double highestY() const;
/**
* Returns the highest rapidity in any parton in the state.
*/
const double lowestY() const;
/**
* Returns the rapidity the state has been evolved to.
*/
const double ymax() const;
/**
* Gets the energy the colliding state will supply.
*/
const Energy collidingEnergy() const;
/**
* Sets the energy the colliding state will supply.
*/
void collidingEnergy(Energy);
/**
* Prints all active dipoles to file.
*/
void printDipoles(ostream &) const;
/**
* Prints all active dipoles to file.
*/
void plotState(bool pause) const;
/**
* Prints info in a format fitting with the mathematica file for the movie.
*/
void printForMovie(double y, int count);
/**
* Just replaces the state with a single dipole of size 1/10. breaks energy cons etc.
*/
void turnIntoDipole();
/**
* Saves the state to a file.
*/
void saveGluonsToFile(double weight);
/**
* Given an output iterator fill the corresponding container with
* final, undecayed dipoles. This is done recursively, ie. if this
* dipole has decayed, the extract methods of the two children are
* called instead.
*/
template <typename OutputIterator>
inline void extract(OutputIterator it) 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();
public:
/**
* Clone this state and all the dipoles and partons in it.
*/
DipoleStatePtr clone();
protected:
/**
* The p+ and p- that the original particle had.
*/
Energy thePlus;
Energy theMinus;
/**
* Keeps track of extra missing p- that is not reflected in the
* valence partons.
*/
Energy theMinusDeficit;
/**
* The controlling DipoleEventHandler.
*/
tcDipoleEventHandlerPtr theHandler;
/**
* The list of initial dipoles.
*/
vector<DipolePtr> theInitialDipoles;
/**
* The active dipoles sorted by colour. /CF
*/
vector< vector<tDipolePtr> > theSwingCandidates;
/**
* Additional info about the wavefunction used to create this state.
*/
WFInfoPtr theWFInfo;
/**
* The weight associated with the generation of this dipole state.
*/
double theWeight;
/**
* If the history should be saved and displayed or not.
*/
bool doTakeHistory;
/**
* The rapidity the state has been evolved to.
*/
double theYmax;
/**
* How much energy the meeting particle have.
*/
Energy theCollidingEnergy;
/**
* The history of this state.
*/
vector<DipoleStatePtr> theHistory;
/**
* The set of all dipoles belonging to this state.
*/
set<DipolePtr> allDipoles;
protected:
/**
* Exception class for badly connected dipoles.
*/
struct DipoleConnectionException: public Exception {};
/**
* Exception class for bad kinematics.
*/
struct DipoleKinematicsException: public Exception {};
/**
* Exception class for bad DGLAP checks.
*/
struct DipoleDGLAPSafeException: public Exception {};
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DipoleState & operator=(const DipoleState &);
};
}
#include "DipoleState.icc"
#endif /* DIPSY_DipoleState_H */
diff --git a/DIPSY/EventFiller.cc b/DIPSY/EventFiller.cc
--- a/DIPSY/EventFiller.cc
+++ b/DIPSY/EventFiller.cc
@@ -1,1554 +1,1665 @@
// -*- 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 "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 <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) {}
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 {
ostream & log = CurrentGenerator::log();
static DebugItem printsteps("DIPSY::PrintSteps", 6);
static DebugItem timesteps("DIPSY::TimeSteps", 6);
if ( printsteps ) log << "<DIPSY> entering event filler." << endl;
lekplatsLeif(dl);
// Get a list of possible dipole-dipole interactions
FList fl = eh.xSecFn().flist(dl, dr, b);
if ( timesteps ) log << "<DIPSY> Extracted possible interactions in "
<< DipoleEventHandler::elapsed() << "s." << endl;
//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 ) {
dl.showHistory();
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();
if ( printsteps ) log << "<DIPSY> entering selectInteractions." << endl;
// Select the interactions which should be performed.
pair<RealPartonStatePtr, RealPartonStatePtr> realStates =
selectInteractions(fl, b, eh.xSecFn());
if ( timesteps ) log << "<DIPSY> Extracted final interactions in "
<< DipoleEventHandler::elapsed() << "s." << endl;
//Statistics for the number of interactions.
histNInt->fill(realStates.first->interactions.size(), currentWeight);
//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());
//Take some (diagnostic) statistics on the event after interactions
//are selected but before virtuals are removed.
//WARNING: This can take a lot of time in an AA event.
//takeBeforeStatistics(dl, dr);
//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);
//Discard event if no final state partons were found.
if ( strings.empty() || ! fillStep(step, inc, strings) ) {
weight = 0.0;
currentWeight = 0.0;
return weight;
}
//Some more statistics.
histBDist->fill(b.bVec().pt()*GeV, currentWeight);
histBDistCount->fill(b.bVec().pt()*GeV, 1.0);
histProbBDist->fill( b.bVec().pt()*GeV, weight );
return weight;
}
void EventFiller::takeBeforeStatistics(const DipoleState & dl, const DipoleState & dr) const {
list<PartonPtr> partonsL = dl.getPartons();
for ( list<PartonPtr>::const_iterator it = partonsL.begin() ; it != partonsL.end(); it++ ) {
histYDistBefore->fill( (*it)->y(), (*it)->pT().pt()/GeV*currentWeight );
if ( !((*it)->interacted()) )
continue;
histYDistRealBefore->fill((*it)->y(), (*it)->pT().pt()/GeV*currentWeight );
histYNumDistRealBefore->fill((*it)->y(), currentWeight );
if ( ((*it)->dipoles().first && (*it)->dipoles().first->interacted()) ||
((*it)->dipoles().second && (*it)->dipoles().second->interacted()) ) {
histYDistIntBefore->fill((*it)->y(), (*it)->pT().pt()/GeV*currentWeight );
histYNumDistIntBefore->fill( (*it)->y(), currentWeight );
}
if ( !((*it)->inInteractingLoop()) ) { //by far not the fastest way to do this.
histYDistLoopBefore->fill((*it)->y(), (*it)->pT().pt()/GeV*currentWeight );
histYNumDistLoopBefore->fill( (*it)->y(), currentWeight );
}
}
}
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;
histNuclParticipantDist->fill( untouchedNucleons, 1.0 );
histNuclColsBDist->fill( b*GeV, untouchedNucleons*currentWeight);
}
/*
* This is the probability that a parton is interacting, given that
* the state interacts, but the partons is not selected as primary.
*
* Please check calculations, but I think this is accurate. /CF
*
* 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 {
DipolePairVector 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<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));
}
//interaction probability (for at least one interaction)
//of the two states.
double totalP = Current<DipoleEventHandler>()->xSecFn().unitarize(sumfij);
if ( Debug::level > 5 ) cout << "1-exp(-F): " << totalP << ", sum(1-exp(-f)): " << sumUP << endl;
//create the real states.
RealPartonStatePtr rrs = new_ptr(RealPartonState());
RealPartonStatePtr lrs = new_ptr(RealPartonState());
if ( Debug::level > 5 ) cout << "created real states, first add valence and save." << endl;
//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();
//Decide if the event should be monitored by debugging.
//By having a monitoring probability proportional to the interaction
//probability, one get's a "minimum bias" sample of events when
//studying/debugging events through the graphical interface.
//This is a great way to understand what is going on in DIPSY.
bool enableMonitoring = false;
double maxWeight = 500.0;
if ( Debug::level > 5 && enableMonitoring && currentWeight > maxWeight*UseRandom::rnd() ) {
if ( currentWeight > maxWeight )
cout << "max weight topped (" << currentWeight
<< " > " << maxWeight << "), change in eventfiller!!" << endl;
cout << "monitor this event" << endl;
lrs->monitored = true;
rrs->monitored = true;
}
if ( Debug::level > 5 ) cout << "now select interactions" << endl;
DipolePairMap potential;
DipolePairMap failedPrims;
bool found = false;
int counter = 0;
while ( !found ) {
potential.clear();
//select a first interaction (since there has to be at least one).
FList::const_iterator prim = sel[UseRandom::rnd()];
if ( Debug::level > 5 ) cout << " selected prime interaction " << endl;
+ //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));
}
+
+ //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;
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 ) {
if ( Debug::level > 5 ) cout << " autofailed dipoles " << it->second->second.first << ", " << it->second->second.second << endl;
continue;
}
+
+ //Try to add the interaction.
if (addInteraction(it->second, lrs, rrs, interactions, b, xSec)) {
found = true;
}
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);
}
+ //We don't really use this mode any longer. Remove?
if ( mode() == 1 ) xSec.performInteractions(interactions, lrs, rrs, b);
histNCandInt->fill(potential.size(), currentWeight);
+ //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);
- } //ta bort, redan gjort?
+ }
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();
- // removeElasticInteractions(interactions, lrs, rrs);
-
+ //Debugging
if ( lrs->monitored ) {
cout << "found ints: " << lrs->interactions.size() << endl;
cout << "tried with " << potential.size() << " interactions" << endl;
lrs->plotBothStates(rrs, false);
- }
-
- if ( lrs->monitored && lrs->interactions.size() > 1 ) {
- cout << "final states:" << endl;
lrs->diagnosis(false);
rrs->diagnosis(false);
}
- else {
- lrs->diagnosis(false);
- rrs->diagnosis(false);
- }
-
+
+ //Diagnostics
for( DipolePairVector::const_iterator it = interactions.begin(); it != interactions.end(); it++ ) {
DipoleXSec::InteractionRecoil rec =
xSec.recoil((*it)->second.first->partons(), (*it)->second.second->partons(), b);
histYDistInt->fill((*it)->second.first->partons().first->y(),
rec.first.first.pt()/GeV*currentWeight);
histYDistInt->fill((*it)->second.first->partons().second->y(),
rec.first.second.pt()/GeV*currentWeight);
histYDistInt->fill(-(*it)->second.second->partons().first->y(),
rec.second.first.pt()/GeV*currentWeight);
histYDistInt->fill(-(*it)->second.second->partons().second->y(),
rec.second.second.pt()/GeV*currentWeight);
}
+ //Debugging and diagnostics.
nTried += potential.size()*currentWeight;
nInt += interactions.size()*currentWeight;
if (interactions.size() == 1) nInt1 += currentWeight;
if ( Debug::level > 5 ) {
cout << "done in select interaction, sum(tij) = " << sumUP << ", interactions are: " << endl;
for( DipolePairVector::const_iterator it = interactions.begin();
it != interactions.end(); it++ ) {
cout << " tij = " << (*it)->first.second << endl;
}
}
+ //return the real states.
return make_pair(lrs, rrs);
}
-void EventFiller::removeElasticInteractions(DipolePairVector interactions,
- RealPartonStatePtr lrs, RealPartonStatePtr rrs) const {
- DipolePairVector elastics;
- do {
- elastics.clear();
- for( DipolePairVector::iterator it = interactions.begin();
- it != interactions.end(); it++ ) {
- if ( (*it)->first.second > 2.0*UseRandom::rnd() )
- elastics.push_back(*it);
- }
- } while(elastics.size() == interactions.size());
- for( DipolePairVector::iterator it = interactions.begin();
- it != interactions.end(); it++ ) {
- lrs->interactions.remove((*it)->second.first);
- rrs->interactions.remove((*it)->second.second);
- }
- cout << "made " << elastics.size() << " elastic interactions out of "
- << interactions.size() << endl;
-}
-bool EventFiller::canInteract(Dipole & di, Dipole & dj) const {
- if ( di.interacted() || dj.interacted() ) return false;
- di.interact(dj);
- dj.interact(di);
- return true;
-}
-bool EventFiller::addInteraction(FList::const_iterator inter, RealPartonStatePtr lrs,
- RealPartonStatePtr rrs, DipolePairVector & inters,
- const ImpactParameters & b, const DipoleXSec & xSec) const {
+bool EventFiller::
+addInteraction(FList::const_iterator inter, RealPartonStatePtr lrs,
+ RealPartonStatePtr rrs, DipolePairVector & inters,
+ const ImpactParameters & b, const DipoleXSec & xSec) const {
+ //diagnostics
if ( Debug::level > 5 ) cout << "entering addInteraction" << endl;
if ( inters.size() > 0 )
nTestedInts += currentWeight;
if ( inter->second.first->interacted() || inter->second.second->interacted() )
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);
+ xSec.doesInt(inter->second.first->partons(),
+ inter->second.second->partons(), b);
+
+ //diagnostics
if ( Debug::level > 5 ) cout << " decided interations: " << doesInt.first.first << ", " << doesInt.first.second << "; " << doesInt.second.first << ", " << doesInt.second.second
<< "----------------------------------------------------" << endl;
- DipoleXSec::InteractionRecoil recs = xSec.recoil(inter->second.first->partons(),
- inter->second.second->partons(), b, doesInt);
- if ( !lrs->fullControlEvolution(inter->second.first, inter->second.second,
- doesInt.first.first, doesInt.first.second,
- recs.first.first.pt(), recs.first.second.pt()) ) {
+ //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()) ) {
+ //Just some diangostics.
histRealIters->fill( lrs->toCheck.size(), lrs->partonCalls);
histRealN->fill( lrs->toCheck.size(), 1.0);
if ( inters.size() > 0 ) {
failedEvo += currentWeight;
}
+
return false;
}
+
+ //diagnostics
if ( Debug::level > 5 ) cout << "found first evo!" << endl;
histRealIters->fill( lrs->toCheck.size(), lrs->partonCalls);
histRealN->fill( lrs->toCheck.size(), 1.0);
- if ( !rrs->fullControlEvolution(inter->second.second, inter->second.first,
- doesInt.second.first, doesInt.second.second,
- recs.second.first.pt(), recs.second.second.pt()) ) {
+
+ //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()) ) {
+ //diagnostics
histRealIters->fill( rrs->toCheck.size(), rrs->partonCalls);
histRealN->fill( rrs->toCheck.size(), 1.0);
+
+ //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);
+
+ //diagnostics.
if ( inters.size() > 0 ) {
failedEvo += currentWeight;
}
+
return false;
}
+
+ //diagnostics
if ( Debug::level > 5 ) cout << "found second evo!" << endl;
histRealIters->fill( rrs->toCheck.size(), rrs->partonCalls);
histRealN->fill( rrs->toCheck.size(), 1.0);
+
+ //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) ) {
+ //diagnostics
if ( Debug::level > 5 ) cout << "reverting with ints: " << inter->second.first << " and "
<< inter->second.second << endl;
+
+ //If failed, remove the interaction from the evolutions.
lrs->revertToPrevious(inter->second.first);
rrs->revertToPrevious(inter->second.second);
inters.pop_back();
+
+ //diagnostics
if ( inters.size() > 0 )
failedRec += currentWeight;
+
return false;
}
+
+ //If the interactions was ok, mark the dipoles as interacting
+ //and save the real states.
inter->second.first->interact(*inter->second.second);
inter->second.second->interact(*inter->second.first);
lrs->saveState();
rrs->saveState();
+
+ //diagnostics
if ( inters.size() > 1 )
foundInt += currentWeight;
if ( Debug::level > 5 )
cout << "++found consistent interaction :)" << endl;
+
return true;
}
vector<EventFiller::String>
EventFiller::extractStrings(DipoleState & dl, DipoleState & dr,
pair<RealPartonStatePtr, RealPartonStatePtr> realStates,
const ImpactParameters & b ) const {
-
+ //diagnostics
ostream & log = CurrentGenerator::log();
static DebugItem printsteps("DIPSY::PrintSteps", 6);
static DebugItem printstepdetails("DIPSY::PrintStepDetails", 6);
static DebugItem timesteps("DIPSY::TimeSteps", 6);
static DebugItem checkkinematics("DIPSY::CheckKinematics", 6);
if ( printsteps ) log << "<DIPSY> Entering extractStrings()" << endl;
+ //Just rename a bit for convenience.
DipoleStatePtr rightState = &dr;
DipoleStatePtr leftState = &dl;
-
RealPartonStatePtr lrs = realStates.first;
RealPartonStatePtr rrs = realStates.second;
+ //Check for errors.
if ( checkkinematics && lrs->checkForNegatives() )
Throw<SpaceLikeGluons>()
<< "negatives found at start of extractstrings!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<SpaceLikeGluons>()
<< "negatives found at start of extractstrings!" << Exception::warning;
- // if (Debug::level > 5 ) cout << "inmcoming:" << endl;
- // if (Debug::level > 5 ) lrs->plotBothStates(rrs, false);
- // if (Debug::level > 5 ) rightState->diagnosis(true);
if ( lrs->interactions.size() != rrs->interactions.size() )
cerr << "not same number of interactions in the right and left state!" << endl;
- //grow back the virtuals that are not part of interacting chain
+ //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);
+ //diagnostics
if ( printstepdetails ) {
log << "<DIPSY> first virtuals removed from states:" << endl;
lrs->plotBothStates(rrs, false);
}
+ //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();
+ //check for errors.
if ( checkkinematics && lrs->checkForNegatives() )
Throw<SpaceLikeGluons>()
<< "negatives found after merged virtuals!" << Exception::warning;
if ( checkkinematics && rrs->checkForNegatives() )
Throw<SpaceLikeGluons>()
<< "negatives found after merged virtuals!" << Exception::warning;
+ //Diagnostics
if ( printstepdetails ) {
log << "<DIPSY> merged virtuals in realstates, and saved them to states:"
<< endl;
lrs->plotBothStates(rrs, false);
}
- //let the recoils hit holes in the parents
+ //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);
- // lrs->diagnosis(false);
- // rrs->diagnosis(false);
-
+ //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);
+ //diagnostics
if ( printstepdetails ) {
log << "<DIPSY> joined states to one." << endl;
lrs->plotBothStates(rrs, false);
}
vector<pair<DipolePtr, DipolePtr> > swinging =
Current<DipoleEventHandler>()->xSecFn().getColourExchanges(lrs, rrs);
+
+ //diagnostics
if ( printstepdetails ) log << "<DIPSY> swinging " << swinging.size()
<< " interaction" << endl;
+
+ //swing the interacting dipoles.
for ( int i = 0; i < int(swinging.size()); i++ ) {
if ( printstepdetails )
log << "<DIPSY> calling reconnect nr " << i << endl;
Current<DipoleEventHandler>()->
xSecFn().reconnect(swinging[i].first, swinging[i].second);
}
+ //diagnostics
if ( printstepdetails ) {
log << "<DIPSY> swinged interactions" << endl;
lrs->plotBothStates(rrs, false);
}
- //reabsorb non_DGLAP pairs.
+ //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);
+ //diagnostics
if ( printsteps )
log << "<DIPSY> removed second round of virtuals" << endl;
- // if (Debug::level > 5 ) lrs->plotBothStates(rrs, false);
- //do colour reconnections
if ( timesteps ) log << "<DIPSY> Created real parton state in "
<< DipoleEventHandler::elapsed() << "s." << endl;
if ( printsteps ) log << "<DIPSY> starting final-state swing." << endl;
+
+ //do final state swings, aka pythias colour reconnection.
double yrange = FSSwingTime();
double ystep = FSSwingTimeStep();
for ( double y = 0.0; y < yrange; y += ystep ) {
if ( printstepdetails ) log << "<DIPSY> evolving FS from " << y
<< " to " << y+ystep << endl;
- // finalState->FSSwing(y, y + ystep);
finalState->swingFS(y, y + ystep);
}
+
+ //diagnostics
if ( printstepdetails ) log << "<DIPSY> done with FSswing" << endl;
if ( timesteps ) log << "<DIPSY> Done with final state swings in "
<< DipoleEventHandler::elapsed() << "s." << endl;
//compensate for the 6 valence charges in the proton
finalState->normaliseValenceCharge(theValenceChargeNormalisation);
+ //diagnostics
if ( printstepdetails ) log << "<DIPSY> normalised valence charge" << endl;
- //make sure the non-particpating nucleons (if AA) have the original colour flow.
+ //make sure the non-particpating nucleons (if AA) have the
+ //original colour flow.
finalState->restoreNonparticipants();
if ( printstepdetails ) log << "<DIPSY> restored nonparticipants." << endl;
- //try to find and fix any problems or inconsistencies that may have popped up.
+ //try to find and fix any problems or inconsistencies that may have
+ //popped up. Error messages are sent if anything is found.
dodgeErrors(finalState);
+ //More diagnostics
if ( printstepdetails ) log << "<DIPSY> done dodging errors." << endl;
-
if ( !finalState->diagnosis(false) && printsteps ) {
log << "<DIPSY> error in finalState!" << endl;
log << " left state diagnosis:" << endl;
lrs->diagnosis(true);
log << " right state diagnosis:" << endl;
rrs->diagnosis(true);
log << " both real states ploted" << endl;
lrs->plotBothStates(rrs, true);
}
else if ( printstepdetails ) log << "<DIPSY> passed error check." << endl;
+ //Prints the final state if in debugging mode. The weight trigger
+ //is good to get minimum bias-like events to study.
double weightmax = 500.0;
if ( currentWeight > weightmax*UseRandom::rnd() && Debug::level > 5 ) {
cout << "weight is " << currentWeight << " max is " << weightmax << endl;
cout << "all done!" << endl;
lrs->plotBothStates(rrs, false);
finalState->diagnosis(true);
}
- // finalState->saveGluonsToFile(currentWeight);
+ //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);
+ //diagnostics
nEvents += currentWeight;
pair<int, int> yesno = lrs->countYESNO();
int yes = yesno.first;
int no = yesno.second;
yesno = rrs->countYESNO();
yes += yesno.first;
no += yesno.second;
nYes += yes*currentWeight;
nNo += no*currentWeight;
avYInEvo += (lrs->avYInEvo()+rrs->avYInEvo())*currentWeight/2.0;
avYInInt += (finalState->avYInInt())*currentWeight;
- bool centralhighpt = false;
-
+ //more diagnostics ans statistics
nLoops += currentWeight*finalState->loops().size();
Energy scalarPT = 0.0*GeV;
Energy maxpT = 0.0*GeV;
list<PartonPtr> partons = finalState->getPartons();
for ( list<PartonPtr>::const_iterator it = partons.begin() ; it != partons.end(); it++ ) {
Energy pT = (*it)->pT().pt();
- if ( pT > 50*GeV && abs((*it)->y()) < 0.5 ) centralhighpt = true;
histYDist->fill( (*it)->y(), pT/GeV*currentWeight );
scalarPT += pT;
if ( pT > maxpT ) maxpT = pT;
histPTSpectrum->fill( pT/GeV, currentWeight);
if ( (*it)->rightMoving() ) {
histYDistDY->fill( (*it)->y(), ((*it)->y() - (*it)->oY())*currentWeight );
histY0DistDY->fill( (*it)->oY(), ((*it)->y() - (*it)->oY())*currentWeight );
histY0NumDist->fill( (*it)->oY(), currentWeight );
if ( (*it)->valence() )
histYDistDYVal->fill( (*it)->y(), ((*it)->y() - (*it)->oY())*currentWeight );
}
else {
histYDistDY->fill( (*it)->y(), -((*it)->y() + (*it)->oY())*currentWeight );
histY0DistDY->fill( -(*it)->oY(), -((*it)->y() + (*it)->oY())*currentWeight );
histY0NumDist->fill( -(*it)->oY(), currentWeight );
if ( (*it)->valence() )
histYDistDYVal->fill( (*it)->y(), -((*it)->y() + (*it)->oY())*currentWeight );
}
if ( (*it)->rightMoving() )
histYDistLeft->fill( (*it)->y(), pT/GeV*currentWeight );
else
histYDistRight->fill( (*it)->y(), pT/GeV*currentWeight );
if ( pT < 5.0*GeV )
histYDist0to5->fill( (*it)->y(), pT/GeV*currentWeight );
else if (pT < 20.0*GeV )
histYDist5to20->fill( (*it)->y(), pT/GeV*currentWeight );
else
histYDist20above->fill( (*it)->y(), pT/GeV*currentWeight );
histYNumDist->fill( (*it)->y(), currentWeight );
if ( (*it)->valence() )
histYNumDistVal->fill( (*it)->y(), currentWeight );
if ( ((*it)->dipoles().first && (*it)->dipoles().first->interacted()) ||
((*it)->dipoles().second && (*it)->dipoles().second->interacted()) ) {
histYNumDistInt->fill( (*it)->y(), currentWeight );
}
}
+ //diagnostics
list<DipolePtr> dips = finalState->getDipoles();
for ( list<DipolePtr>::const_iterator it = dips.begin() ; it != dips.end(); it++ ) {
histR->fill((*it)->size()*GeV, currentWeight);
}
for ( map<tPartonPtr,RealPartonPtr>::iterator it = lrs->partons.begin();
it != rrs->partons.end(); it++ ) {
if ( it == lrs->partons.end() ) it = rrs->partons.begin();
tPartonPtr p = it->first;
tRealPartonPtr rp = it->second;
if ( rp->keep != RealParton::YES ) continue;
for ( list< pair<tRealPartonPtr, pair<TransverseMomentum, Energy> > >::iterator jt = rp->recoils.begin();
jt != rp->recoils.end(); jt++ ) {
FillRecRange(p->y(), jt->first->theParton->y(), jt->second.first.pt());
- // if ( Debug::level > 5 ) cout << p->y() << " recoiled with " << jt->first->theParton->y() << " for pt " << jt->second.first.pt()/GeV << endl;
}
}
- // if ( Debug::level > 5 ) lrs->plotBothStates(rrs, true);
-
+
+ //debugging
if ( Debug::level > 5 && maxpT > 100.0*GeV ) {
cout << "paused since max pT is " << maxpT/GeV << endl;
lrs->plotBothStates(rrs, false);
finalState->diagnosis(true);
finalState->showHistory();
}
+ //diagnostics
histTotalScalarPT->fill( scalarPT/GeV, currentWeight );
//Sort the remaining partons in strings.
if ( printsteps )
log << "<DIPSY> Extracting final ThePEG particles." << endl;
vector<String> ret = finalState->strings();
if ( timesteps ) log << "<DIPSY> Done rest of extractStrings() in "
<< DipoleEventHandler::elapsed() << "s." << endl;
return ret;
}
void EventFiller::FillRecRange(double y1, double y2, Energy pt) const {
if ( abs(y1+0.5) < 0.5 )
histRecRange0to1->fill(y2, pt/GeV*currentWeight);
if ( abs(y2+0.5) < 0.5 )
histRecRange0to1->fill(y1, pt/GeV*currentWeight);
if ( abs(y1+2.5) < 0.5 )
histRecRange2to3->fill(y2, pt/GeV*currentWeight);
if ( abs(y2+2.5) < 0.5 )
histRecRange2to3->fill(y1, pt/GeV*currentWeight);
if ( abs(y1+5.0) < 1.0 )
histRecRange4to6->fill(y2, pt/GeV*currentWeight);
if ( abs(y2+5.0) < 1.0 )
histRecRange4to6->fill(y1, pt/GeV*currentWeight);
}
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 ( 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() ) 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
- for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = lrs->partons.begin();
- it != lrs->partons.end(); it++) {
+
+ //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++) {
+ 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"
<< "rightPlus: " << rightPlus/GeV
<< "\nrightMinus: " << rightMinus/GeV
<< "\nleftPlus: " << leftPlus/GeV
<< "\nleftMinus: " << leftMinus/GeV
<< "\nA: " << A
<< "\nB: " << B
<< Exception::warning;
return;
}
- //loop again, making left plus go to x*plus (and minus go to minus/x), and right minus go to y*minus (and plus go to plus/y)
+ //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::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 {
ostream & log = CurrentGenerator::log();
static DebugItem checkkinematics("DIPSY::CheckKinematics", 6);
static DebugItem printstepdetails("DIPSY::PrintStepDetails", 6);
- if ( mode() == 1 ) {
- return xSec.checkInteractions(sel, lrs, rrs, b);
- }
+ //We don't really use this mode any longer. Remove?
+ if ( mode() == 1 ) return xSec.checkInteractions(sel, lrs, rrs, b);
+ //Diagnostics
if ( printstepdetails ) log << "<DIPSY> now do interactions!" << endl;
- // if ( Debug::level > 5 ) lrs->plotBothStates(rrs, true);
- list<pair<bool, bool> >::const_iterator leftDoesInt = lrs->doesInts.begin();
- list<pair<bool, bool> >::const_iterator rightDoesInt = rrs->doesInts.begin();
+
+ //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++ ) {
+ //diagnostics
if ( checkkinematics ) {
log << "<DIPSY> summed momenta before interaction:" << endl;
lrs->checkPlusMinus();
rrs->checkPlusMinus();
}
+
+ //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;
- if ( printstepdetails )
- log << "<DIPSY> removed total recoil from left: ("
- << (recoil.first.first + recoil.first.second).x()/GeV << ", "
- << (recoil.first.first + recoil.first.second).y()/GeV << ")" << endl;
- if ( printstepdetails )
- log << "<DIPSY> removed total recoil from right: ("
- << (recoil.second.first + recoil.second.second).x()/GeV << ", "
- << (recoil.second.first + recoil.second.second).y()/GeV
- << ")" << endl;
+
+ //diagnostics
+ if ( printstepdetails )
+ log << "<DIPSY> removed total recoil from left: ("
+ << (recoil.first.first + recoil.first.second).x()/GeV << ", "
+ << (recoil.first.first + recoil.first.second).y()/GeV << ")" << endl;
+ if ( printstepdetails )
+ log << "<DIPSY> removed total recoil from right: ("
+ << (recoil.second.first + recoil.second.second).x()/GeV << ", "
+ << (recoil.second.first + recoil.second.second).y()/GeV
+ << ")" << endl;
if ( printstepdetails ) log << "<DIPSY> failed interaction!" << endl;
- // if ( Debug::level > 5 ) lrs->plotBothStates(rrs, true);
return false;
}
+
+ //diagnostics
if ( checkkinematics ) {
log << "<DIPSY> summed momenta after interaction:" << endl;
lrs->checkPlusMinus();
rrs->checkPlusMinus();
}
leftDoesInt++;
rightDoesInt++;
}
+
+ //diagnostics
if ( Debug::level > 5 ) cout << "done with interactions, success!" << endl;
- // if ( Debug::level > 5 ) lrs->plotBothStates(rrs, true);
if ( lrs->monitored )
cout << " passed!" << endl;
- // lrs->diagnosis(false);
- // rrs->diagnosis(false);
return true;
}
void EventFiller::removeVirtuals(DipoleStatePtr state) const {
+ //debugging
ostream & log = CurrentGenerator::log();
static DebugItem printstepdetails("DIPSY::PrintStepDeatils", 6);
list<PartonPtr> vP = state->getPartons();
+ //diagnostics
if ( printstepdetails )
log << "<DIPSY> checking " << vP.size() << " for virtuals." << endl;
+
+ //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 ) {
+ //diagnostics
if ( printstepdetails )
log << "<DIPSY> only one on shell in chain, swing." << endl;
+
+ //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);
}
}
+
+ //diagnostics
if ( printstepdetails ) log << "<DIPSY> done in remove virtuals" << endl;
}
// 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() {
if ( nEvents > 0 ) {
generator()->log()
<< "number of tried ints on passing try: " << nTried/nEvents << endl
<< "number of interactions: " << nInt/nEvents << endl
<< "fraction of single interactions: " << nInt1/nEvents << endl
<< "ratio of ints failed at evo: " << failedEvo/nTestedInts << endl
<< "ratio of ints failed at int: " << failedRec/nTestedInts << endl
<< "average YES and NO partons: " << nYes/nEvents
<< ", " << nNo/nEvents << endl
<< "average rapidity distance in evolution and int: " << avYInEvo/nEvents
<< ", " << avYInInt/nEvents << endl;
}
generator()->histogramFactory()->histogramFactory().
divide("/realItersNorm", *histRealIters, *histRealN);
generator()->histogramFactory()->histogramFactory().
divide("/probBDistNorm", *histProbBDist, *histBDistCount);
generator()->histogramFactory()->histogramFactory().
divide("/dipColsBDistNorm", *histDipColsBDist, *histBDist);
generator()->histogramFactory()->histogramFactory().
divide("/nuclColsBDistNorm", *histNuclColsBDist, *histBDist);
generator()->histogramFactory()->histogramFactory().
divide("/nuclCols2BDistNorm", *histNuclCols2BDist, *histBDist);
generator()->histogramFactory()->histogramFactory().
divide("/yDistDYNorm", *histYDistDY, *histYNumDist);
generator()->histogramFactory()->histogramFactory().
divide("/y0DistDYNorm", *histY0DistDY, *histY0NumDist);
generator()->histogramFactory()->histogramFactory().
divide("/yDistDYValNorm", *histYDistDYVal, *histYNumDistVal);
generator()->histogramFactory()->histogramFactory().
divide("/swingProb", *histSwings, *histEmissions);
histMaxPT->scale( generator()->histogramScale()/(millibarn*4.0) );
histInteractionPT->scale( generator()->histogramScale()/(millibarn*4.0) );
histMinR->scale( generator()->histogramScale()/(millibarn*4.0) );
histYDist->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDistDY->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histY0DistDY->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDistDYVal->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDistLeft->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDistRight->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDistBefore->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDistRealBefore->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDist0to5->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDist5to20->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDist20above->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDistInt->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDistIntBefore->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYDistLoopBefore->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYNumDistIntBefore->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYNumDistLoopBefore->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYNumDist->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histY0NumDist->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYNumDistVal->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYNumDistRealBefore->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histYNumDistInt->scale( generator()->histogramScale()/(generator()->integratedXSec()*0.2) );
histEvolutionPT->scale( generator()->histogramScale()/(millibarn*4.0) );
histTotalScalarPT->scale( generator()->histogramScale()/(millibarn*4.0) );
histNInt->scale( generator()->histogramScale()/generator()->integratedXSec() );
histNCandInt->scale( generator()->histogramScale()/generator()->integratedXSec() );
histSecTries->scale( generator()->histogramScale()/generator()->integratedXSec() );
histR->scale( generator()->histogramScale()/generator()->integratedXSec()*20.0 );
histRInt->scale( generator()->histogramScale()/generator()->integratedXSec()*20.0 );
histRIntDist->scale( generator()->histogramScale()/generator()->integratedXSec()*20.0 );
// Leifs lekplats
histRoverR0->scale( 20.0/generator()->N() );
histRoverR0MaxN->scale( 20.0/generator()->N() );
histRoverR0MinN->scale( 20.0/generator()->N() );
generator()->histogramFactory()->histogramFactory().
divide("/RmeanN", *histRavN, *histN);
generator()->histogramFactory()->histogramFactory().
divide("/PTmeanN", *histPTavN, *histN);
generator()->histogramFactory()->histogramFactory().
divide("/RlmeanN", *histRlavN, *histN);
// Slut på Leifs lekplats
}
void EventFiller::doinitrun() {
nEvents = 0.0;
nInt = 0.0;
nInt1 = 0.0;
nTried = 0.0;
nTestedInts = 0.0;
foundInt = 0.0;
failedEvo = 0.0;
rescatter = 0.0;
failedRec = 0.0;
overTenEvents = 0.0;
rejectedEvents = 0.0;
nLoops = 0.0;
nYes = 0.0;
nNo = 0.0;
avYInEvo = 0.0;
avYInInt = 0.0;
generator()->histogramFactory()->initrun();
histPTSpectrum = generator()->histogramFactory()->createHistogram1D
("PTSpectrum",200,0.0,100.0);
histSwings = generator()->histogramFactory()->createHistogram1D
("swings",110,-10.0,1.0);
histEmissions = generator()->histogramFactory()->createHistogram1D
("emissions",110,-10.0,1.0);
histRealIters = generator()->histogramFactory()->createHistogram1D
("realIters",101,-0.5,100.5);
histRealN = generator()->histogramFactory()->createHistogram1D
("realN",101,-0.5,100.5);
histNCandInt = generator()->histogramFactory()->createHistogram1D
("nCandInt",100,-0.5,99.5);
histSecTries = generator()->histogramFactory()->createHistogram1D
("secTries",100,-0.5,99.5);
histRecRange0to1 = generator()->histogramFactory()->createHistogram1D
("recRange0to1",200,-10.0,10.0);
histRecRange1to2 = generator()->histogramFactory()->createHistogram1D
("recRange1to2",200,-10.0,10.0);
histRecRange2to3 = generator()->histogramFactory()->createHistogram1D
("recRange2to3",200,-10.0,10.0);
histRecRange3to4 = generator()->histogramFactory()->createHistogram1D
("recRange3to4",200,-10.0,10.0);
histRecRange4to6 = generator()->histogramFactory()->createHistogram1D
("recRange4to6",200,-10.0,10.0);
histRecRange6to8 = generator()->histogramFactory()->createHistogram1D
("recRange6to8",200,-10.0,10.0);
histBDist = generator()->histogramFactory()->createHistogram1D
("bDist",200,0.0,200.0);
histBDistCount = generator()->histogramFactory()->createHistogram1D
("bDistCount",200,0.0,200.0);
histProbBDist = generator()->histogramFactory()->createHistogram1D
("probBDist",200,0.0,200.0);
histDipColsBDist = generator()->histogramFactory()->createHistogram1D
("dipColsBDist",200,0.0,200.0);
histNuclColsBDist = generator()->histogramFactory()->createHistogram1D
("nuclColsBDist",200,0.0,200.0);
histNuclCols2BDist = generator()->histogramFactory()->createHistogram1D
("nuclCols2BDist",200,0.0,200.0);
histNuclParticipantDist = generator()->histogramFactory()->createHistogram1D
("nuclParticipantDist",400,-0.5,399.5);
histYShift = generator()->histogramFactory()->createHistogram1D
("YShift",100,0.0,1.0);
histInteractionPT = generator()->histogramFactory()->createHistogram1D
("interactionPT",50,0.0,200.0);
histMaxPT = generator()->histogramFactory()->createHistogram1D
("maxPT",50,0.0,200.0);
histMinR = generator()->histogramFactory()->createHistogram1D
("minR",50,0.0,200.0);
histYDist = generator()->histogramFactory()->createHistogram1D
("yDist",100,-10.0,10.0);
histYDistDY = generator()->histogramFactory()->createHistogram1D
("yDistDY",100,-10.0,10.0);
histY0DistDY = generator()->histogramFactory()->createHistogram1D
("yDistDY0",100,-10.0,10.0);
histYDistDYVal = generator()->histogramFactory()->createHistogram1D
("yDistDYVal",100,-10.0,10.0);
histYDistLeft = generator()->histogramFactory()->createHistogram1D
("yDistLeft",100,-10.0,10.0);
histYDistRight = generator()->histogramFactory()->createHistogram1D
("yDistRight",100,-10.0,10.0);
histYDistBefore = generator()->histogramFactory()->createHistogram1D
("yDistBefore",100,-10.0,10.0);
histYDistRealBefore = generator()->histogramFactory()->createHistogram1D
("yDistRealBefore",100,-10.0,10.0);
histYDist0to5 = generator()->histogramFactory()->createHistogram1D
("yDist0to5",100,-10.0,10.0);
histYDist5to20 = generator()->histogramFactory()->createHistogram1D
("yDist5to20",100,-10.0,10.0);
histYDist20above = generator()->histogramFactory()->createHistogram1D
("yDist20above",100,-10.0,10.0);
histYDistInt = generator()->histogramFactory()->createHistogram1D
("yDistInt",100,-10.0,10.0);
histYDistIntBefore = generator()->histogramFactory()->createHistogram1D
("yDistIntBefore",100,-10.0,10.0);
histYDistLoopBefore = generator()->histogramFactory()->createHistogram1D
("yDistLoopBefore",100,-10.0,10.0);
histYNumDistIntBefore = generator()->histogramFactory()->createHistogram1D
("yNumDistIntBefore",100,-10.0,10.0);
histYNumDistLoopBefore = generator()->histogramFactory()->createHistogram1D
("yNumDistLoopBefore",100,-10.0,10.0);
histYNumDist = generator()->histogramFactory()->createHistogram1D
("yNumDist",100,-10.0,10.0);
histY0NumDist = generator()->histogramFactory()->createHistogram1D
("y0NumDist",100,-10.0,10.0);
histYNumDistVal = generator()->histogramFactory()->createHistogram1D
("yNumDistVal",100,-10.0,10.0);
histYNumDistRealBefore = generator()->histogramFactory()->createHistogram1D
("yNumDistRealBefore",100,-10.0,10.0);
histYNumDistInt = generator()->histogramFactory()->createHistogram1D
("yNumDistInt",100,-10.0,10.0);
histEvolutionPT = generator()->histogramFactory()->createHistogram1D
("evolutionPT",50,0.0,200.0);
histTotalScalarPT = generator()->histogramFactory()->createHistogram1D
("totalScalarPT",50,0.0,200.0);
histNInt = generator()->histogramFactory()->createHistogram1D
("nInt",51.0,-0.5,50.5);
histR = generator()->histogramFactory()->createHistogram1D
("r",200,0.0,10.0);
histRInt = generator()->histogramFactory()->createHistogram1D
("rInt",200,0.0,10.0);
histRIntDist = generator()->histogramFactory()->createHistogram1D
("rIntDist",200,0.0,10.0);
histr32 = generator()->histogramFactory()->createHistogram1D
("r32",200,0.0,10.0);
histrint2 = generator()->histogramFactory()->createHistogram1D
("rint2",200,0.0,10.0);
histfij = generator()->histogramFactory()->createHistogram1D
("fij",200,-0.0025,0.9975);
// Leifs lekplats
histRoverR0 = generator()->histogramFactory()->createHistogram1D
("RoverR0",100,0.0,5.0);
histRoverR0MaxN = generator()->histogramFactory()->createHistogram1D
("RoverR0MaxN",100,0.0,5.0);
histRoverR0MinN = generator()->histogramFactory()->createHistogram1D
("RoverR0MinN",100,0.0,5.0);
histRavN = generator()->histogramFactory()->createHistogram1D
("RavN",100,0.0,100.0);
histN = generator()->histogramFactory()->createHistogram1D
("RavNref",100,0.0,100.0);
histPTavN = generator()->histogramFactory()->createHistogram1D
("PTavN",100,0.0,100.0);
histRlavN = generator()->histogramFactory()->createHistogram1D
("RlavN",100,0.0,100.0);
// Slut på Leifs lekplats
HandlerBase::doinitrun();
}
void EventFiller::lekplatsLeif(DipoleState & ds) const {
list<DipolePtr> dips = ds.getDipoles();
InvEnergy r0 = ds.WFInfo()->r();
InvEnergy rsum = 0.0*InvGeV;
Energy ptsum = 0.0*GeV;
double rprod = 1.0;
for ( list<DipolePtr>::iterator it = dips.begin(); it != dips.end(); ++it ) {
histRoverR0->fill((**it).size()/r0, 1.0);
if ( dips.size() > 40 )
histRoverR0MaxN->fill((**it).size()/r0, 1.0);
else
histRoverR0MinN->fill((**it).size()/r0, 1.0);
rsum += (**it).size();
ptsum += 1.0/(**it).size();
rprod *= (**it).size()*GeV;
}
rsum /= dips.size();
ptsum /= dips.size();
rprod = pow(rprod, 1.0/dips.size());
histRlavN->fill(dips.size(), rprod);
histRavN->fill(dips.size(), rsum/InvGeV);
histPTavN->fill(dips.size(), ptsum/GeV);
histN->fill(dips.size(), 1.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;
}
LorentzMomentum p1save = str[i1]->momentum();
LorentzMomentum p2save = str[i2]->momentum();
LorentzMomentum p3save = str[i3]->momentum();
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);
}
void EventFiller::persistentInput(PersistentIStream & is, int) {
is >> theAbsorber >> theRecoilScheme >> theMode >> theSingleMother
>> theDGLAPinPT >> theEffectiveWeights >> theFSSwingTime
>> theFSSwingTimeStep >> theValenceChargeNormalisation
>> iunit(thePTCut, GeV);
}
DescribeClass<EventFiller,HandlerBase>
describeDIPSYEventFiller("DIPSY::EventFiller", "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. Reccomended for pp and anything except heavy ions. Possibly high energy pp can be a bit slow.",
+ "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.",
+ "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.",
+ "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 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);
}
diff --git a/DIPSY/EventFiller.h b/DIPSY/EventFiller.h
--- a/DIPSY/EventFiller.h
+++ b/DIPSY/EventFiller.h
@@ -1,569 +1,561 @@
// -*- 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;
/**
* 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;
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;
/**
* 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;
- /**
- * Takes the interactions in interact and removes some
- * of them from the interactions in the real parton states.
- **/
- virtual void removeElasticInteractions(DipolePairVector interactions,
- RealPartonStatePtr lrs, RealPartonStatePtr rrs) const;
- /**
- * Check if two dipoles can interact and, in that case, mark them as
- * interacted with eachother.
- */
- virtual bool canInteract(Dipole & di, Dipole & dj) 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 ;
/**
* Extract all strings after according to the given interactions.
*/
virtual vector<String>
extractStrings(DipoleState & dl, DipoleState & dr,
pair<RealPartonStatePtr, RealPartonStatePtr>, 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;
/**
* 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;
}
/**
* Take statistics on an emission that has swinged or not before emission.
*/
inline void swingedEmission(double y) const {
histSwings->fill(y,1.0);
histEmissions->fill(y,1.0);
}
inline void nonSwingedEmission(double y) const {
histEmissions->fill(y,1.0);
}
//@}
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:
+ //All histograms diagnostics. Some educational information, but nothing
+ //that can't easily be put back in, so feel free to remove.
FactoryBase::tH1DPtr histSwings;
FactoryBase::tH1DPtr histEmissions;
FactoryBase::tH1DPtr histRealIters;
FactoryBase::tH1DPtr histRealN;
FactoryBase::tH1DPtr histNCandInt;
FactoryBase::tH1DPtr histSecTries;
FactoryBase::tH1DPtr histRecRange0to1;
FactoryBase::tH1DPtr histRecRange1to2;
FactoryBase::tH1DPtr histRecRange2to3;
FactoryBase::tH1DPtr histRecRange3to4;
FactoryBase::tH1DPtr histRecRange4to6;
FactoryBase::tH1DPtr histRecRange6to8;
FactoryBase::tH1DPtr histBDist;
FactoryBase::tH1DPtr histBDistCount;
FactoryBase::tH1DPtr histProbBDist;
FactoryBase::tH1DPtr histDipColsBDist;
FactoryBase::tH1DPtr histNuclColsBDist;
FactoryBase::tH1DPtr histNuclCols2BDist;
FactoryBase::tH1DPtr histNuclParticipantDist;
FactoryBase::tH1DPtr histYShift;
FactoryBase::tH1DPtr histInteractionPT;
FactoryBase::tH1DPtr histMaxPT;
FactoryBase::tH1DPtr histMinR;
FactoryBase::tH1DPtr histPTSpectrum;
FactoryBase::tH1DPtr histYDist;
FactoryBase::tH1DPtr histYDistDY;
FactoryBase::tH1DPtr histYDistDYVal;
FactoryBase::tH1DPtr histY0DistDY;
FactoryBase::tH1DPtr histYDistLeft;
FactoryBase::tH1DPtr histYDistRight;
FactoryBase::tH1DPtr histYDistBefore;
FactoryBase::tH1DPtr histYDistRealBefore;
FactoryBase::tH1DPtr histYDist0to5;
FactoryBase::tH1DPtr histYDist5to20;
FactoryBase::tH1DPtr histYDist20above;
FactoryBase::tH1DPtr histYDistInt;
FactoryBase::tH1DPtr histYDistIntBefore;
FactoryBase::tH1DPtr histYDistLoopBefore;
FactoryBase::tH1DPtr histYNumDistIntBefore;
FactoryBase::tH1DPtr histYNumDistLoopBefore;
FactoryBase::tH1DPtr histYNumDist;
FactoryBase::tH1DPtr histY0NumDist;
FactoryBase::tH1DPtr histYNumDistVal;
FactoryBase::tH1DPtr histYNumDistRealBefore;
FactoryBase::tH1DPtr histYNumDistInt;
FactoryBase::tH1DPtr histEvolutionPT;
FactoryBase::tH1DPtr histTotalScalarPT;
FactoryBase::tH1DPtr histNInt;
FactoryBase::tH1DPtr histR;
FactoryBase::tH1DPtr histRInt;
FactoryBase::tH1DPtr histRIntDist;
FactoryBase::tH1DPtr histr32;
FactoryBase::tH1DPtr histrint2;
FactoryBase::tH1DPtr histfij;
// Leifs lekplats
FactoryBase::tH1DPtr histRoverR0;
FactoryBase::tH1DPtr histRoverR0MaxN;
FactoryBase::tH1DPtr histRoverR0MinN;
FactoryBase::tH1DPtr histRavN;
FactoryBase::tH1DPtr histN;
FactoryBase::tH1DPtr histPTavN;
FactoryBase::tH1DPtr histRlavN;
void lekplatsLeif(DipoleState &) const;
// Slut på Leifs lekplats
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;
/**
* 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;
/**
* Takes statistics on the states before the absorbtion or collision.
*/
void takeBeforeStatistics(const DipoleState & dl, const DipoleState & dr) const;
/**
* Take statistics on a recoil of pt between particles at y1 and y2.
**/
void FillRecRange(double y1, double y2, Energy pt) 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;
public:
/**
* Exception class for space-like gluon momenta.
*/
struct SpaceLikeGluons: 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/LHCpPb.in b/DIPSY/LHCpPb.in
new file mode 100644
--- /dev/null
+++ b/DIPSY/LHCpPb.in
@@ -0,0 +1,72 @@
+cd /DIPSY
+
+#To get output to a file, uncomment
+#finalState->saveGluonsToFile(currentWeight)
+#in EventFiller.cc. The events are saved to files
+#~/DIPSYevents/filename/eventN.dat
+#Where the directory ~/DIPSYevents/filename must exist before the run.
+#Otherwise the directory can be set in DipoleState::saveGluonsToFile
+#Filename is the name of the runfile, so LHCpPb or LHCpPbB01 in
+#this case.
+
+#Probably want to switch off FSR and hadronisation as well, as
+#the event is saved before that.
+
+set /Defaults/Particles/u:NominalMass 0.33
+set /Defaults/Particles/d:NominalMass 0.33
+set /Defaults/Particles/s:NominalMass 0.5
+set /Defaults/Particles/c:NominalMass 1.5
+set /Defaults/Particles/b:NominalMass 4.8
+
+read CurrentTune.in
+
+read RivetAnalyses.in
+
+create ThePEG::ProgressLog Logger ProgressLog.so
+set Logger:Interval 600
+
+set LHCpPbEventHandler:PreSamples 0
+set LHCpPbEventHandler:ConsistencyLevel 0
+set LHCpPbEventHandler:DecayHandler:MaxLifeTime 10
+set LHCpPbEventHandler:XSecFn:CheckOffShell false
+insert LHCpPbGenerator:AnalysisHandlers[0] Logger
+#insert LHCpPbGenerator:AnalysisHandlers[0] HIAna
+set LHCGenerator:EventHandler:EventFiller:PTCut 0.6
+
+#Not sure what CoM energy they will run at. 40TeV is about same
+#nucleon-nucleon as PbPb.
+set LHCpPbLumi:Energy 400.0
+
+#Run with fix impact parameter, like in AA.
+#This generates events in a ring. The weight
+#(with interaction probability of 1) will automatically
+#be set to 1, 3, 5, 7 in the B-ranges 0-1, 1-2, 2-3, etc, so
+#nothing should have to be micromanaged with the weights,
+#as long as the events are run in the below B rings.
+#With equal amount of events in all rings (out to where int prob
+#is very low), one can then find expectation values by average over
+#all events, with the provided weights.
+create DIPSY::FixedImpactGenerator FixedB FixedImpactGenerator.so
+set LHCpPbGenerator:EventHandler:BGen FixedB
+
+#AA has 1000 events per ring.
+#pPb will run much faster, but maybe not take that much less memory.
+#See how much time/space you got.
+set LHCpPbGenerator:NumberOfEvents 1000
+
+set FixedB:MinB = 0.0
+set FixedB:MaxB = 1.0
+saverun LHCpPbB01 LHCpPbGenerator
+
+set FixedB:MinB = 1.0
+set FixedB:MaxB = 2.0
+saverun LHCpPbB12 LHCpPbGenerator
+
+set FixedB:MinB = 2.0
+set FixedB:MaxB = 3.0
+saverun LHCpPbB12 LHCpPbGenerator
+
+#etc.
+#I ran Au-Au up to B=15, and PbPb to 19.
+#So up to 15 should be enough here.
+#See where the weight start dropping of too much.
\ No newline at end of file
diff --git a/DIPSY/RealPartonState.cc b/DIPSY/RealPartonState.cc
--- a/DIPSY/RealPartonState.cc
+++ b/DIPSY/RealPartonState.cc
@@ -1,2270 +1,2308 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the RealPartonState class.
//
#include "RealPartonState.h"
#include "RealParton.h"
#include "Dipole.h"
#include "DipoleState.h"
#include "DipoleEventHandler.h"
#include "ThePEG/Utilities/Current.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/DebugItem.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include <iostream>
using namespace DIPSY;
RealPartonState::RealPartonState(): consistent(true), calls(0) {}
RealPartonState::~RealPartonState() {}
void RealPartonState::newInteraction(tDipolePtr intDip, tDipolePtr otherIntDip,
bool firstInt, bool secondInt,
Energy rec1, Energy rec2) {
+ //The ancestors of the interacting partons, ie the partons that have
+ //been involved in making the interacting dipole.
RealPartonSet ancestors;
+
+ //diagnostics
if ( Debug::level > 5 )
cout << "adding interaction " << intDip << " to " << this << endl;
+
+ //clear the partons to be checked in this interaction.
+ //Only partons involved in this interaction are checked,
+ //ther others are assumed to still be ok since last check.
toCheck.clear();
+
+ //Add interacting dipoles and which partons are interacting.
interactions.push_back(intDip);
doesInts.push_back(make_pair(firstInt, secondInt));
+
+ //Set up the real partons.
+ //Some settings only add one of partons in the interacting dipole,
+ //hence all the "if (firstInt )" etc.
+ //Note that addParton recursively calls the parents as well,
+ //so all new realparton instances needed are created.
RealPartonPtr rp1, rp2;
if ( firstInt ) {
if ( Debug::level > 5 ) cout << "calling first addparton for interacting parton " << intDip->partons().first->oY() << endl;
addParton(intDip->partons().first);
rp1 = getReal(intDip->partons().first);
}
if ( secondInt ) {
if ( Debug::level > 5 ) cout << "calling second addparton for interacting parton " << intDip->partons().second->oY() << endl;
addParton(intDip->partons().second);
rp2 = getReal(intDip->partons().second);
}
+
+ //Find the interaction distance of the interacting partons.
+ //This is needed when checking pT-ordering (or almost equivalently dipole-
+ //size ordering) later, and the last-before-interaction dipole needs
+ //to be checked. The picture here is to compare along the colour chains
+ //AFTER the interacting dipole has swinged, which is a more frame-
+ //independent approach, as otherwise we would have to compare the
+ //last-before-interacting dipole with a dipole that is not present in
+ //the final state.
+ //Also, since the interacting dipole has (roughly) the weight of a
+ //dipole that has emitted, it would mess up the weights if it entered
+ //the pT-max reweighting scheme.
if ( rp1 ) {
if ( !rp1->interacting ) {
if ( rec1 == ZERO ) rec1 = rec2;
if ( rec1 == ZERO ) cerr << "zero recoil in realpartonstate::newinteraction" << endl;
rp1->intDist = rp1->theParton->pTScale()/rec1;
}
else {
if ( rec1 == ZERO ) rec1 = rec2;
if ( rec1 == ZERO ) cerr << "zero recoil in realpartonstate::newinteraction" << endl;
rp1->intDist = min(rp1->intDist, rp1->theParton->pTScale()/rec1);
}
+
+ //book keeping.
rp1->interacting = intDip;
rp1->secondInt = true;
RealPartonSet a1 = rp1->ancestors();
ancestors.insert(a1.begin(), a1.end());
- // rp1->addInteraction(intDip);
}
if ( rp2 ) {
if ( !rp2->interacting ) {
if ( rec2 == ZERO ) rec2 = rec1;
if ( rec2 == ZERO ) cerr << "zero recoil in realpartonstate::newinteraction" << endl;
rp2->intDist = rp2->theParton->pTScale()/rec2;
}
else {
if ( rec2 == ZERO ) rec2 = rec1;
if ( rec2 == ZERO ) cerr << "zero recoil in realpartonstate::newinteraction" << endl;
rp2->intDist = min(rp2->intDist, rp2->theParton->pTScale()/rec2);
}
+
+ //book keeping for future referens
rp2->interacting = intDip;
rp2->firstInt = true;
RealPartonSet a2 = rp2->ancestors();
ancestors.insert(a2.begin(), a2.end());
- // rp2->addInteraction(intDip);
}
+ //and make all the ancestors remember that they contributed.
for ( RealPartonSet::iterator it = ancestors.begin(); it != ancestors.end(); it++ ) {
(*it)->interactions.insert(intDip);
}
}
tRealPartonPtr RealPartonState::getReal(tPartonPtr p) {
+ //If the parton already exists in the real state, just retur
+ //the real parton.
map<tPartonPtr,RealPartonPtr>::iterator it = partons.lower_bound(p);
- if ( it->first == p ) {
- return it->second;
- }
+ if ( it->first == p ) return it->second;
+
+ //Otherwise create a new realparton for the parton.
RealPartonPtr rp = new_ptr(RealParton(p));
rp->realState = this;
+
+ //fluct -1 means that they are not a fluctuation,
+ //ie they are not (yet) set up to be merged with an other partons.
rp->fluct = -1;
rp->cFluct = -1;
+
+ //insert in the right place in the ordered map of parton -> realparton.
if ( it != partons.begin() ) partons.insert(--it, make_pair(p, rp));
else partons.insert(make_pair(p, rp));
return rp;
}
tRealPartonPtr RealPartonState::addParton(tPartonPtr p) {
if ( !p ) return RealPartonPtr();
if ( partons.lower_bound(p)->first == p ) return getReal(p);
RealPartonPtr rp = getReal(p);
rp->nMothers = 2;
if ( !(p->valence()) && rp->interactions.size() == 0 ) {
if ( Current<DipoleEventHandler>()->eventFiller().singleMother() &&
!p->swingedEmission() ) {
rp->nMothers = 1;
double P1 = p->dist2(*p->parents().second)/
(p->dist2(*p->parents().second) + p->dist2(*p->parents().first));
if ( P1 > UseRandom::rnd() ) {//pick feynman diagram
rp->setOMother(addParton(p->parents().first));
}
else
rp->setOMother(addParton(p->parents().second));
}
else {
rp->setFirstOMother(addParton(p->parents().first));
rp->setSecondOMother(addParton(p->parents().second));
}
}
if ( p->valence() && rp->interactions.size() == 0 ) {
valence.insert(rp);
rp->setYES();
}
else toCheck.insert(rp);
return rp;
}
void RealPartonState::addValence(DipoleState & state) {
if ( Debug::level > 5 ) cout << "entering addvalence" << endl;
plus = state.plus();
minus = state.minus();
if ( Debug::level > 5 ) cout << " total plus and minus set to " << plus/GeV << ", " << minus/GeV << endl;
monitored = false;
totalRecoil = TransverseMomentum(ZERO, ZERO);
if ( Debug::level > 5 ) cout << " total recoil set to 0" << endl;
vector<DipolePtr> valenceDip = state.initialDipoles();
for ( int i = 0; i < int(valenceDip.size()); i++) {
tRealPartonPtr p1 = addParton(valenceDip[i]->partons().first);
tRealPartonPtr p2 = addParton(valenceDip[i]->partons().second);
p1->cValence = true;
p1->cValencePlus = p1->theParton->valencePlus();
p1->cValencePT = p1->theParton->valencePT();
p2->cValence = true;
p2->cValencePlus = p2->theParton->valencePlus();
p2->cValencePT = p2->theParton->valencePT();
oValence.insert(p1);
oValence.insert(p2);
}
}
void RealPartonState::merge(RealPartonPtr rp1, RealPartonPtr rp2) {
if ( Debug::level > 5 )
cout << "merging partons at y = " << rp1->oY() << ", " << rp2->oY() << " or ("
<< rp1->theParton->position().x()*GeV << ", " << rp1->theParton->position().y()*GeV << "), ("
<< rp2->theParton->position().x()*GeV << ", " << rp2->theParton->position().y()*GeV << ")\n";
if ( rp1->fluct != -1 && rp2->fluct != -1 ) {
//merge flucts if different (can they already be in same?)
if ( rp1->fluct != rp2->fluct ) {
flucts[rp1->fluct].insert(flucts[rp2->fluct].begin(), flucts[rp2->fluct].end());
int old = rp2->fluct;
for ( RealPartonSet::iterator it = flucts[old].begin(); it != flucts[old].end(); it++ )
(*it)->fluct = rp1->fluct;
flucts[old].clear();
}
}
else if ( rp1->fluct != -1 ) {
rp2->fluct = rp1->fluct;
flucts[rp1->fluct].insert(rp2);
}
else if ( rp2->fluct != -1 ) {
rp1->fluct = rp2->fluct;
flucts[rp2->fluct].insert(rp1);
}
else {
flucts.push_back(RealPartonSet());
int i = flucts.size() - 1;
flucts[i].insert(rp1);
flucts[i].insert(rp2);
rp1->fluct = i;
rp2->fluct = i;
}
// redistributePlus(flucts[rp1->fluct]);
//dela upp p+ och pT jamnt mellan de virtuella
makeCollinear(flucts[rp1->fluct]);
if ( Debug::level > 5 ) cout << " done merging" << endl;
}
void RealPartonState::makeCollinear(const RealPartonSet & rps) {
if ( rps.empty() ) return;
Energy avPlus = ZERO;
TransverseMomentum avPT = TransverseMomentum();
bool interacting = false;
bool valence = false;
for ( RealPartonSet::const_iterator it = rps.begin(); it != rps.end(); it++ ) {
if ( Debug::level > 5 )
cout << " from " << (*it)->oY() << ": (" << (*it)->pT.x()/GeV << ", "
<< (*it)->pT.y()/GeV << ")\n";
if ( (*it)->interacting ) interacting = true;
if ( (*it)->theParton->valence() ) valence = true;
avPlus += (*it)->plus;
avPT += (*it)->pT;
}
if ( Debug::level > 5 ) cout << " total (" << avPT.x()/GeV << ", " << avPT.y()/GeV << ")\n";
avPlus /= rps.size();
avPT /= rps.size();
if ( Debug::level > 5 ) cout << " each get (" << avPT.x()/GeV << ", " << avPT.y()/GeV << ")\n";
for ( RealPartonSet::const_iterator it = rps.begin();
it != rps.end(); it++ ) {
(*it)->plus = avPlus;
(*it)->pT = avPT;
(*it)->updateYMinus();
}
}
void RealPartonState::redistributePlus(const RealPartonSet & rps) {
if ( rps.empty() ) return;
Energy totalPlus = ZERO;
double totalWeight = ZERO;
for ( RealPartonSet::const_iterator it = rps.begin();
it != rps.end(); it++ ) {
totalWeight += exp(-(*it)->oY());
totalPlus += (*it)->plus;
}
for ( RealPartonSet::const_iterator it = rps.begin();
it != rps.end(); it++ ) {
(*it)->plus = totalPlus*exp(-(*it)->oY())/totalWeight;
(*it)->updateYMinus();
// cout << " redistributed plus at " << (*it)->oY() << " to " << (*it)->plus/GeV
// << ", relative y is " << -log(exp(-(*it)->oY())/totalWeight) << endl;
}
}
void RealPartonState::splitFluct(RealPartonPtr rp1, RealPartonPtr rp2) {
// cout << "entered splitFluct with rps " << rp1->oY() << ", " << rp2->oY() << "!\n";
// diagnosis(true);
if ( rp1->fluct == -1 || rp2->fluct == -1 ) {
// cout << "asked to split a fluct where the partons werent in a fluct" << endl;
return;
}
else if ( rp1->fluct != rp2->fluct ) {
// cout << "asked to split a fluct that is already split" << endl;
return;
}
int i = rp1->fluct;
//create two new flucts with rp1 and rp2
flucts.push_back(RealPartonSet());
int i1 = flucts.size() - 1;
flucts[i1].insert(rp1);
flucts.push_back(RealPartonSet());
int i2 = flucts.size() - 1;
flucts[i2].insert(rp2);
// cout << " created new flucts " << i1 << " and " << i2 << ", old was " << i << endl;
//move all other partons to the fluct belonging to the closest parton.
for ( RealPartonSet::iterator it = flucts[i].begin(); it != flucts[i].end(); it++) {
if ((*it)->theParton->dist2(*rp1->theParton) < (*it)->theParton->dist2(*rp2->theParton) ) {
// cout << " move " << (*it)->oY() << " to " << i1 << endl;
(*it)->fluct = i1;
flucts[i1].insert(*it);
}
else {
// cout << " move " << (*it)->oY() << " to " << i2 << endl;
(*it)->fluct = i2;
flucts[i2].insert(*it);
}
}
// cout << " empty old fluct" << endl;
flucts[i].clear();
//find highest oY parton.
RealPartonPtr firstRP = (rp1->oY() < rp2->oY() ) ? rp1:rp2;
RealPartonPtr secondRP = (rp1->oY() < rp2->oY() ) ? rp2:rp1;
//let the higher oY recoil to the lower if they are connected. pT only, no plus transfer.
if ( secondRP->mothers.first == firstRP ) {
// cout << firstRP->oY() << " is first mother of " << secondRP << endl;
secondRP->pT += secondRP->theParton->recoil(firstRP->theParton);
secondRP->updateYMinus();
InvEnergy range = sqrt(min(secondRP->theParton->dist2(*firstRP->theParton),
firstRP->theParton->dist2(*secondRP->mothers.second->theParton)/4.0));
secondRP->doEffectiveRecoil(firstRP, range, true, ZERO,
secondRP->theParton->recoil(firstRP->theParton));
}
else if ( secondRP->mothers.second == firstRP ) {
secondRP->pT += secondRP->theParton->recoil(firstRP->theParton);
secondRP->updateYMinus();
InvEnergy range = sqrt(min(secondRP->theParton->dist2(*firstRP->theParton),
firstRP->theParton->dist2(*secondRP->mothers.first->theParton)/4.0));
secondRP->doEffectiveRecoil(firstRP, range, false, ZERO,
secondRP->theParton->recoil(firstRP->theParton));
}
makeCollinear(flucts[i1]);
makeCollinear(flucts[i2]);
// diagnosis(true);
//given yminus should be ok I hope?
}
bool RealPartonState::isMotherStepUp(tRealPartonPtr mom, bool firstSide, InvEnergy scale) {
if ( mom->theParton->valence() ) return false;
tRealPartonPtr grandMom = (firstSide ? mom->mothers.first:mom->mothers.second);
if ( mom->nMothers == 1 ) grandMom = mom->mother;
if ( firstSide && (!mom->children.first.empty() || !mom->firstInt) ) return false;
if ( !firstSide && (!mom->children.second.empty() || !mom->secondInt) ) return false;
if ( mom->theParton->dist2(*grandMom->theParton) > sqr(scale) ) return false;
else return true;
}
bool RealPartonState::inFSRRegion(tRealPartonPtr rp) {
//valence are always ok
if ( rp->theParton->valence() ) return false;
//if parton is outside, but not to an interaction
if ( isOutside(rp, true) ) {
RealPartonPtr mom = (rp->nMothers == 1 ? rp->mother:rp->mothers.first);
//merged partons will share the momentum, and thus only have 1/N of the pt.
int nMerged = 1;
if ( rp->fluct != -1 ) nMerged = flucts[rp->fluct].size();
//check ordering with outside partons without effective partons
//dont chek if already merged, and account for that plus and minus is shared
//between merged partons
if ( !(mom->fluct != -1 && mom->fluct == rp->fluct) ) {
if ( Debug::level > 5 ) cout << "checking FSR PS for " << rp->oY() << ", mom: "
<< mom->oY() << endl;
int nMomMerged = 1;
if ( mom->fluct != -1 ) nMomMerged = flucts[mom->fluct].size();
if ( mom->minus*double(nMomMerged) > rp->minus*double(nMerged) ) return true;
}
//dont bother to check child if rp is interacting since
//we dont know where the other parton is yet
if ( !rp->secondInt && !rp->children.second.empty() ) {
RealPartonPtr child = *rp->children.second.rbegin();
if ( !(child->fluct != -1 && child->fluct == rp->fluct) ) {
if ( Debug::level > 5 ) cout << "checking FSR PS for " << rp->oY()
<< ", child: " << child->oY() << endl;
int nChildMerged = 1;
if ( child->fluct != -1 ) nChildMerged = flucts[child->fluct].size();
if ( child->plus*double(nChildMerged) > rp->plus*double(nMerged) ) return true;
}
}
}
//and check for outside on the other side in the same way.
if ( isOutside(rp, false) ) {
RealPartonPtr mom = (rp->nMothers == 1 ? rp->mother:rp->mothers.second);
//merged partons will share the momentum, and thus only have 1/N of the pt.
int nMerged = 1;
if ( rp->fluct != -1 ) nMerged = flucts[rp->fluct].size();
//check ordering with outside partons without effective partons
//dont chek if already merged, and account for that plus and minus is shared
//between merged partons
if ( !(mom->fluct != -1 && mom->fluct == rp->fluct) ) {
if ( Debug::level > 5 ) cout << "checking FSR PS for " << rp->oY() << ", mom: "
<< mom->oY() << endl;
int nMomMerged = 1;
if ( mom->fluct != -1 ) nMomMerged = flucts[mom->fluct].size();
if ( mom->minus*double(nMomMerged) > rp->minus*double(nMerged) ) return true;
}
//dont bother to check child if rp is interacting since
//we dont know where the other parton is yet
if ( !rp->firstInt && !rp->children.first.empty() ) {
RealPartonPtr child = *rp->children.first.rbegin();
if ( !(child->fluct != -1 && child->fluct == rp->fluct) ) {
if ( Debug::level > 5 ) cout << "checking FSR PS for " << rp->oY()
<< ", child: " << child->oY() << endl;
int nChildMerged = 1;
if ( child->fluct != -1 ) nChildMerged = flucts[child->fluct].size();
if ( child->plus*double(nChildMerged) > rp->plus*double(nMerged) ) return true;
}
}
}
return false;
}
void RealPartonState::fixFSRDoubleCounting(tRealPartonPtr rp) {
RealPartonPtr mom;
if ( rp->nMothers == 1 ) mom = rp->mother;
else if ( isOutside(rp, true) ) mom = rp->mothers.first;
else mom = rp->mothers.second;
if ( rp->fluct != -1 && rp->fluct == mom->fluct ) return;
if ( Debug::level > 5 ) cout << rp->oY() << " is FSR Double Counted with respect to "
<< mom->oY() << ". Merge." << endl;
if ( !rp->recoils.empty() ) {
TransverseMomentum pT;
if ( rp->nMothers == 1 ) pT = rp->theParton->recoil(mom->theParton);
else pT = rp->theParton->recoil(rp->mothers.first->theParton) +
rp->theParton->recoil(rp->mothers.second->theParton);
if ( Debug::level > 5 ) cout << "removing pT (" << pT.x()/GeV << ", " << pT.y()/GeV << ")" << endl;
rp->pT -= pT;
rp->plus -= pT.pt()*exp(-rp->oY());
rp->undoRecoils();
if ( rp->plus < ZERO ) fixNegativePlus(rp);
}
merge(rp, mom);
if ( Debug::level > 5 ) cout << "fixed FSR double counting." << endl;
}
Energy RealPartonState::fixNegativePlus(RealPartonPtr rp) {
//if already nonnegative, nothing has to be done.
if ( rp->plus >= ZERO ) return rp->plus;
//valence partons cant be helped, since there are no parents.
if ( rp->theParton->valence() ) {
cerr << "valence parton with negative plus in RealPartonState::fixNegativePlus" << endl;
return rp->plus;
}
if ( rp->nMothers == 1 ) {
InvEnergy range = sqrt(rp->theParton->dist2(*rp->mother->theParton));
Energy totalPlus = rp->mother->effectivePlusMinus(range, true).first + rp->plus;
//check that there is enough plus in parents, otherwise return false.
if ( totalPlus < ZERO ) {
Throw<RealPartonKinematicException>()
<< "not enough plus in single parent to balance negative plus in "
<< "RealPartonState::fixNegativePlus" << Exception::warning;
return rp->plus;
}
//if enough, divide up so that the parton takes 10% of the effective parents plus.
rp->doEffectiveRecoil(rp->mother, range, true,
-rp->plus + totalPlus/10.0, TransverseMomentum());
rp->plus = totalPlus/10.0;
return totalPlus/10.0;
}
else if ( rp->nMothers == 2 ) {
RealPartonPtr mom1 = rp->mothers.first;
RealPartonPtr mom2 = rp->mothers.second;
InvEnergy range1 = sqrt(min(rp->theParton->dist2(*mom1->theParton),
mom2->theParton->dist2(*mom1->theParton)/4.0));
InvEnergy range2 = sqrt(min(rp->theParton->dist2(*mom2->theParton),
mom1->theParton->dist2(*mom2->theParton)/4.0));
Energy plus1 = mom1->effectivePlusMinus(range1, true).first;
Energy plus2 = mom2->effectivePlusMinus(range2, false).first;
Energy parentPlus = plus1 + plus2;
Energy totalPlus = parentPlus + rp->plus;
//check that there is enough plus in parents, otherwise return false.
if ( parentPlus + rp->plus < ZERO ) {
Throw<RealPartonKinematicException>()
<< "not enough plus in parents to balance negative plus in "
<< "RealPartonState::fixNegativePlus" << Exception::warning;
return rp->plus;
}
//if enough, give enough to put plus to 0
rp->doEffectiveRecoil(mom1, range1, true,
-rp->plus*plus1/parentPlus + totalPlus*plus1/parentPlus/10.0,
TransverseMomentum());
rp->doEffectiveRecoil(mom2, range2, false,
-rp->plus*plus2/parentPlus + totalPlus*plus2/parentPlus/10.0,
TransverseMomentum());
rp->plus = totalPlus/10.0;
return totalPlus/10.0;
}
return rp->plus;
}
bool RealPartonState::isOutside(tRealPartonPtr rp, bool firstSide) {
//to be outside, there can be neither kids not an interaction on one side
//the other side has to have either an interaction or a child. also valence dont count.
if ( (firstSide || rp->nMothers == 1) && rp->children.first.empty() && !rp->firstInt &&
( rp->secondInt || !rp->children.second.empty() ) && !rp->theParton->valence() )
return true;
if ( (!firstSide || rp->nMothers == 1) && rp->children.second.empty() && !rp->secondInt &&
( rp->firstInt || !rp->children.first.empty() ) && !rp->theParton->valence() )
return true;
return false;
}
InvEnergy RealPartonState::childScale(tRealPartonPtr rp, bool firstSide) {
InvEnergy ret = ZERO;
if ( !rp->interacting ) { //if no interaction, then take last child
if ( firstSide ? rp->children.second.empty():rp->children.first.empty() ) return ZERO;
RealPartonPtr child = (firstSide ? *rp->children.second.rbegin():
*rp->children.first.rbegin());
if ( Debug::level > 5 ) cout << " child is " << child->oY() << endl;
ret = sqrt(rp->theParton->dist2(*child->theParton));
if ( Debug::level > 5 ) cout << " child dist is " << ret*GeV << endl;
if (rp->fluct != -1 && child->fluct == rp->fluct) {
ret = childScale(child, firstSide);
if ( Debug::level > 5 ) cout << " in same fluct, recur " << endl;
}
else {
InvEnergy childDist = min(child->intDist, Current<DipoleEventHandler>()->coherenceRange());
if ( child->interacting && child->resolutionScales.lower_bound(rp)->first == rp &&
child->resolutionScales[rp] <
sqr(childDist)*Current<DipoleEventHandler>()->alphaS(childDist) ) {
ret = child->intDist;
if ( Debug::level > 5 ) cout << " child is merged int, new scale: " << ret*GeV << endl;
}
}
}
else {
ret = rp->intDist;
if ( Debug::level > 5 ) cout << " intdist is " << ret*GeV << endl;
}
return ret;
}
bool RealPartonState::checkFixSingleDGLAP(RealPartonPtr rp, InvEnergy scale) {
if ( rp->DGLAPchecked ) return true;
rp->DGLAPchecked = true;
tRealPartonPtr mom = rp->mother;
if ( Debug::level > 5 ) cout << rp->oY() << " is outside (single). mother: " << mom->oY() << endl;
//if the parton is already tagged to merge with the mother, then dont do anything
if ( rp->fluct != -1 && mom->fluct == rp->fluct ) {
if ( Debug::level > 5 ) cout << " already merged, nm." << endl;
return false;
}
if ( isOutside(rp, true) ) {
//if no scale (distance to mother and child) was supplied, then calculate them
if ( scale == ZERO ) scale = childScale(rp, true);
// *** ATTENTION *** not used: InvEnergy motherScale = sqrt(rp->theParton->dist2(*(mom->theParton)));
//check if mom is a further step up in pt, then check her first, with this scale
if ( isMotherStepUp(mom, true, scale) || isMotherStepUp(mom, false, scale) ) {
//if the higher backwards pt max is safe, then this is also ok.
if ( checkFixDGLAP(mom, scale) ) {
if ( Debug::level > 5 ) cout << mom->oY() << " is safe. then so is " << rp->oY() << endl;
return true;
}
}
//compare to mother and fix if not safe
if ( !rp->DGLAPSafe(scale) ) {
fixDGLAP(rp);
return false;
}
}
if ( isOutside(rp, false) ) {
//if no scale (distance to mother and child) was supplied, then calculate them
if ( scale == ZERO ) scale = childScale(rp, false);
// *** ATTENTION *** not used: InvEnergy motherScale = sqrt(rp->theParton->dist2(*(mom->theParton)));
//check if mom is a further step up in pt, then check her first, with this scale
if ( isMotherStepUp(mom, true, scale) || isMotherStepUp(mom, false, scale) ) {
//if the higher backwards pt max is safe, then this is also ok.
if ( checkFixDGLAP(mom, scale) ) {
if ( Debug::level > 5 ) cout << mom->oY() << " is safe. then so is " << rp->oY() << endl;
return true;
}
}
//compare to mother and fix if not safe
if ( !rp->DGLAPSafe(scale) ) {
fixDGLAP(rp);
return false;
}
}
return true;
}
void RealPartonState::fixSingleDGLAP(RealPartonPtr rp) {
if ( rp->fluct != -1 && rp->fluct == rp->mother->fluct )
return; //already merged with the mother.
else if ( rp->fluct != -1 && rp->recoils.empty() ) {
if ( !rp->interacting )
merge(rp, rp->mother); //no recoils that has to be undone.
return;
}
TransverseMomentum pT = rp->theParton->recoil(rp->mother->theParton);
rp->pT -= pT; //remove rec from itself
Energy originalPlus = (pT).pt()*exp(-rp->oY());
rp->plus -= originalPlus;
rp->undoRecoils();
//now readd some recoil limited by the DGLAP scale.
InvEnergy r = ZERO;
if ( rp->interacting ) r = rp->intDist;
else if ( !rp->children.first.empty() )
r = sqrt((*rp->children.first.rbegin())->theParton->dist2(*rp->theParton));
else if ( !rp->children.second.empty() )
r = sqrt((*rp->children.second.rbegin())->theParton->dist2(*rp->theParton));
else cerr << "fixing single DGLAP, but no children, or interaction scales!" << endl;
//find from the resolutionscale (alpha_s(r)*r^2) what the minimum r couldve been.
//since alpha_s(r) has no analytic inverse, do 3 recursive steps from cohernce range.
if ( r > Current<DipoleEventHandler>()->coherenceRange() )
r = Current<DipoleEventHandler>()->coherenceRange();
InvEnergy resScale = sqrt(rp->resolutionScales[rp->mother]/
Current<DipoleEventHandler>()->alphaS(r));
resScale = sqrt(rp->resolutionScales[rp->mother]/
Current<DipoleEventHandler>()->alphaS(resScale));
resScale = sqrt(rp->resolutionScales[rp->mother]/
Current<DipoleEventHandler>()->alphaS(resScale));
double penalty = resScale/r;
if ( !rp->interacting || true )
penalty = 0.0;
else if ( Debug::level > 5 ) {
cout << " SINGLE: penalty for " << rp->oY() << " is " << penalty << ", PT1 to " << rp->mother->oY()
<< " goes from " << pT.pt()/GeV << " to "
<< penalty*pT.pt()/GeV << endl;
}
Energy plusRec = (penalty*pT).pt()*exp(-rp->oY());
rp->pT += penalty*pT;
rp->plus += plusRec;
rp->doEffectiveRecoil(rp->mother, sqrt(rp->theParton->dist2(*rp->mother->theParton)),
true, plusRec, penalty*pT);
if ( !rp->interacting || true )
merge(rp, rp->mother);
else if ( rp->plus < ZERO || rp->mother->plus < ZERO ) {
merge(rp, rp->mother);
}
if ( (rp->plus < ZERO || rp->mother->plus < ZERO) )
cerr << "negative p+ after fixsingleDGLAP! :o" << endl;
}
bool RealPartonState::checkFixDGLAP(RealPartonPtr rp, InvEnergy scale) {
if ( rp->DGLAPchecked ) return true;
rp->DGLAPchecked = true;
//DGLAP should be checked only for partons on the "outside" of a chain.
//if it is
if ( isOutside(rp, true) ) { //if kids/interaction on only second side
tRealPartonPtr mom;
if ( rp->nMothers == 1 ) mom = rp->mother;
else mom = rp->mothers.first;
if ( Debug::level > 5 ) cout << rp->oY() << " is outside 1. mother: " << mom->oY() << endl;
//if the parton is already tagged to merge with the mother, then dont do anything
if ( rp->fluct != -1 && mom->fluct == rp->fluct ) {
if ( Debug::level > 5 ) cout << rp->oY() << " already merged, nm." << endl;
return false;
}
//if no scale (distance to mother and child) was supplied, then calculate them
if ( scale == ZERO ) scale = childScale(rp, true);
InvEnergy motherScale = sqrt(rp->theParton->dist2(*(mom->theParton)));
if ( Debug::level > 5 ) cout << rp->oY() << " scale: " << scale*GeV
<< ", motherscale: " << motherScale*GeV << endl;
//if grandmother is a further step up in pt, then go there and check first
//then compare to THIS childs scale, to get full suppression
if ( isMotherStepUp(mom, true, motherScale) || isMotherStepUp(mom, false, motherScale) ) {
if ( Debug::level > 5 ) cout << " mother is step up. recur to " << mom->oY() << endl;
//if the higher backwards pt max is safe, then this is also ok.
if ( checkFixDGLAP(mom, scale) ) {
if ( Debug::level > 5 ) cout << mom->oY() << " is safe. then so is " << rp->oY() << endl;
return true;
}
//if the backwards max is removed, then this is the new max, and should be checked
}
if ( Debug::level > 5 ) cout << " checkDGLAPsafe." << endl;
//check if the recoil should be removed or not
if ( !rp->checkDGLAPSafe(mom, rp, scale) ) {
fixDGLAP(rp);
return false;
}
}
//the same thing, if kids/interaction on only first side
if ( isOutside(rp, false) ) {
tRealPartonPtr mom;
if ( rp->nMothers == 1 ) mom = rp->mother;
else mom = rp->mothers.second;
if ( Debug::level > 5 ) cout << rp->oY() << " is outside 2. mother: " << mom->oY() << endl;
if ( rp->fluct != -1 && mom->fluct == rp->fluct ) {
if ( Debug::level > 5 ) cout << rp->oY() << " already merged, nm." << endl;
return false;
}
if ( scale == ZERO ) scale = childScale(rp, false);
InvEnergy motherScale = sqrt(rp->theParton->dist2(*(mom->theParton)));
if ( Debug::level > 5 ) cout << rp->oY() << " scale: " << scale*GeV
<< ", motherscale: " << motherScale*GeV << endl;
if ( isMotherStepUp(mom, true, motherScale) || isMotherStepUp(mom, false, motherScale) ) {
if ( Debug::level > 5 ) cout << " mother is step up. recur to " << mom->oY() << endl;
if ( checkFixDGLAP(mom, scale) ) {
if ( Debug::level > 5 ) cout << mom->oY() << " is safe. then so is " << rp->oY() << endl;
return true;
}
}
if ( Debug::level > 5 ) cout << " checkDGLAPsafe." << endl;
if ( !rp->checkDGLAPSafe(mom, rp, scale) ) {
fixDGLAP(rp);
return false;
}
}
if ( Debug::level > 5 ) cout << rp->oY() << " is not outside" << endl;
return true;
}
void RealPartonState::fixDGLAP(RealPartonPtr rp) {
if ( Debug::level > 5 && monitored ) cout << rp->oY() << " is unordered. Remove." << endl;
if ( rp->nMothers == 1 ) {
fixSingleDGLAP(rp);
return;
}
RealPartonPtr mother = ((rp->children.first.empty() && !rp->firstInt) ? rp->mothers.first:rp->mothers.second);
if ( rp->fluct != -1 && flucts[rp->fluct].find(mother) != flucts[rp->fluct].end() ) {
return; //already merged with the mother.
}
else if ( rp->fluct != -1 && rp->recoils.empty() ) {
if ( !rp->interacting )
merge(rp, mother); //recoils already undone, dont reundo them.
return;
}
if ( (rp->plus < ZERO || mother->plus < ZERO) && Debug::level > 5 )
cerr << "negative p+ before fixdglap" << endl;
//remove recoil between itself and other mother, but save recs from future emissions.
//recoil with merging mother doesnt matter, since averaged anyways.
RealPartonPtr otherMother = (mother == rp->mothers.first) ?
rp->mothers.second:rp->mothers.first;
TransverseMomentum pT1 = rp->theParton->recoil(mother->theParton);
TransverseMomentum pT2 = rp->theParton->recoil(otherMother->theParton);
rp->pT -= pT1 + pT2; //remove rec from itself
rp->plus -= (pT1 + pT2).pt()*exp(-rp->oY());
rp->undoRecoils();
//now readd some recoil limited by the DGLAP scale.
InvEnergy r2;
if ( rp->interacting ) r2 = rp->intDist;
else {
if ( mother == rp->mothers.first )
r2 = sqrt((*rp->children.second.rbegin())->theParton->dist2(*rp->theParton));
else
r2 = sqrt((*rp->children.first.rbegin())->theParton->dist2(*rp->theParton));
}
if ( r2 > Current<DipoleEventHandler>()->coherenceRange() )
r2 = Current<DipoleEventHandler>()->coherenceRange();
InvEnergy resScale = sqrt(rp->resolutionScales[mother]/
Current<DipoleEventHandler>()->alphaS(r2));
resScale = sqrt(rp->resolutionScales[mother]/
Current<DipoleEventHandler>()->alphaS(resScale));
resScale = sqrt(rp->resolutionScales[mother]/
Current<DipoleEventHandler>()->alphaS(resScale));
double penalty = resScale/r2;
if ( !rp->interacting )
penalty = 0.0;
else if ( Debug::level > 5 ) {
cout << " penalty for " << rp->oY() << " is " << penalty << ", PT1 to " << mother->oY()
<< " goes from " << pT1.pt()/GeV << " to "
<< penalty*pT1.pt()/GeV << endl;
}
double P1 = rp->theParton->dist2(*otherMother->theParton)/
(rp->theParton->dist2(*mother->theParton) + rp->theParton->dist2(*otherMother->theParton));
double P2 = 1.0 - P1;
Energy plusRec = (penalty*pT1 + pT2).pt()*exp(-rp->oY());
rp->pT += penalty*pT1 + pT2;
rp->plus += plusRec;
if ( Debug::level > 5 ) {
cout << " readding pT of (" << (penalty*pT1 + pT2).x()/GeV << ", "
<< (penalty*pT1 + pT2).y()/GeV << ")" << endl;
cout << " pT is (" << rp->pT.x()/GeV << ", " << rp->pT.y()/GeV << endl;
cout << "second plus recoil is " << P2*plusRec/GeV << endl;
}
rp->doEffectiveRecoil(mother, sqrt(min(rp->theParton->dist2(*mother->theParton), mother->theParton->dist2(*otherMother->theParton)/4.0)), mother == rp->mothers.first, P1*plusRec, penalty*pT1);
rp->doEffectiveRecoil(otherMother, sqrt(min(rp->theParton->dist2(*otherMother->theParton), mother->theParton->dist2(*otherMother->theParton)/4.0)), otherMother == rp->mothers.first, P2*plusRec, pT2);
if ( rp->plus + mother->plus < ZERO ) { //handle this case. should be kindof rare.
// if ( Debug::level > 5 ) cout << "returned too much plus in fixDGLAP. taking some back again... :/" << endl;
// if ( Debug::level > 5 ) cout << " plus is " << rp->plus/GeV << endl;
rp->doEffectiveRecoil(otherMother, sqrt(min(rp->theParton->dist2(*otherMother->theParton), mother->theParton->dist2(*otherMother->theParton)/4.0)), otherMother == rp->mothers.first, -rp->plus, TransverseMomentum());
rp->plus = ZERO;
}
if ( !rp->interacting || true )
merge(rp, mother);
else if ( rp->plus < ZERO || mother->plus < ZERO ) {
merge(rp, mother);
if ( (rp->plus < ZERO || mother->plus < ZERO) )
Throw<RealPartonKinematicException>()
<< "even after refix negative p+ after fixdglap with interacting!" << Exception::warning;
}
if ( otherMother->plus < ZERO ) {
merge(rp, mother);
merge(rp, otherMother);
}
if ( otherMother->plus < ZERO )
Throw<RealPartonKinematicException>()
<< "other mother in RealPartonState::fixDGLAP() has negative p+!" << Exception::warning;
}
void RealPartonState::fixUnOrdered(RealPartonPtr rp, bool forced) {
if ( !rp ) {
cerr << "fixUnOrdered called without a realparton!" << endl;
}
if ( rp->theParton->valence() ) return;
RealPartonPtr closestMother;
if ( rp->nMothers == 1 ) closestMother = rp->mother;
else
closestMother = ( (rp->theParton->dist2(*rp->mothers.first->theParton)
< rp->theParton->dist2(*rp->mothers.second->theParton))
? rp->mothers.first:rp->mothers.second);
//dont merge interacting unless forced
if ( rp->interacting && !forced && rp->nMothers == 2 ) {
//force the unordered parents to merge
if ( !rp->isFirstOrdered() && !rp->mothers.first->theParton->valence() )
fixUnOrdered(rp->mothers.first, true);
if ( !rp->isSecondOrdered() && !rp->mothers.second->theParton->valence() )
fixUnOrdered(rp->mothers.second, true);
return;
}
else if ( rp->interacting && !forced && rp->nMothers == 1 ) {
if ( !rp->isSingleOrdered() && !rp->mother->theParton->valence() )
fixUnOrdered(rp->mother, true);
return;
}
if ( rp->fluct != -1 && rp->fluct == closestMother->fluct ) {
if ( Debug::level > 5 && closestMother->theParton->valence() )
cout << "WARNING, recuring back to valence in RPS::fixUnOrdered" << endl;
fixUnOrdered(closestMother, forced);
return;
}
//remove recoil between itself and other mother, but save recs from future emissions.
//recoil with merging mother doesnt matter, since averaged anyways.
RealPartonPtr otherMother = (closestMother == rp->mothers.first) ?
rp->mothers.second:rp->mothers.first;
TransverseMomentum pT1 = rp->theParton->recoil(closestMother->theParton);
TransverseMomentum pT2;
if ( rp->nMothers == 1 ) pT2 = TransverseMomentum();
else pT2 = rp->theParton->recoil(otherMother->theParton);
rp->pT -= pT1 + pT2; //remove rec from itself
rp->plus -= (pT1 + pT2).pt()*exp(-rp->oY());
rp->undoRecoils();
merge(rp, closestMother);
}
bool RealPartonState::checkPlusMinus() {
Energy totalPlus = ZERO;
Energy totalMinus = ZERO;
Energy onShellMinus = ZERO;
for ( map<tPartonPtr,RealPartonPtr>::iterator it = partons.begin();
it != partons.end(); it++ ) {
tRealPartonPtr rp = (*it).second;
if ( rp->keep == RealParton::NO ) continue;
totalPlus += rp->plus;
totalMinus += rp->minus;
onShellMinus += rp->givenMinus;
if ( rp->theParton->valence() )
onShellMinus += sqr(rp->theParton->valencePT().pt())/rp->theParton->valencePlus();
}
if ( Debug::level > 5 ) cout << " total plus: " << totalPlus/GeV << endl
<< " total minus: " << totalMinus/GeV << endl
<< "onshell minus: " << onShellMinus/GeV << endl;
return true;
}
bool RealPartonState::checkForNegatives() {
bool ret = false;
for ( map<tPartonPtr,RealPartonPtr>::iterator it = partons.begin();
it != partons.end(); it++ ) {
tRealPartonPtr rp = (*it).second;
if ( rp->keep == RealParton::NO ) continue;
if ( rp->minus < ZERO ) {
// Throw<RealPartonKinematicException>()
// << "negative minus in real parton at oY " << rp->oY() << Exception::warning;
ret = true;
}
if ( rp->plus < ZERO ) {
// Throw<RealPartonKinematicException>()
// << "negative plus in real parton at oY " << rp->oY() << Exception::warning;
ret = true;
}
}
return ret;
}
bool RealPartonState::checkFlucts() {
bool ok = true;
for ( vector< RealPartonSet>::iterator it = flucts.begin(); it != flucts.end(); it++ ) {
RealPartonSet & fluct = *it;
if ( fluct.empty() ) continue;
Energy plus = (*fluct.begin())->plus;
Energy minus = (*fluct.begin())->minus;
Energy ptx = (*fluct.begin())->pT.x();
Energy pty = (*fluct.begin())->pT.y();
double oy = (*fluct.begin())->oY();
for( RealPartonSet::iterator jt = ++fluct.begin(); jt != fluct.end(); jt++) {
RealPartonPtr rp = *jt;
if ( abs((rp->plus - plus)/plus) > 0.01 ) {
ok = false;
if ( Debug::level > 5 )
cout << "real parton at " << rp->oY() << " has plus " << rp->plus/GeV << ", while "
<< oy << " in the same fluct has " << plus/GeV << endl;
}
if ( abs((rp->minus - minus)/minus) > 0.01 ) {
ok = false;
if ( Debug::level > 5 )
cout << "real parton at " << rp->oY() << " has minus " << rp->minus/GeV << ", while "
<< oy << " in the same fluct has " << minus/GeV << endl;
}
if ( abs((rp->pT.x() - ptx)/ptx) > 0.01 ) {
ok = false;
if ( Debug::level > 5 )
cout << "real parton at " << rp->oY() << " has ptx " << rp->pT.x()/GeV << ", while "
<< oy << " in the same fluct has " << ptx/GeV << endl;
}
if ( abs((rp->pT.y() - pty)/pty) > 0.01 ) {
ok = false;
if ( Debug::level > 5 )
cout << "real parton at " << rp->oY() << " has pty " << rp->pT.y()/GeV << ", while "
<< oy << " in the same fluct has " << pty/GeV << endl;
}
}
}
return ok;
}
void RealPartonState::mergeVirtuals() {
if ( !checkFlucts() ) {
Throw<RealPartonKinematicException>()
<< "the fluctuations are not collinear in mergevirtuals!!" << Exception::warning;
}
for ( vector< RealPartonSet>::iterator it = flucts.begin(); it != flucts.end(); it++ ) {
RealPartonSet & fluct = *it;
if ( fluct.empty() ) continue;
RealPartonPtr mother = *fluct.begin();
if ( Debug::level > 5 ) cout << "merging a fluct from " << mother->oY() << " with "
<< fluct.size() << " partons" << endl;
for( RealPartonSet::iterator jt = ++fluct.begin(); jt != fluct.end(); jt++) {
RealPartonPtr rp = *jt;
if ( rp->theParton->valence() ) {
if ( Debug::level > 5 ) cout << " dont merge valence partons!" << endl;
continue;
}
if ( Debug::level > 5 ) cout << " removing " << rp->oY() << endl;
if ( abs(mother->y - rp->y) > 0.00001 ) {
cout << " different y!!!! " << endl;
plotState(true);
}
mother->plus += rp->plus;
mother->pT += rp->pT;
rp->theParton->onShell(false);
rp->keep = RealParton::NO;
for ( RealPartonSet::iterator kt = rp->children.first.begin();
kt != rp->children.first.end(); kt++ ) {
if ( (*kt)->nMothers == 2 ) {
(*kt)->mothers.second = mother;
mother->children.first.insert(*kt);
}
if ( (*kt)->nMothers == 1 ) {
(*kt)->mother = mother;
mother->children.first.insert(*kt);
}
}
for ( RealPartonSet::iterator kt = rp->children.second.begin();
kt != rp->children.second.end(); kt++ ) {
if ( (*kt)->nMothers == 2 ) {
(*kt)->mothers.first = mother;
mother->children.second.insert(*kt);
}
if ( (*kt)->nMothers == 1 ) {
(*kt)->mother = mother;
mother->children.second.insert(*kt);
}
}
if ( !rp->theParton->valence() ) {
if ( rp->nMothers == 2 ) {
rp->mothers.first->children.second.erase(rp);
rp->mothers.second->children.first.erase(rp);
rp->mothers.first = RealPartonPtr();
rp->mothers.second = RealPartonPtr();
}
if ( rp->nMothers == 1 ) {
rp->mother->children.second.erase(rp);
rp->mother->children.first.erase(rp);
rp->mother = RealPartonPtr();
}
}
}
mother->updateYMinus();
// plotState(true);
}
// cout << "done merging\n";
}
void RealPartonState::doEvolution() {
// cout << "entering doEvolution" << endl;
for ( RealPartonSet::iterator it = toCheck.begin(); it != toCheck.end(); it++ ) {
RealPartonPtr rp = *it;
currentY = rp->theParton->oY();
if ( !rp->quickSetYES() ) {
rp->checkEmissionProblem();
break;
}
rp->checkEmissionProblem();
}
}
bool RealPartonState::findConsistentEvolution(RealPartonSet::iterator it) {
if ( it == toCheck.end() ) {
return true;}
partonCalls++;
if ( partonCalls > 1000000 ) {return false;}
tRealPartonPtr rp = *it;
currentY = rp->theParton->oY();
if ( !(rp->setYES()) ) {
unOrdered++;
if ( monitored && interactions.size() > 1 )
cout << "unordered" << endl;
// cout << rp->oY() << " unordered" << endl;
return rp->setNO();
}
if ( !rp->theParton->valence() && (!rp->mothers.first || !rp->mothers.second) ) {
if ( monitored && interactions.size() > 1 )
cout << "non-valence without mothers, NO" << endl;
// cout << rp->oY() << " mother structure error" << endl;
return rp->setNO();
}
if ( rp->mothers.first && rp->mothers.first == rp->mothers.second ) {
if ( monitored && interactions.size() > 1 )
cout << "samme mothers" << endl;
// cout << rp->oY() << " mother structure error" << endl;
return rp->setNO();
}
// cout << rp->oY() << " passed first tests. check." << endl;
bool existHistory = false;
while ( !existHistory ) {
existHistory = findConsistentEvolution(++it);
currentY = rp->theParton->oY();
if ( existHistory && !(rp->DGLAPSafe()) ) {
if ( monitored && interactions.size() > 1 )
cout << "unDGLAP" << endl;
// cout << rp->oY() << ", (" << rp->theParton->position().x()*GeV << ", "
// << rp->theParton->position().y()*GeV << ") is unDGLAP" << endl;
unDGLAP++;
rp->unDGLAP++;
existHistory = false;
}
if ( existHistory && !(rp->interacting) && !(rp->theParton->valence()) && !(rp->hasChild()) )
existHistory = false;
// if ( !existHistory && it != toCheck.end() ) {
// cout << rp->oY() << " found no history at " << (*it)->oY() << ", set future NO" << endl;
// // plotState(true);
// }
if ( !existHistory ) setFutureNO(rp); //behovs denna verkligen??
if ( !existHistory && (it == toCheck.end() || (*it)->interacting) ) return rp->setNO();
if ( !existHistory && (*it)->theParton->valence() ) return rp->setNO();
}
return existHistory;
}
bool RealPartonState::controlEvolution(tDipolePtr intDip, tDipolePtr otherIntDip) {
if ( failedInts.find(intDip) != failedInts.end() ) return false;
//think through how to do this without redoing the evolution!
if ( intDip->interacted() ) {
rescatter = true;
// return false; //removes rescattering
}
else rescatter = false;
unOrdered = 0;
unDGLAP = 0;
partonCalls = 0;
newInteraction(intDip, otherIntDip, true, true);
// checkOnlyNew(intDip);
checkHistory(getReal(intDip->partons().first));
checkHistory(getReal(intDip->partons().second));
checkFuture();
cout << "new evo------------------- " << toCheck.size() << " to check ---------------------" << endl;
cout << "-- interaction partons are " << getReal(intDip->partons().first)->oY() << " and "
<< getReal(intDip->partons().second)->oY() << endl;
// checkInteraction(intDip);
// cout << "starting real evo" << endl;
// diagnosis(true);
if ( (getReal(intDip->partons().first)->interactions.size() > 1 &&
getReal(intDip->partons().first)->keep == RealParton::NO) ||
(getReal(intDip->partons().second)->interactions.size() > 1 &&
getReal(intDip->partons().second)->keep == RealParton::NO) ) {
cout << "interacting valence NO" << endl;
revertToPrevious(intDip);
failedInts.insert(intDip);
return false;
}
if ( !toCheck.empty() && !rescatter ) { //think through rescatter. will intpT still be ok?
for ( RealPartonSet::iterator it = toCheck.begin(); !findConsistentEvolution(it); it++ ) {
if ( (*it)->interacting || (*it)->theParton->valence() ) {
cout << "failed evo, unordered: " << unOrdered << ", unDGLAP: " << unDGLAP << endl;
revertToPrevious(intDip);
failedInts.insert(intDip);
return false;
}
}
}
if ( !toCheck.empty() )
currentY = (**(--toCheck.end())).theParton->oY();
cout << "found evo" << endl;
// for ( RealPartonSet::iterator it = toCheck.begin(); it != toCheck.end(); it++ )
// (*it)->checkInteractionRecoils();
return true;
}
bool RealPartonState::singleControlEvolution(tDipolePtr intDip, tDipolePtr otherIntDip,
bool firstInt, bool secondInt,
Energy rec1, Energy rec2) {
ostream & log = CurrentGenerator::log();
static DebugItem printstepdetails("DIPSY::PrintStepDetails", 6);
static DebugItem checkkinematics("DIPSY::CheckKinematics", 6);
if ( failedInts.find(intDip) != failedInts.end() ) return false;
if ( intDip->interacted() ) {
if ( monitored )
cout << " rescatter, autopass evolution" << endl;
rescatter = true;
}
else rescatter = false;
if ( printstepdetails )
log << "<DIPSY> enter singlecontrolevolution." << endl;
unOrdered = 0;
unDGLAP = 0;
partonCalls = 0;
newInteraction(intDip, otherIntDip, firstInt, secondInt, rec1, rec2);
reset();
checkAllInteracting();
for ( RealPartonSet::iterator it = toCheck.begin(); it != toCheck.end(); it++ )
if ( !(*it)->setYES() ) fixUnOrdered(*it);
if ( printstepdetails ) log << "<DIPSY> done with forward sweep." << endl;
if ( monitored ) plotState(true);
for ( RealPartonSet::reverse_iterator it = toCheck.rbegin(); it != toCheck.rend(); it++ )
checkFixDGLAP(*it);
if ( printstepdetails )
log << "<DIPSY> done with DGLAP corrections." << endl;
if ( monitored ) plotState(true);
for ( RealPartonSet::iterator it = toCheck.begin(); it != toCheck.end(); it++ ) {
if ( inFSRRegion(*it) ) {
if ( printstepdetails )
log << "<DIPSY> FSR conflict found (check) at " << (*it)->oY() << endl;
fixFSRDoubleCounting(*it);
if ( printstepdetails ) log << " fixed, check." << endl;
}
}
if ( printstepdetails ) log << "<DIPSY> done with FSR matching" << endl;
if ( monitored ) plotState(true);
//make sure all fluctuations are collinear. In general not needed.
for ( int i = 0; i < int(flucts.size()); i++ ) {
makeCollinear(flucts[i]);
}
if ( checkkinematics && checkForNegatives() )
Throw<RealPartonKinematicException>()
<< "negatives found at end of cascade evo!" << Exception::warning;
return true;
}
bool RealPartonState::nonRecursiveControlEvolution(tDipolePtr intDip, tDipolePtr otherIntDip) {
if ( failedInts.find(intDip) != failedInts.end() ) return false;
if ( intDip->interacted() ) {
if ( monitored )
cout << " rescatter, autopass evolution" << endl;
rescatter = true;
}
else rescatter = false;
// cout << "entering nonrecursivecontrol\n";
// diagnosis(false);
unOrdered = 0;
unDGLAP = 0;
partonCalls = 0;
newInteraction(intDip, otherIntDip, true, true);
reset();
checkAllInteracting();
// cout << "enter loop to find evo------------------------------------------" << endl;
while ( true ) {
undoEmissions();
doEvolution();
// cout << "did evo" << endl;
RealPartonPtr worstProblem = findWorstProblem();
// if ( worstProblem )
// cout << "worst problem at " << worstProblem->oY() << ", ("
// << worstProblem->theParton->position().x()*GeV << ", "
// << worstProblem->theParton->position().y()*GeV << ")" << endl;
// else {
// cout << "no problem! evo ok for " << this << "" << endl;
// // diagnosis(false);
// }
if ( !worstProblem ) return true;
if ( !fix(worstProblem) ) {
// cout << "failed fixing theproblem, interaction failed" << endl;
revertToPrevious(intDip);
failedInts.insert(intDip);
return false;
}
}
}
bool RealPartonState::fullControlEvolution(tDipolePtr intDip, tDipolePtr otherIntDip,
bool firstInt, bool secondInt,
Energy rec1, Energy rec2) {
if ( intDip->dipoleState().handler().eventFiller().mode() == 1 )
return controlEvolution(intDip, otherIntDip);
if ( intDip->dipoleState().handler().eventFiller().mode() == 2 )
return singleControlEvolution(intDip, otherIntDip, firstInt, secondInt, rec1, rec2);
if ( intDip->dipoleState().handler().eventFiller().mode() == 3 )
return nonRecursiveControlEvolution(intDip, otherIntDip);
if ( failedInts.find(intDip) != failedInts.end() ) return false;
//think through how to do this without redoing the evolution!
if ( intDip->interacted() ) {
if ( monitored )
cout << " rescatter, autopass evolution" << endl;
rescatter = true;
cout << "rescatter!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n";
// return false; //removes rescattering
}
else rescatter = false;
unOrdered = 0;
unDGLAP = 0;
partonCalls = 0;
newInteraction(intDip, otherIntDip, true, true);
reset();
checkAllInteracting();
cout << "new evo---------------------------------------------------------" << endl;
for ( RealPartonSet::iterator it = toCheck.begin(); !findConsistentEvolution(it); it++ ) {
if ( (*it)->interacting || (*it)->theParton->valence() ) {
// cout << "failed evo, unordered: " << unOrdered << ", unDGLAP: " << unDGLAP << endl;
// for ( RealPartonSet::iterator jt = toCheck.begin(); jt != toCheck.end(); jt++ ) {
// (*jt)->setYES();
// }
// plotState(true);
revertToPrevious(intDip);
failedInts.insert(intDip);
return false;
}
}
cout << " found consistent evolution!" << endl;
// diagnosis(false);
// cout << "checked" << endl;
currentY = (**(--toCheck.end())).theParton->oY();
// cout << " did " << partonCalls << " calls on " << toCheck.size() << " partons.";
return true;
}
void RealPartonState::reset() {
for ( RealPartonSet::iterator it = oValence.begin(); it != oValence.end(); it++ ) {
(*it)->setYES();
}
for ( map<tPartonPtr,RealPartonPtr>::iterator it = partons.begin();
it != partons.end(); it++ ) {
(*it).second->givenMinus = ZERO;
(*it).second->fluct = -1;
(*it).second->DGLAPchecked = false;
if ( (*it).second->keep == RealParton::YES && !(*it).second->theParton->valence() ) {
(*it).second->setNO();
}
(*it).second->intRecoil = TransverseMomentum();
}
// for ( RealPartonSet::iterator it = valence.begin(); it != valence.end(); it++ ) {
// (*it)->setYES();
// }
totalRecoil = TransverseMomentum();
backwardsPartons.clear();
flucts.clear();
suspects.clear();
toCheck.clear();
minusDeficit = (*interactions.begin())->dipoleState().minusDeficit();
}
void RealPartonState::undoEmissions() {
for ( RealPartonSet::reverse_iterator it = toCheck.rbegin();
it != toCheck.rend(); it++ ) {
RealPartonPtr rp = *it;
if ( !rp->theParton->valence() )
rp->setNO();
else {
if ( rp->mothers.first ) rp->mothers.first->children.second.erase(rp);
if ( rp->mothers.second ) rp->mothers.second->children.first.erase(rp);
rp->quickSetYES();
}
}
}
void RealPartonState::checkOnlyNew(tDipolePtr intDip) {
cout << "implement checkonlynew if you really want to use it!!" << endl;
// toCheck.clear();
// //interacting partons cannot be left NO, even if earlier evos did so.
// if ( getReal(intDip->partons().first)->keep == RealParton::NO )
// toCheck.insert(getReal(intDip->partons().first));
// if ( getReal(intDip->partons().second)->keep == RealParton::NO )
// toCheck.insert(getReal(intDip->partons().second));
// for ( map<tPartonPtr,RealPartonPtr>::iterator it = partons.begin();
// it != partons.end(); it++ ) {
// if ( ((*it).second->interactions.size() == 1) &&
// (*it).second->interactions.find(intDip) != (*it).second->interactions.end() ) {
// toCheck.insert(it->second);
// }
// }
// currentY = 0.0;
}
void RealPartonState::checkInteraction(tDipolePtr intDip) {
toCheck.clear();
RealPartonPtr rp1 = getReal(intDip->partons().first);
checkHistory(rp1);
RealPartonPtr rp2 = getReal(intDip->partons().second);
checkHistory(rp2);
}
void RealPartonState::checkHistory(tRealPartonPtr rp) {
toCheck.insert(rp);
if ( rp->oMothers.first )
checkHistory(rp->oMothers.first);
if ( rp->oMothers.second )
checkHistory(rp->oMothers.second);
}
void RealPartonState::checkFuture() {
RealPartonSet future;
for ( RealPartonSet::iterator it = toCheck.begin(); it != toCheck.end(); it++ ) {
RealPartonPtr rp = *it;
if ( !rp->future.empty() )
future.insert(rp->future.begin(), rp->future.end());
}
if ( !future.empty() )
toCheck.insert(future.begin(), future.end());
}
void RealPartonState::checkAllInteracting() { //goes through all partons.
for ( map<tPartonPtr,RealPartonPtr>::iterator it = partons.begin();
it != partons.end(); it++ ) {
if ( !((*it).second->interactions.empty()) ) {
toCheck.insert(it->second);
}
if ( it->first->valence() ) it->second->setValenceMomentum();
}
currentY = 0.0;
}
void RealPartonState::setFutureNO(tRealPartonPtr rp) { //remade to go from high y to low/
for ( RealPartonSet::reverse_iterator it = toCheck.rbegin(); *it != rp; it++) {
if ( (*it)->keep == RealParton::YES && !(*it)->theParton->valence() )
(*it)->setNO();
}
}
void RealPartonState::setNO(const RealPartonSet & rps) {
for ( RealPartonSet::const_reverse_iterator it = rps.rbegin();
it != rps.rend(); it++)
if ( (*it)->keep == RealParton::YES &&
!(*it)->theParton->valence() ) (*it)->setNO();
}
RealPartonPtr RealPartonState::findWorstProblem() {
// cout << "entering worst problem" << endl;
Energy largestScale = ZERO;
RealPartonPtr ret;
for ( RealPartonSet::iterator it = toCheck.begin(); it != toCheck.end(); it++) {
if ( (*it)->keep == RealParton::NO ) break;
Energy scale = (*it)->problemScale();
// cout << "got scale " << scale/GeV << " at " << (*it)->oY() << ", ("
// << (*it)->theParton->position().x()*GeV << ", "
// << (*it)->theParton->position().y()*GeV << ")" << endl;
if ( scale > largestScale ) {
largestScale = scale;
ret = *it;
}
}
return ret;
}
bool RealPartonState::fix(RealPartonPtr problem) {
// cout << "entered fix" << endl;
RealPartonPtr cause = problem->findCause();
if ( !cause ) return false;
cause->setNO();
toCheck.erase(cause);
return true;
}
void RealPartonState::saveState() {
if ( Debug::level > 5 ) cout << "entering RPS::saveState" << endl;
for ( map<tPartonPtr, RealPartonPtr>::iterator it = partons.begin();
it != partons.end(); it++) {
(*it).second->saveState();
}
for ( RealPartonSet::iterator it = valence.begin(); //can be removed?
it != valence.end(); it++) {
(*it)->saveState();
}
failedInts.clear(); //when interaction found, the fails may now be ok.
if ( Debug::level > 5 ) cout << " failedInts cleared" << endl;
cTotalRecoil = totalRecoil;
cFlucts = flucts;
}
void RealPartonState::revertToPrevious(DipolePtr intDip) {
if ( Debug::level > 5 ) cout << "revert to previous called with intdip " << intDip << endl;
// if ( rescatter ) {
// //if rescattering, some extra care has to be taken in
// //removing the interaction tags from the partons,
// //so dont remove any tags the normal way, and call removeLastRescatter instead.
// removeLastRescatter(intDip);
// intDip = DipolePtr();
// }
for ( map<tPartonPtr,RealPartonPtr>::iterator it = partons.begin();
it != partons.end(); it++) {
it->second->revert(intDip);
}
interactions.pop_back();
doesInts.pop_back();
totalRecoil = cTotalRecoil;
flucts = cFlucts;
}
void RealPartonState::removeLastRescatter(tDipolePtr intDip) {
if ( Debug::level > 5 ) cout << "entering removeLastRescatter for intdip " << intDip << endl;
//first check which of the two partons has been interacting in previous interactions
bool firstOldInt = false;
bool secondOldInt = false;
//loop over all but the last interaction, to not include this last rescattering
if ( Debug::level > 5 ) cout << " number of interactions total: " << interactions.size() << endl;
list<pair<bool, bool> > ::iterator jt = doesInts.begin();
for ( list<tDipolePtr>::iterator it = interactions.begin(); it != --interactions.end();
it++, jt++ ) {
if ( Debug::level > 5 ) cout << " old int " << *it << endl;
if ( *it == intDip ) {
if ( Debug::level > 5 ) cout << " found old rescattering" << endl;
if ( jt->first ) firstOldInt = true;
if ( jt->second ) secondOldInt = true;
}
}
//if one of the partons was not tagged before, but got tagged by this last rescatter,
//then those new tags should be removed, otherwise nothing has to be done.
tRealPartonPtr oldInt;
if ( firstOldInt == false && doesInts.rbegin()->first ) {
if ( Debug::level > 5 ) cout << " reverting tags" << endl;
//the tags should be only from the second parton, so remove first, and readd second
//only removing ancestors of first does not work, since they may share ancestors
removeInt(getReal(intDip->partons().first), intDip);
addInt(getReal(intDip->partons().second), intDip);
}
else if ( secondOldInt == false && doesInts.rbegin()->second ) {
if ( Debug::level > 5 ) cout << " reverting tags" << endl;
removeInt(getReal(intDip->partons().second), intDip);
addInt(getReal(intDip->partons().first), intDip);
}
if ( Debug::level > 5 ) cout << "done in removeLastRescatter" << endl;
}
void RealPartonState::removeInt(tRealPartonPtr rp, DipolePtr intDip) {
rp->interactions.erase(intDip);
if ( rp->oMothers.first ) removeInt(rp->oMothers.first, intDip);
if ( rp->oMothers.second ) removeInt(rp->oMothers.second, intDip);
if ( rp->oMother ) removeInt(rp->oMother, intDip);
}
void RealPartonState::addInt(tRealPartonPtr rp, DipolePtr intDip) {
rp->interactions.insert(intDip);
if ( rp->oMothers.first ) addInt(rp->oMothers.first, intDip);
if ( rp->oMothers.second ) addInt(rp->oMothers.second, intDip);
if ( rp->oMother ) addInt(rp->oMother, intDip);
}
Energy RealPartonState::neededValenceMinus() {
if ( interactions.empty() ) cerr << "neededvalenceminus called for a state "
" without interaction!" << endl;
Energy ret = minusDeficit;
minusDeficit = ZERO;
for ( RealPartonSet::iterator it = valence.begin();
it != valence.end(); it++ ) {
ret += (**it).giveMinus();
}
return ret;
}
void RealPartonState::setOnShell(tDipolePtr intDip) {
getReal(intDip->partons().first)->setOnShell();
getReal(intDip->partons().second)->setOnShell();
}
void RealPartonState::removeBackwards() {
for ( RealPartonSet::const_iterator it = backwardsPartons.begin();
it != backwardsPartons.end(); it++ ) {
}
}
Energy RealPartonState::totalMinus() {
Energy ret = ZERO;
for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = partons.begin();
it != partons.end(); it++) {
tRealPartonPtr rp = (*it).second;
if ( rp->keep == RealParton::YES ) {
ret += rp->minus;
}
}
return ret;
}
Energy RealPartonState::totalPlus() {
Energy ret = ZERO;
for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = partons.begin();
it != partons.end(); it++) {
tRealPartonPtr rp = (*it).second;
if ( rp->keep == RealParton::YES ) {
ret += rp->plus;
}
}
return ret;
}
void RealPartonState::addRecoilers() {
if ( Debug::level > 5 ) {
cout << "before recoilers" << endl;
// plotState(true);
}
//loop over all active partons
for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = partons.begin();
it != partons.end(); it++) {
if ( (*it).second->keep == RealParton::NO ) continue;
//hit holes in the parton from the recoils
addRecoiler2((*it).second);
}
if ( Debug::level > 5 ) {
cout << "with recoilers" << endl;
// plotState(true);
}
}
void RealPartonState::addRecoiler2(RealPartonPtr rp) {
TransverseMomentum opt = rp->theParton->valencePT();
if ( rp->mothers.first ) opt += rp->theParton->recoil(rp->mothers.first->theParton);
if ( rp->mothers.second ) opt += rp->theParton->recoil(rp->mothers.second->theParton);
if ( rp->nMothers == 1 && rp->mother ) opt = rp->theParton->valencePT() +
rp->theParton->recoil(rp->mother->theParton);
if ( opt.pt() == ZERO ) opt = rp->theParton->pT();
TransverseMomentum rec = rp->pT - opt;
if ( rec.pt() > 2.0*opt.pt() ) {
double plusRatio = opt.pt()/rec.pt();
rp->emitRecoiler(rec, plusRatio);
if ( Debug::level > 5 ) {
cout << "added recoiler for " << rp->oY() << endl;
// plotState(true);
}
}
}
void RealPartonState::addRecoiler(RealPartonPtr rp) {
PartonPtr p = rp->theParton;
//loop over all children.
//if any first children, create extra gluon in both real and dipole state
if ( !rp->children.first.empty() || rp->firstInt ) {
if ( Debug::level > 5 ) {
cout << "try to create first recoiler to " << rp->oY() << endl;
}
//set up
DipolePtr dip;
if ( p->dipoles().first ) dip = p->dipoles().first;
else {
Throw<RealPartonKinematicException>()
<< "Quark wants to create Recoiler on the wrong side!" << Exception::warning;
return;
}
TransverseMomentum pt = TransverseMomentum();
double ymax = rp->oY();
//move pt from rp to the new recoiler
if ( !rp->children.first.empty() ) {
for ( RealPartonSet::iterator it = rp->children.first.begin();
it != rp->children.first.end(); it++ ) {
RealPartonPtr child = *it;
//loop through "recoils" and move all pt recoils to the recoiler
for ( list< RealParton::Recoil >::iterator jt = child->recoils.begin();
jt != child->recoils.end(); jt++ ) {
cout << child->oY() << " has a recoil of " << jt->second.first.pt()/GeV
<< " to " << jt->first->oY() << endl;
if ( jt->first != rp ) continue;
if ( ymax < (*it)->y ) ymax = (*it)->y;
pt -= jt->second.first;
if ( Debug::level > 5 )
cout << child->oY() << " gives a pt kick of " << jt->second.first.pt()/GeV
<< " to " << rp->oY() << endl;
if ( Debug::level > 5 )
cout << child->oY() << " recoilers pt is now " << pt.pt()/GeV << endl;
}
// if ( child->nMothers == 1 ) child->mother = realRecoiler;
// if (child->nMothers == 2 ) child->mothers.second = realRecoiler;
}
}
if ( rp->firstInt ) {
vector<RealPartonPtr> ints;
//do the same, but use intDist to set the recoil
for ( list< RealParton::Recoil >::iterator it = rp->recoils.begin();
it != rp->recoils.end(); it++ ) {
if ( it->first->realState != this ) {
ints.push_back(it->first);
}
for (int i = 0; i < int(ints.size()); i++ ) {
//loop through "recoils" and move all pt recoils to the recoiler
for ( list< RealParton::Recoil >::iterator jt = ints[i]->recoils.begin();
jt != ints[i]->recoils.end(); jt++ ) {
cout << -ints[i]->oY() << " has an int recoil of " << jt->second.first.pt()/GeV
<< " to " << jt->first->oY() << endl;
if ( jt->first != rp ) continue;
if ( ymax < -ints[i]->y ) ymax = -ints[i]->y;
pt -= jt->second.first;
if ( Debug::level > 5 )
cout << ints[i]->oY() << " from other state gives a pt kick of " << jt->second.first.pt()/GeV
<< " to " << rp->oY() << endl;
if ( Debug::level > 5 )
cout << ints[i]->oY() << " recoilers pt is now " << pt.pt()/GeV << endl;
}
}
}
//maybe take direction from current pt?
//can "recoils" be used? does it store int recoils?
}
//if the hole is smaller than the mother, kick out a recoiler
if ( pt.pt() > (rp->pT - pt).pt() ) {
//create gluons
PartonPtr recoiler = new_ptr(Parton());
recoiler->onShell(true);
dip->generatedGluon(recoiler);
//fix colour flow for dipole state (let it be connected to the original one)
dip->splitDipole(0.5); //use better approximation for colour choice?
//initialise real recoiler
RealPartonPtr realRecoiler = getReal(recoiler);
realRecoiler->nMothers = 1;
realRecoiler->setOMother(rp);
realRecoiler->mother = rp;
realRecoiler->keep = RealParton::YES;
//set momentum of recoiler
realRecoiler->pT = pt;
rp->pT -= pt;
//otherwise use ymax to limit how far the recoiler can go in rapidity
double plusRatio = rp->pT.pt()/(realRecoiler->pT.pt() + rp->pT.pt());
realRecoiler->plus = plusRatio*rp->plus;
rp->plus = (1.0 - plusRatio)*rp->plus;
realRecoiler->updateYMinus();
rp->updateYMinus();
recoiler->oY(realRecoiler->y);
rp->saveState();
realRecoiler->saveState();
if ( Debug::level > 5 )
cout << "recoiler takes plus ratio " << plusRatio << endl;
//set position of recoiler
recoiler->position(p->position() + p->pTScale()*pt/sqr(pt.pt()));
//set mother structure
if ( !rp->firstInt ) {
RealPartonPtr lastChild = *rp->children.first.rbegin();
if ( lastChild->nMothers == 1 ) lastChild->mother = realRecoiler;
if ( lastChild->nMothers == 2 ) lastChild->mothers.second = realRecoiler;
realRecoiler->children.first.insert(lastChild);
}
else {
realRecoiler->interacting = rp->interacting;
realRecoiler->intDist = rp->intDist;
realRecoiler->firstInt = true;
}
if ( Debug::level > 5 ) {
cout << "recoiler places at transverse distance " << p->pTScale()*rp->pT.pt()/sqr(rp->pT.pt())*GeV << endl;
cout << "recoiler got pt " << realRecoiler->pT.pt()/GeV << endl;
cout << "done with first recoiler of " << rp->oY() << endl;
// dip->dipoleState().plotState(false);
// plotState(true);
}
}
else {
if ( Debug::level > 5 ) {
cout << "not enough pt, dont do recoiler" << endl;
cout << "recoiler pt was " << pt.pt()/GeV << ", canceled parent pt "
<< (rp->pT - pt).pt()/GeV << endl;
// plotState(true);
}
}
}
//same for other side
//if any first children, create extra gluon in both real and dipole state
if ( !rp->children.second.empty() || rp->secondInt ) {
if ( Debug::level > 5 ) {
cout << "try to create second recoiler to " << rp->oY() << endl;
}
//set up
DipolePtr dip;
if ( p->dipoles().second ) dip = p->dipoles().second;
else {
Throw<RealPartonKinematicException>()
<< "Quark wants to create Recoiler on the wrong side!" << Exception::warning;
return;
}
TransverseMomentum pt = TransverseMomentum();
double ymax = rp->oY();
//move pt from rp to the new recoiler
if ( !rp->children.second.empty() ) {
for ( RealPartonSet::iterator it = rp->children.second.begin();
it != rp->children.second.end(); it++ ) {
RealPartonPtr child = *it;
if ( rp->children.first.find(child) != rp->children.first.end() ) {
cout << "found second child at " << child->oY()
<< " that is also first. skip!" << endl;
plotState(true);
continue;
}
//loop through "recoils" and move all pt recoils to the recoiler
for ( list< RealParton::Recoil >::iterator jt = child->recoils.begin();
jt != child->recoils.end(); jt++ ) {
cout << child->oY() << " has a recoil of " << jt->second.first.pt()/GeV
<< " to " << jt->first->oY() << endl;
if ( jt->first != rp ) continue;
if ( ymax < (*it)->y ) ymax = (*it)->y;
pt -= jt->second.first;
if ( Debug::level > 5 )
cout << child->oY() << "gives a pt kick of " << jt->second.first.pt()/GeV
<< " to " << rp->oY() << endl;
if ( Debug::level > 5 )
cout << child->oY() << "recoilers pt is now " << pt.pt()/GeV << endl;
}
// if ( child->nMothers == 1 ) child->mother = realRecoiler;
// if (child->nMothers == 2 ) child->mothers.second = realRecoiler;
}
}
if ( rp->secondInt ) {
//do the same, but use intDist to set the recoil
//maybe take direction from current pt?
//can "recoils" be used? does it store int recoils?
}
//if the hole is smaller than the mother, kick out a recoiler
if ( pt.pt() > (rp->pT - pt).pt() ) {
//create gluons
PartonPtr recoiler = new_ptr(Parton());
recoiler->onShell(true);
dip->generatedGluon(recoiler);
//fix colour flow for dipole state (let it be connected to the original one)
dip->splitDipole(0.5); //use better approximation for colour choice?
//initialise real recoiler
RealPartonPtr realRecoiler = getReal(recoiler);
realRecoiler->nMothers = 1;
realRecoiler->setOMother(rp);
realRecoiler->mother = rp;
realRecoiler->keep = RealParton::YES;
//set momentum of recoiler
realRecoiler->pT = pt;
rp->pT -= pt;
//otherwise use ymax to limit how far the recoiler can go in rapidity
double plusRatio = rp->pT.pt()/(realRecoiler->pT.pt() + rp->pT.pt());
realRecoiler->plus = plusRatio*rp->plus;
rp->plus = (1.0 - plusRatio)*rp->plus;
realRecoiler->updateYMinus();
rp->updateYMinus();
recoiler->oY(realRecoiler->y);
rp->saveState();
realRecoiler->saveState();
if ( Debug::level > 5 )
cout << "recoiler takes plus ratio " << plusRatio << endl;
//set position of recoiler
recoiler->position(p->position() + p->pTScale()*pt/sqr(pt.pt()));
//set mother structure
if ( !rp->secondInt ) {
RealPartonPtr lastChild = *rp->children.second.rbegin();
if ( lastChild->nMothers == 1 ) lastChild->mother = realRecoiler;
if ( lastChild->nMothers == 2 ) lastChild->mothers.first = realRecoiler;
realRecoiler->children.second.insert(lastChild);
}
else {
realRecoiler->interacting = rp->interacting;
realRecoiler->intDist = rp->intDist;
realRecoiler->secondInt = true;
}
if ( Debug::level > 5 ) {
cout << "recoiler places at transverse distance " << p->pTScale()*rp->pT.pt()/sqr(rp->pT.pt())*GeV << endl;
cout << "recoiler got pt " << realRecoiler->pT.pt()/GeV << endl;
cout << "done with first recoiler of " << rp->oY() << endl;
// dip->dipoleState().plotState(false);
// plotState(true);
}
}
else {
if ( Debug::level > 5 ) {
cout << "not enough pt, dont do recoiler" << endl;
cout << "recoiler pt was " << pt.pt()/GeV << ", canceled parent pt "
<< (rp->pT - pt).pt()/GeV << endl;
// plotState(true);
}
}
}
cout << "done with " << rp->oY() << endl;
plotState(true);
}
void RealPartonState::plotState(bool pause) const {
cout << "print state to realState.dat" << endl;
ofstream file ("realState.dat");
plotStateToFile(file, false);
file.close();
if ( pause && Debug::level > 5 ) {
cout << "press any key to continue..." << endl;
getchar();
}
}
void RealPartonState::plotBothStates(RealPartonStatePtr rrs, bool pause) const {
cout << "print state to realState.dat" << endl;
ofstream file ("realState.dat");
plotStateToFile(file, false);
rrs->plotStateToFile(file, true);
file.close();
if ( pause && Debug::level > 5 ) {
cout << "press any key to continue..." << endl;
getchar();
}
}
void RealPartonState::plotStateToFile(ostream & file, bool m) const {
double ptmax = -100;
double ptmin = 100000;
double ymax = 0.0;
PartonPtr p = partons.begin()->second->theParton;
if ( p->dipoles().first ) ymax = p->dipoles().first->dipoleState().ymax();
if ( p->dipoles().second ) ymax = p->dipoles().second->dipoleState().ymax();
for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = partons.begin();
it != partons.end(); it++) {
RealPartonPtr rp = it->second;
if ( rp->pT.pt()/GeV > ptmax ) ptmax = rp->pT.pt()/GeV;
if ( rp->pT.pt()/GeV < ptmin ) ptmin = rp->pT.pt()/GeV;
bool val = rp->theParton->valence();
TransverseMomentum opt = rp->theParton->valencePT();
if ( rp->mothers.first ) opt += rp->theParton->recoil(rp->mothers.first->theParton);
if ( rp->mothers.second ) opt += rp->theParton->recoil(rp->mothers.second->theParton);
if ( rp->nMothers == 1 && rp->mother ) opt = rp->theParton->valencePT() +
rp->theParton->recoil(rp->mother->theParton);
if ( opt.pt() == ZERO ) opt = rp->theParton->pT();
//the main real partons
file << ((!m) ? rp->theParton->oY():-rp->theParton->oY()) << '\t'
<< rp->theParton->position().y()*GeV << '\t'
<< rp->keep << '\t'
<< ((!m) ? rp->y:-rp->y) << '\t'
<< val << '\t'
<< opt.pt()/GeV << '\t'
<< bool(toCheck.find(rp) != toCheck.end()) << '\t'
<< rp->theParton->position().x()*GeV << endl << endl;
//info on the consistent state, N/A if state never saved.
if ( interactions.size() > 1 || cTotalRecoil == totalRecoil ) {
file << ((!m) ? rp->theParton->oY():-rp->theParton->oY()) << '\t'
<< opt.pt()/GeV << '\t'
<< 3 << '\t' << endl;
file << ((!m) ? rp->theParton->y():abs(rp->theParton->y())) << '\t'
<< rp->theParton->pT().pt()/GeV << '\t'
<< 4 << '\t' << endl << endl;
}
//the interacting marks
if ( rp->interacting )
file << ((!m) ? rp->theParton->oY():-rp->theParton->oY()) << '\t'
<< rp->theParton->position().y()*GeV << '\t'
<< 2 << '\t'
<< ((!m) ? rp->y:-rp->y) << '\t'
<< val << '\t'
<< opt.pt()/GeV << '\t'
<< rp->theParton->position().x()*GeV << endl
<< ymax << '\t'
<< rp->theParton->position().y()*GeV << '\t'
<< 2 << '\t'
<< ((!m) ? rp->y:-rp->y) << '\t'
<< val << '\t'
<< rp->theParton->pTScale()/rp->intDist/GeV << '\t'
<< rp->theParton->position().x()*GeV << endl << endl;
//line to current state
file << ((!m) ? rp->theParton->oY():-rp->theParton->oY()) << '\t'
<< rp->theParton->position().y()*GeV << '\t'
<< rp->keep << '\t'
<< 7 << '\t'
<< opt.pt()/GeV << endl
<< ((!m) ? rp->y:-rp->y) << '\t'
<< rp->theParton->position().y()*GeV << '\t'
<< rp->keep << '\t'
<< 8 << '\t'
<< rp->pT.pt()/GeV << endl << endl;
//the pT
file << ((!m) ? rp->theParton->oY():-rp->theParton->oY()) << '\t'
<< rp->theParton->position().y()*GeV << '\t'
<< ((rp->keep == RealParton::YES) ? 9:10) << '\t'
<< rp->theParton->position().x()*GeV << '\t'
<< 0.0 << '\t'
<< 0.0 << endl
<< ((!m) ? rp->theParton->oY():-rp->theParton->oY()) << '\t'
<< rp->theParton->position().y()*GeV << '\t'
<< ((rp->keep == RealParton::YES) ? 9:10) << '\t'
<< rp->theParton->position().x()*GeV << '\t'
<< rp->pT.x()/GeV << '\t'
<< rp->pT.y()/GeV << endl << endl;
vector<RealPartonPtr> mots;
mots.push_back(rp->oMothers.first);
mots.push_back(rp->oMothers.second);
mots.push_back(rp->mothers.first);
mots.push_back(rp->mothers.second);
mots.push_back(rp->oMother);
mots.push_back(rp->mother);
//lines to mothers and original mothers
for ( int i = 0; i < int(mots.size()); i++) {
RealPartonPtr rp2 = rp;
if ( mots[i] ) rp2 = mots[i];
TransverseMomentum opt2 = rp2->theParton->valencePT();
if ( rp2->mothers.first ) opt2 += rp2->theParton->recoil(rp2->mothers.first->theParton);
if ( rp2->mothers.second ) opt2 += rp2->theParton->recoil(rp2->mothers.second->theParton);
if ( rp2->nMothers == 1 && rp2->mother ) opt2 = rp2->theParton->valencePT() +
rp2->theParton->recoil(rp2->mother->theParton);
if ( opt2.pt() == ZERO ) opt2 = rp2->theParton->pT();
file << ((!m) ? rp->theParton->oY():-rp->theParton->oY()) << '\t'
<< rp->theParton->position().y()*GeV << '\t'
<< rp->keep << '\t'
<< i+3 + 100 << '\t'
<< opt.pt()/GeV << '\t'
<< rp->theParton->position().x()*GeV << '\t'
<< rp->nMothers + 10 << endl
<< ((!m) ? rp2->theParton->oY():-rp2->theParton->oY()) << '\t'
<< rp2->theParton->position().y()*GeV << '\t'
<< rp2->keep << '\t'
<< i+3 + 100 << '\t'
<< opt2.pt()/GeV << '\t'
<< rp2->theParton->position().x()*GeV << '\t'
<< rp->nMothers + 10 << endl;
}
}
file << ymax << '\t' << ptmax << '\t' << 10 << '\t' << 10 << endl
<< ymax << '\t' << ptmin << '\t' << 10 << '\t' << 10 << endl;
file << ((!m) ? currentY:-currentY) << '\t' << ptmax << '\t' << 10 << '\t' << 9 << endl
<< ((!m) ? currentY:-currentY) << '\t' << ptmin << '\t' << 10 << '\t' << 9 << endl;
}
pair<int, int> RealPartonState::countYESNO() const {
int yes = 0;
int no = 0;
for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = partons.begin();
it != partons.end(); it++) {
if ( (*it).second->keep == RealParton::YES ) yes++;
else no++;
}
return make_pair(yes, no);
}
double RealPartonState::avYInEvo() const {
double number = 0.0;
double length = 0.0;
for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = partons.begin();
it != partons.end(); it++) {
RealPartonPtr rp = (*it).second;
if ( rp->keep == RealParton::YES && !(rp->theParton->valence()) ) {
number++;
if ( rp->nMothers == 2 )
length += (abs(rp->y - rp->mothers.first->y) + abs(rp->y - rp->mothers.second->y))/2.0;
if ( rp->nMothers == 1 )
length += abs(rp->y - rp->mother->y);
}
}
if (number == 0.0 ) return 0.0;
return length/number;
}
bool RealPartonState::diagnosis(bool pause) const {
calls++;
bool ok = true;
// cout << "diagnosis call number " << calls << " for " << this << endl;
TransverseMomentum pT = TransverseMomentum();
Energy dplus = ZERO;
Energy dminus = ZERO;
TransverseMomentum cpT = TransverseMomentum();
Energy cplus = ZERO;
Energy cminus = ZERO;
Energy gminus = ZERO;
int yes = 0;
int consInt = 0;
int vals = 0;
int mothers = 0;
int kids = 0;
int exchanges = 0;
int interactingPartons = 0;
TransverseMomentum valPT = TransverseMomentum();
TransverseMomentum intPT = TransverseMomentum();
Energy valPlus = ZERO;
Energy valMinus = ZERO;
for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = partons.begin();
it != partons.end(); it++) {
tRealPartonPtr rp = (*it).second;
if ( !rp->exchanges.empty() )
exchanges += rp->exchanges.size();
if ( rp->cInteracting ) consInt++;
if ( rp->theParton->valence() ) {
vals++;
valPT += rp->theParton->valencePT();
valPlus += rp->theParton->valencePlus();
if ( rp->movedValence )
valMinus += rp->theParton->valencePlus()*exp(2.0*rp->movedValence->oY());
else
valMinus += rp->theParton->valencePlus()*exp(2.0*rp->oY());
if ( rp->keep == RealParton::NO ) {
pT += rp->theParton->valencePT();
dplus += rp->theParton->valencePlus();
dminus += sqr(rp->theParton->valencePT().pt())/rp->theParton->valencePlus();
cout << "* valence parton that is NO!!" << endl;
ok = false;
}
}
if ( rp->keep == RealParton::YES ) {
yes++;
pT += rp->pT;
dplus += rp->plus;
dminus += rp->minus;
gminus += rp->givenMinus;
if ( rp->mothers.first ) mothers++;
if ( rp->mothers.second ) mothers++;
if ( rp->mother ) mothers++;
kids += rp->children.first.size();
kids += rp->children.second.size();
if ( rp->interacting ) interactingPartons++;
}
if ( rp->cKeep == RealParton::YES ) {
cpT += rp->theParton->pT();
cplus += rp->theParton->plus();
cminus += rp->theParton->minus();
}
if ( rp->interacting ) {
intPT += rp->intRecoil;
}
// if ( rp->interacting && toCheck.find(rp) == toCheck.end() ) {
// cout << "interacting parton at " << rp->theParton->oY() << " is not in tocheck" << endl;
// ok = false;
// }
}
for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = partons.begin();
it != partons.end(); it++) {
tRealPartonPtr rp = (*it).second;
if ( rp->keep != RealParton::NO && rp->keep != RealParton::YES )
Throw<RealPartonKinematicException>()
<< "Real parton at oY " << rp->oY() << " is neither YES nor NO..." << Exception::warning;
if ( rp->keep == RealParton::YES ) {
// if ( !rp->checkMomentum() && !(rp->interacting && totalRecoil.pt() != ZERO )
// && rp->theParton->oY() <= currentY )
// ok = false;
if ( rp->plus < ZERO ) {
Throw<RealPartonKinematicException>()
<< "negative plus in real parton diagnosis at oY " << rp->oY() << Exception::warning;
ok = false;
}
if ( rp->minus < ZERO ) {
Throw<RealPartonKinematicException>()
<< "negative minus in real parton diagnosis at oY " << rp->oY() << Exception::warning;
ok = false;
}
for ( RealPartonSet::const_iterator it = rp->children.first.begin();
it != rp->children.first.end(); it++ ) {
if ( (*it)->mothers.second != rp && (*it)->mother != rp ) {
cout << "\n* parton at " << (*it)->theParton->oY() << " doesnt have " << rp->theParton->oY()
<< " as second mother! child-mother broken!" << endl;
ok = false;
}
if ( (*it)->keep == RealParton::NO ) {
cout << "\n* parton at " << rp->theParton->oY() << "has a NO-child at "
<< (*it)->theParton->oY() << endl;
ok = false;
}
}
for ( RealPartonSet::const_iterator it = rp->children.second.begin();
it != rp->children.second.end(); it++ ) {
if ( (*it)->mothers.first != rp && (*it)->mother != rp ) {
cout << "\n* parton at " << (*it)->theParton->oY() << " doesnt have " << rp->theParton->oY()
<< " as first mother! child-mother broken!" << endl;
ok = false;
}
if ( (*it)->keep == RealParton::NO ) {
cout << "\n* parton at " << rp->theParton->oY() << "has a NO-child at "
<< (*it)->theParton->oY() << endl;
ok = false;
}
}
if ( rp->mothers.first && (rp->mothers.first->children.second.find(rp)
== rp->mothers.first->children.second.end()) ) {
cout << "\n* first mother of " << rp->oY() << " doesnt have this as second child!" << endl;
ok = false;
}
if ( rp->mothers.second && (rp->mothers.second->children.first.find(rp)
== rp->mothers.second->children.first.end()) ) {
cout << "\n* second mother of " << rp->oY() << " doesnt have this as first child!" << endl;
ok = false;
}
if ( rp->mother &&
(rp->mother->children.first.find(rp) == rp->mother->children.first.end()) &&
(rp->mother->children.second.find(rp) == rp->mother->children.second.end()) ) {
cout << "\n* single mother of " << rp->oY()
<< " doesnt have this as first or second child!" << endl;
ok = false;
}
if ( !rp->interacting && !rp->theParton->valence() &&
rp->children.first.empty() && rp->children.second.empty() ) {
cout << "\n* " << rp->oY() << " has no kids, and isnt interacting! where to get p-??\n";
ok = false;
}
}
}
if ( (intPT - totalRecoil).pt() > 0.000000001*GeV ) {
cout << "* sum of intPT != totalRecoil" << endl;
ok = false;
}
if ( (pT - totalRecoil).pt() > 0.000000001*GeV ) {
cout << "* total pT != totalRecoil" << endl;
ok = false;
}
// if ( valence.size() != 3 ) {
// cout << "* size of valence is not 3" << endl;
// ok = false;
// }
if ( valPT.pt() > 0.000000001*GeV ) {
cout << "* valencePT does dot add up to zero!" << endl;
ok = false;
}
if ( vals != int(valence.size()) ) {
cout << "* there are " << vals << " valence partons, but size of valence is "
<< valence.size() << endl;
ok = false;
}
if ( kids != mothers ) {
cout << "* not same number of kids as mothers!" << endl;
ok = false;
}
if ( abs(valPlus - dplus) > 0.00000001*GeV && totalRecoil.pt() == ZERO ) {
cout << "* initial state plus not conserved" << endl;
ok = false;
}
if ( dplus == ZERO ) {
cout << "* no plus!" << endl;
ok = false;
}
if ( isnan(dplus/GeV) || isnan(dminus/GeV) ) {
cout << "* plus or minus is nan!!" << endl;
ok = false;
}
if ( !ok && Debug::level > 5 ) {
cout << "*\n* REAL STATE NOT OK!!!\n*\n";
}
if ((!ok || pause) && Debug::level > 5 ) {
cout << setprecision(10);
cout << "********* Diagnosis for " << this << " ******* call nr " << calls << " ********" << endl;
cout << "* number of partons: " << partons.size() << endl;
cout << "* size of valence: " << valence.size() << endl;
cout << "* number of valence part: " << vals << endl;
cout << "* number of toCheck: " << toCheck.size() << endl;
cout << "* number of interacting: " << interacting.size() << endl;
cout << "* number of int rps: " << interactingPartons << endl;
cout << "* number of cInteracting: " << consInt << endl;
cout << "* number of interactions: " << interactions.size() << endl;
cout << "* number gluon exchanges: " << exchanges << endl;
cout << "* consistent: " << consistent << endl;
cout << "* total no of kids = " << kids << endl;
cout << "* total no mothers = " << mothers << endl;
cout << "* number of YES or INT: " << yes << endl
<< "* " << endl;
cout << "* total transverse = (" << pT.x()/GeV << ", " << pT.y()/GeV << ")" << endl;
cout << "* total recoil = (" << totalRecoil.x()/GeV << ", "
<< totalRecoil.y()/GeV << ")" << endl;
cout << "* sum of int recoil = (" << intPT.x()/GeV << ", "
<< intPT.y()/GeV << ")" << endl;
cout << "* original plus = " << plus/GeV << endl;
cout << "* total plus = " << dplus/GeV << endl;
cout << "* 1800 - total plus + " << 1800.0 - dplus/GeV << endl;
cout << "* original minus = " << minus/GeV << endl;
cout << "* total minus = " << dminus/GeV << endl;
cout << "* given +val minus = " << (gminus + valMinus)/GeV << endl;
cout << "* valence minus = " << valMinus/GeV << endl;
cout << "* total ctransverse = (" << cpT.x()/GeV << ", " << cpT.y()/GeV << ")" << endl;
cout << "* total cplus = " << cplus/GeV << endl;
cout << "* total cminus = " << cminus/GeV << endl;
cout << "* total valencePT = (" << valPT.x()/GeV << ", " << valPT.y()/GeV << ")" << endl;
cout << "* total valencePlus = " << valPlus/GeV << endl;
cout << "* number of flucts = " << flucts.size() << endl;
for ( map<tPartonPtr,RealPartonPtr>::const_iterator it = partons.begin();
it != partons.end(); it++) {
tRealPartonPtr rp = (*it).second;
if ( rp->fluct != -1 ) {
cout << "* " << rp->oY() << " belongs to fluct no " << rp->fluct << endl;
}
// cout << "* rp at " << rp->theParton->oY() << " has vetoed p+1: " << rp->plus1veto
// << ", p+2: " << rp->plus2veto << ", p-1: " << rp->minus1veto << ", p-2: " << rp->minus2veto
// << ", DGLAP: " << rp->unDGLAP << endl;
// if ( rp->theParton->valence() )
// cout << "* valencepT: (" << rp->theParton->valencePT().x()/GeV << ", "
// << rp->theParton->valencePT().y()/GeV
// << "), valenceplus: " << rp->theParton->valencePlus()/GeV << endl;
}
cout << "*********************************************************" << endl;
plotState(pause || !ok);
}
if ( Debug::level > 5 ) cout << "real state tested. ok: " << ok << endl;
return ok;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 2:36 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804771
Default Alt Text
(266 KB)

Event Timeline