Page MenuHomeHEPForge
Files F8723647

No OneTemporary
Actions

View Options

diff --git a/.hgignore b/.hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -1,46 +1,51 @@
.*/Makefile
.*/Makefile\.in
.*/\.deps
.*/\.libs
.*/.*\.l[ao]
.*/.*\.so\.*
.*/.*\.o
.*~
.*/done-all-links
lib/AriadneDefaults.rpo
(DIPSY|src)/.*\.(run|tex|out|log|rpo|spc|top|dump|dot|aux|pdf|ps|png|svg|hepmc|dat|aida|rz|eps|root|tuple|gp|exe)
src/done-all-links
autom4te.cache
aclocal.m4
libtool
Makefile
Makefile.in
config.log
config.status
configure
Config/config.h
Config/config.h.in
Config/config.sub
Config/depcomp
Config/install-sh
Config/missing
Config/stamp-h.
Config/config.guess
include/done-all-links
include/Ariadne
DIPSY/NuclearDistribution-1.C
DIPSY/NuclearDistribution.C
DIPSY/plots
src/3jet.lhe.gz
src/MadGraph_DIS
src/mgdis.in
src/runThePEG
DIPSY/runThePEG
lib/Ariadne5Defaults.rpo
DIPSY/fsswing.txt
DIPSY/rphi.pl
DIPSY1991
DIPSY/mcnetRivetbin
DIPSY/pythia8
DIPSY/exps
DIPSY/temp.*
+DIPSY/95
+DIPSY/.*\.yoda
+DIPSY/.*\.info
+DIPSY/.*\.plot
+
diff --git a/Cascade/AriadneHandler.cc b/Cascade/AriadneHandler.cc
--- a/Cascade/AriadneHandler.cc
+++ b/Cascade/AriadneHandler.cc
@@ -1,653 +1,697 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the AriadneHandler class.
//
#include "AriadneHandler.h"
#include "ConsistencyChecker.h"
#include "QCDDipoleFinder.h"
#include "EMDipoleFinder.h"
#include "ReweightBase.h"
#include "DipoleState.h"
#include "EmitterBase.h"
#include "BornCheckerBase.h"
#include "ScaleSetter.h"
#include "DISFinder.h"
#include "ResonanceFinder.h"
#include "Ariadne/Cascade/Models/RemnantModel.h"
#include "Ariadne/Cascade/Models/ColourResonanceModel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Collision.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Handlers/Hint.h"
#include "ThePEG/Handlers/XComb.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Utilities/DescriptionList.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "ThePEG/Utilities/Current.h"
+#include "ThePEG/Utilities/DebugItem.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/DescribeClass.h"
+#include "Models/FSGluonEmission.h"
+#include "Models/FSQQEmission.h"
+#include "Models/DipoleSwing.h"
+#include "ThePEG/Analysis/FactoryBase.h"
+#ifndef LWH_AIAnalysisFactory_H
+#ifndef LWH
+#define LWH ThePEGLWH
+#endif
+#include "ThePEG/Analysis/LWH/AnalysisFactory.h"
+#endif
+
+
using namespace Ariadne5;
AriadneHandler::AriadneHandler()
: theRunningCoupling(simpleRunning), theAlpha0(0.2), theLambdaQCD(0.22*GeV),
scaleFactor(1.0), thePTCut(0.6*GeV), theAlphaEM0(1.0/137.0),
thePTCutEM(0.6*GeV), theNCol(8), theNFlav(5), thePhotonEmissions(false),
theSoftMu(0.6*GeV), theSoftAlpha(1.0), theHardAlpha(-1.0), theBeta(2.0),
theMaxEmissions(0), suspendConsistencyChecks(false) {}
AriadneHandler::~AriadneHandler() {}
IBPtr AriadneHandler::clone() const {
return new_ptr(*this);
}
IBPtr AriadneHandler::fullclone() const {
return new_ptr(*this);
}
bool AriadneHandler::preInitialize() const {
if ( CascadeHandler::preInitialize() ) return true;
return runningCoupling() == internalRunning && !internalAlphaS();
}
void AriadneHandler::doinit() throw(InitException) {
CascadeHandler::doinit();
if ( runningCoupling() != internalRunning || internalAlphaS() ) return;
theInternalAlphaS = dynamic_ptr_cast<Ptr<AlphaSBase>::pointer>
(generator()->preinitCreate("ThePEG::O1AlphaS", fullName() + "/AlphaS",
"O1AlphaS.so"));
ostringstream os;
os << lambdaQCD()/GeV;
generator()->preinitInterface(theInternalAlphaS,
"LambdaQCD", "set", os.str());
generator()->preinitInterface(theInternalAlphaS, "LambdaFlav", "set", "5");
theInternalAlphaS->update();
// Remove duplicate emitters.
vector<EmitterPtr> cleaned;
set<const ClassDescriptionBase *> cleanset;
for ( int i = 0, N = emitters().size(); i < N; ++i ) {
const ClassDescriptionBase * cd =
DescriptionList::find(typeid(*emitters()[i]));
if ( cleanset.find(cd) != cleanset.end() ) continue;
cleaned.push_back(emitters()[i]);
cleanset.insert(cd);
}
theEmitters.swap(cleaned);
}
-inline void AriadneHandler::doinitrun() {
+void AriadneHandler::doinitrun() {
+ static DebugItem histem("Ariadne5::HistEm", 6);
CascadeHandler::doinitrun();
theFlavourThresholds.clear();
vector<Energy2> thrsh;
switch ( runningCoupling() ) {
case noRunning:
return;
case simpleRunning:
case externalRunning:
thrsh = SM().alphaSPtr()->flavourThresholds();
break;
case internalRunning:
thrsh = internalAlphaS()->flavourThresholds();
}
for ( int i = 0, N = thrsh.size(); i < N; ++i ) thrsh[i] /= scaleFactor;
theFlavourThresholds.insert(thrsh.begin(), thrsh.end());
-
+
+ if ( histem ) {
+ generator()->histogramFactory()->initrun();
+ generator()->histogramFactory()->registerClient(this);
+ generator()->histogramFactory()->mkdir("/Ariadne5Debug");
+ histall = generator()->histogramFactory()->createHistogram1D
+ ("/Ariadne5Debug/All", 100, -1.0, 9.0);
+ histswing = generator()->histogramFactory()->createHistogram1D
+ ("/Ariadne5Debug/Swings", 100, -1.0, 9.0);
+ histglue = generator()->histogramFactory()->createHistogram1D
+ ("/Ariadne5Debug/Glue", 100, -1.0, 9.0);
+ histqq = generator()->histogramFactory()->createHistogram1D
+ ("/Ariadne5Debug/QQ", 100, -1.0, 9.0);
+ histlam = generator()->histogramFactory()->createHistogram1D
+ ("/Ariadne5Debug/QQ", 100, -1.0, 9.0);
+ }
}
void AriadneHandler::cascade() {
+ static DebugItem histem("Ariadne5::HistEm", 6);
Current<AriadneHandler> current(this);
DipoleStatePtr state;
tSubProPtr sub;
tPVector final;
Energy rhomax = ZERO;
Energy rhomin = pTCut();
bool perf = false;
int emnbr = 0;
// Check if the state has already been generated in reweightCKKW
if ( theCKKWMap.count(lastXCombPtr()) ){
CKKWState ckkw = theCKKWMap[lastXCombPtr()];
sub = subProcess();
state = ckkw.history->state();
checkState(*state);
rhomax = startingScale(*state);
while ( emnbr == 0 &&
( rhomax = state->select(rhomin, rhomax) ) > rhomin ) {
SaveDipoleState backup(state);
if ( state->perform() ) {
if ( !ckkw.maxMult && checkTreeState(*state) && passCuts(state) ) {
theCKKWMap.clear();
throw Veto();
}
perf = true;
emnbr++;
} else {
state = backup.revert();
state->selected()->dipole->touch();
}
}
}
else {
tCollPtr coll = eventHandler()->currentCollision();
state = new_ptr(DipoleState(coll->incoming()));
sub = coll->primarySubProcess();
if ( sub->decayed() ) sub = tSubProPtr();
// Decide whether we are cascadeing a whole sub-process or if we
// only need to deal with final-state shower.
if ( ( hint().tagged() || !sub ) && !tagged().empty() ) {
final = tagged();
}
else if ( sub ) {
state->setup(*sub);
}
else
return;
if ( !final.empty() ) {
state->setup(set<tPPtr>(final.begin(), final.end()));
}
checkState(*state);
if ( !state->checkIntegrity() )
Throw<Exception>()
<< "Ariadne failed to setup dipole system. This is a serious error,"
<< "Please inform the author." << Exception::runerror;
rhomax = startingScale(*state);
}
theCKKWMap.clear();
while ( ( rhomax = state->select(rhomin, rhomax) ) > rhomin &&
( maxEmissions() == 0 || emnbr < maxEmissions() ) ) {
SaveDipoleState backup(state);
if ( state->perform() ) {
perf = true;
emnbr++;
+ if ( histem ) {
+ double w = 1.0/double(state->activeDipoles().size());
+ double lrho = log10(state->selected()->rho/GeV);
+ histall->fill(lrho, w);
+ histlam->fill(lrho, state->lambdaMeasure());
+ if ( dynamic_cast<DipoleSwing*>((Emission*)(state->selected())) )
+ histswing->fill(lrho, w);
+ else if ( dynamic_cast<FSGluonEmission*>((Emission*)(state->selected())) )
+ histglue->fill(lrho, w);
+ else if ( dynamic_cast<FSQQEmission*>((Emission*)(state->selected())) )
+ histqq->fill(lrho, w);
+ else
+ state->selected()->debug();
+ }
} else {
state = backup.revert();
state->selected()->dipole->touch();
}
if ( !state->checkIntegrity() ) {
throw IntegretyException() << "AriadneHandler::cascade "
<< "The dipole state is not self consistent."
<< Exception::eventerror;
}
}
if ( final.empty() ) {
if ( perf ) state->fill(*newStep());
sub->decayed(true);
} else {
if ( perf ) state->fill(*newStep());
}
}
double AriadneHandler::reweightCKKW(int minMult, int maxMult) {
if(minMult == maxMult){
return 1.0;
}
CKKWState ckkw;
DipoleStatePtr state = new_ptr(DipoleState());
tXCPtr lastXC = lastXCombPtr();
int outgoing = lastXC->subProcess()->outgoing().size();
if ( outgoing < minMult || outgoing > maxMult ) {
throw CKKWMultiplicityException()
<< "Ariadne::CascadeHandler::reweightCKKW "
<< "Number of outgoing particles out of range."
<< Exception::eventerror;
}
int steps = outgoing - minMult;
ckkw.maxMult = (outgoing == maxMult);
state->setup(*(lastXC->subProcess()));
checkState(*state);
if ( !state->checkIntegrity() )
Throw<Exception>()
<< "Ariadne failed to setup dipole system. This is a serious error,"
<< "Please inform the author." << Exception::runerror;
ckkw.history = HistoryPtr(new History(steps, state));
if ( !ckkw.history->select() ) return 0.0;
double w = ckkw.history->weight(lastXC->lastAlphaS(), lastXC->lastScale());
if ( w > 0.0 ) theCKKWMap[lastXC] = ckkw;
return w;
}
bool AriadneHandler::passCuts(tcDipoleStatePtr state) {
/* *** ATTENTION ***
tCutsPtr cuts = lastCutsPtr();
Lorentz5Momentum ph = state->hardFS().momentum();
cuts->initSubProcess(ph.mass2(), ph.rapidity());
tcPDVector pdata;
vector< LorentzMomentum > p;
typedef HardSubSys::PartonSet PartonSet;
const PartonSet & active = state->hardSubSys().coloured();
const PartonSet & produced = state->hardSubSys().produced();
for(PartonSet::const_iterator it = active.begin(); it != active.end();
it++){
pdata.push_back((*it)->dataPtr());
p.push_back((*it)->momentum());
}
for(PartonSet::const_iterator it = produced.begin(); it != produced.end();
it++){
pdata.push_back((*it)->dataPtr());
p.push_back((*it)->momentum());
}
LorentzRotation R(0.0, 0.0, - ph.z()/ph.e());
Utilities::transform(p, R);
return cuts->passCuts(pdata, p, state->particles().first->dataPtr(),
state->particles().second->dataPtr());
*/ return false;
}
double AriadneHandler::alphaS(Energy2 scale) const {
scale *= scaleFactor;
int Nf = SM().Nf(scale);
switch ( runningCoupling() ) {
case noRunning:
return alpha0();
case simpleRunning:
return 12.0*Constants::pi/
((33.0 - 2.0*min(int(SM().Nf(scale)),5))*log(scale/sqr(lambdaQCD())));
case internalRunning:
return internalAlphaS()->value(scale, SM());
case externalRunning:
return SM().alphaS(scale);
}
return ZERO*Nf;
}
Energy AriadneHandler::checkBornState(const DipoleState & ds) const {
Energy scale = ZERO;
for ( int i = 0, N = theBornCheckers.size(); i < N; ++i ) {
Energy mu = theBornCheckers[i]->check(ds);
if ( mu > ZERO ) return mu;
if ( scale == ZERO ) scale = mu;
}
if ( scale == ZERO )
Throw<CKKWBornException>()
<< "Could not reconstruct the given event in the CKKW-L algorithm. "
<< *(eventHandler()->currentEvent()) << Exception::runerror;
return scale;
}
bool AriadneHandler::checkTreeState(const DipoleState & ds) const {
for ( int i = 0, N = theBornCheckers.size(); i < N; ++i )
if ( theBornCheckers[i]->checkTree(ds) ) return true;
return false;
}
vector<tQCDPtr> AriadneHandler::findQCDDipoles(DipoleState & state) const {
return theQCDFinder->findDipoles(state);
}
vector<tEMDipPtr> AriadneHandler::findEMDipoles(DipoleState & state) const {
return theEMFinder? theEMFinder->findDipoles(state): vector<tEMDipPtr>();
}
Energy AriadneHandler::startingScale(const DipoleState & state) const {
return theScaleSetter->scale(state);
}
bool AriadneHandler::checkState(DipoleState & state, tcEmPtr e) {
if ( !consistency ) return true;
if ( !e ) suspendConsistencyChecks = false;
else if ( suspendConsistencyChecks ) return true;
if ( consistency->check(state, e) ) return true;
if ( !e ) suspendConsistencyChecks = true;
return false;
}
pair<tRemParPtr,tRemParPtr>
AriadneHandler::findDISLeptons(SubProcess & sub, DipoleState & state) const {
return theDISFinder? theDISFinder->findDISLeptons(sub, state):
pair<tRemParPtr,tRemParPtr>();
}
pair<tRemParPtr,tRemParPtr>
AriadneHandler::findDISQuarks(pair<tRemParPtr,tRemParPtr> leptons,
SubProcess & sub, DipoleState & state) const {
return theDISFinder? theDISFinder->findDISQuarks(leptons, sub, state):
pair<tRemParPtr,tRemParPtr>();
}
tPVector AriadneHandler::resonances(SubProcess & sub) const {
return theResonanceFinder->resonances(sub);
}
void AriadneHandler::persistentOutput(PersistentOStream & os) const {
os << oenum(theRunningCoupling) << theAlpha0 << ounit(theLambdaQCD, GeV)
<< theInternalAlphaS << scaleFactor << theFlavourThresholds.size()
<< ounit(thePTCut, GeV) << theAlphaEM0 << ounit(thePTCutEM, GeV)
<< theNCol << theNFlav << thePhotonEmissions << ounit(theSoftMu, GeV)
<< theSoftAlpha << theHardAlpha << theBeta << theReweighters
<< theMaxEmissions << theEmitters << theBornCheckers << theScaleSetter
<< theDISFinder << theResonanceFinder << theQCDFinder << theEMFinder
<< theColourResonanceModel << theRemnantModel
<< consistency << suspendConsistencyChecks;
for ( set<Energy2>::const_iterator it = theFlavourThresholds.begin();
it != theFlavourThresholds.end(); ++it ) os << ounit(*it, GeV2);
}
void AriadneHandler::persistentInput(PersistentIStream & is, int) {
int size = 0;
is >> ienum(theRunningCoupling) >> theAlpha0 >> iunit(theLambdaQCD, GeV)
>> theInternalAlphaS >> scaleFactor >> size
>> iunit(thePTCut, GeV) >> theAlphaEM0 >> iunit(thePTCutEM, GeV)
>> theNCol >> theNFlav >> thePhotonEmissions >> iunit(theSoftMu, GeV)
>> theSoftAlpha >> theHardAlpha >> theBeta >> theReweighters
>> theMaxEmissions >> theEmitters >> theBornCheckers >> theScaleSetter
>> theDISFinder >> theResonanceFinder >> theQCDFinder >> theEMFinder
>> theColourResonanceModel >> theRemnantModel
>> consistency >> suspendConsistencyChecks;
theFlavourThresholds.clear();
Energy2 t = ZERO;
for ( int i = 0; i < size; ++i ) {
is >> iunit(t, GeV2);
theFlavourThresholds.insert(t);
}
}
DescribeClass<AriadneHandler,ThePEG::CascadeHandler>
describeAriadne5("Ariadne5::AriadneHandler", "libAriadne5.so");
Energy AriadneHandler::minPTCut() const {
return ( runningCoupling() == simpleRunning ||
( runningCoupling() == internalRunning && !theInternalAlphaS ) )?
lambdaQCD()/sqrt(scaleFactor): 0.0*GeV;
}
Energy AriadneHandler::maxLambdaQCD() const {
return ( runningCoupling() == simpleRunning ||
( runningCoupling() == internalRunning && !theInternalAlphaS ) )?
pTCut()*sqrt(scaleFactor): Constants::MaxEnergy;
}
void AriadneHandler::Init() {
static ClassDocumentation<AriadneHandler> documentation
("The AriadneHandler class administers the Ariadne dipole "
"cascade.",
"Parton cascades performed by Ariadne\\cite{Lav05} according to the "
"Dipole Cascade Model\\cite{Gustafson:1986db,Gustafson:1988rq,"
"Andersson:1989gp,Andersson:1990ki}.",
"\\bibitem{Lav05}"
"Nils Lavesson and Leif L\\\"onnblad, Preprint in preparation.\n"
"\\bibitem{Gustafson:1986db}"
"G\\\"osta Gustafson, Phys.~Lett.~{\\bf B175} (1986) 453.\n"
"\\bibitem{Gustafson:1988rq}"
"G\\\"osta Gustafson and Ulf Pettersson, "
"Nucl.~Phys.~{\\bf B306} (1988) 746.\n"
"\\bibitem{Andersson:1989gp}"
"Bo Andersson, et al., Z.~Phys.~{\\bf C43} (1989) 625.\n"
"\\bibitem{Andersson:1990ki}"
"Bo Andersson, et al., Nucl.~Phys.~{\\bf B339} (1990) 393.");
static Switch<AriadneHandler,RunningOption> interfaceRunningCoupling
("RunningCoupling",
"Strategy for handling \\f$\\alpha_S\\f$.",
&AriadneHandler::theRunningCoupling, simpleRunning, true, false);
static SwitchOption interfaceRunningCouplingSimpleRunning
(interfaceRunningCoupling,
"SimpleRunning",
"Use a one loop running coupling with \\f$\\Lambda_{QCD}\\f$ given "
"by <interface>LambdaQCD</interface>.",
simpleRunning);
static SwitchOption interfaceRunningCouplingConstant
(interfaceRunningCoupling,
"Constant",
"Use a constant coupling given by <interface>Alpha0</interface>.",
noRunning);
static SwitchOption interfaceRunningCouplingExternalRunning
(interfaceRunningCoupling,
"ExternalRunning",
"Use whatever coupling is specified by the current StandardModelBase "
"object.",
externalRunning);
static SwitchOption interfaceRunningCouplingInternalRUnning
(interfaceRunningCoupling,
"InternalRunning",
"Use the coupling specified by <interface>InternalAlphaS</interface>.",
internalRunning);
static Parameter<AriadneHandler,double> interfaceAlpha0
("Alpha0",
"The constant \\f$\\alpha_S\\f$ to use if "
"<interface>RunningCoupling</interface> is set to Constant.",
&AriadneHandler::theAlpha0, 0.2, 0.0, 0,
true, false, Interface::lowerlim);
static Parameter<AriadneHandler,double> interfaceAlphaEM0
("AlphaEM0",
"The constant \\f$\\alpha_{EM}\\f$ to use. If zero, use whatever "
"is specified in the current StandardModelBase object.",
&AriadneHandler::theAlphaEM0, 1.0/137.0, 0.0, 0,
true, false, Interface::lowerlim);
static Parameter<AriadneHandler,Energy> interfaceLambdaQCD
("LambdaQCD",
"The \\f$\\Lambda_{QCD}\\f$ to use in the one loop running "
"\\f$\\alpha_S\\f$ if <interface>RunningCoupling</interface> "
"is set to Running.",
&AriadneHandler::theLambdaQCD, GeV, 0.22*GeV, 0.0*GeV,
Constants::MaxEnergy,
true, false, Interface::limited,
0, 0, 0, &AriadneHandler::maxLambdaQCD, 0);
static Parameter<AriadneHandler,double> interfaceScaleFactor
("ScaleFactor",
"Scale factor used to multiply the emission scales in the argument of "
"\\f$\alpha_S\\f$.",
&AriadneHandler::scaleFactor, 1.0, 0.0, 0,
true, false, Interface::lowerlim);
static Parameter<AriadneHandler,Energy> interfacePTCut
("PTCut",
"The cutoff in invariant transverse momentum for QCD emissions.",
&AriadneHandler::thePTCut, GeV, 0.6*GeV, 0.0*GeV, 0*GeV,
true, false, Interface::lowerlim,
0, 0, &AriadneHandler::minPTCut, 0, 0);
static Parameter<AriadneHandler,Energy> interfacePTCutEM
("PTCutEM",
"The cutoff in invariant transverse momentum for QED emissions.",
&AriadneHandler::thePTCutEM, GeV, 0.6*GeV, 0.0*GeV, 0*GeV,
true, false, Interface::lowerlim);
static Parameter<AriadneHandler,int> interfaceDipoleColours
("DipoleColours",
"The number of differently coloured dipoles possible. This should "
"normally be 8, but may be varied to check the effects of colour "
"reconnections.",
&AriadneHandler::theNCol, 8, 3, 0,
true, false, Interface::lowerlim);
static Parameter<AriadneHandler,int> interfaceNFlav
("NFlav",
"The number of possible flavours in a \f$g\to q\bar{q}\f$ splitting.",
&AriadneHandler::theNFlav, 5, 0, 8,
true, false, Interface::lowerlim);
static Switch<AriadneHandler,bool> interfacePhotonEmissions
("PhotonEmissions",
"Switches photon emission in the cascade on and off.",
&AriadneHandler::thePhotonEmissions, false, true, false);
static SwitchOption interfacePhotonEmissionsOn
(interfacePhotonEmissions,
"On",
"Switch photon emission on",
true);
static SwitchOption interfacePhotonEmissionsOff
(interfacePhotonEmissions,
"Off",
"Switch photon emission off.",
false);
static Parameter<AriadneHandler,Energy> interfaceSoftMu
("SoftMu",
"The inverse extension of a hadron remnant used in the soft "
"suppression mechanism. See also <interface>SoftAlphs</interface> and "
"<interface>Beta</interface>",
&AriadneHandler::theSoftMu, GeV, 0.6*GeV, 0.0*GeV, 0*GeV,
true, false, Interface::lowerlim);
static Parameter<AriadneHandler,double> interfaceSoftAlpha
("SoftAlpha",
"The dimension of the extension of a hadron remnant used in the "
"soft-suppression mechanism. See also <interface>SoftMu</interface> "
"and <interface>Beta</interface>.",
&AriadneHandler::theSoftAlpha, 1.0, 0.0, 0,
true, false, Interface::lowerlim);
static Parameter<AriadneHandler,double> interfaceHardAlpha
("HardAlpha",
"The dimension of the extension of a hard remnant used in the "
"soft-suppression mechanism for radiation off a scattered quark in "
"DIS at small \\f$Q^2\\f$. If set negative, the value of "
"<interface>SoftAlpha</interface> will be used instead. See also "
"<interface>Beta</interface>.",
&AriadneHandler::theHardAlpha, -1.0, -1.0, 0,
true, false, Interface::lowerlim);
static Parameter<AriadneHandler,double> interfaceBeta
("Beta",
"The power in the suppression of radiation from extended dipoles. "
"The original soft suppression model used a sharp cutoff in the "
"transverse-momentum--rapidity space of an emitted gluon. This parameter "
"is used to allow emissions with larger transverse momentum according to "
"\\f$P(p_\\perp^2>p_{\\perp cut}^2="
"\\left(\\frac{p_{\\perp cut}^2}{p_\\perp^2}\\right)^\\beta\\f$. if "
"negative, the sharp cutoff is retained.",
&AriadneHandler::theBeta, 2.0, -1.0, 0,
true, false, Interface::lowerlim);
static RefVector<AriadneHandler,Ariadne5::ReweightBase> interfaceReweighters
("Reweighters",
"A vector of objects implementing reweightings of basic dipole "
"emissions. Each dipole emission will be reweighted.",
&AriadneHandler::theReweighters, -1, true, false, true, false, false);
static Parameter<AriadneHandler,int> interfaceMaxEmissions
("MaxEmissions",
"This number specifies the maximum number of emissions from the "
"cascade. If it is set to zero an unlimited number is allowed. "
"This parameter should only be used for debugging purposes.",
&AriadneHandler::theMaxEmissions, 0, 0, 0,
true, false, Interface::lowerlim);
static Reference<AriadneHandler,AlphaSBase> interfaceInternalAlphaS
("InternalAlphaS",
"An internal AlphaSBase object to be used if "
"<interface>RunningCoupling</interface> is set to InternalRunning. "
"If no such object is given, A O1AlphaS object will be created and "
"assigned in the initialization with <interface>LambdaQCD</interface> "
"used as lambda for five flavours.",
&AriadneHandler::theInternalAlphaS, true, false, true, true, false);
static RefVector<AriadneHandler,EmitterBase> interfaceEmitters
("Emitters",
"The vector of EmittorBase objects responsible for generatong and "
"performing emissions according to the Dipole Cascade Model and its "
"extentions. For each dipole, all Emitters are tried in turn if to see "
"if they are able to model emissions. Only one EmittorBase object of "
"each SubClass will be used (earlier ones will override later ones "
"in the vector).",
&AriadneHandler::theEmitters, -1, true, false, true, false, false);
static RefVector<AriadneHandler,BornCheckerBase> interfaceBornCheckers
("BornCheckers",
" A vector of BornCheckerBase objects which are used in the CKKW-L "
"algorithm to check if a reclustered dipole state corresponds to a "
"reasonable Born-level state (the lowest jet-multiplicity state in a "
"CKKW-L merging). At least one object must be assigned in order for "
"the CKKW-L algorithm to work.",
&AriadneHandler::theBornCheckers, -1, true, false, true, false, false);
static Reference<AriadneHandler,ScaleSetter> interfaceScaleSetter
("ScaleSetter",
"The object responsible for choosing the starting scale of the "
"Ariadne dipole shower.",
&AriadneHandler::theScaleSetter, true, false, true, false, true);
static Reference<AriadneHandler,DISFinder> interfaceDISFinder
("DISFinder",
"The object responsible for identifying DIS-like event.",
&AriadneHandler::theDISFinder, true, false, true, true, false);
static Reference<AriadneHandler,ResonanceFinder> interfaceResonanceFinder
("ResonanceFinder",
"The object responsible for identifying DIS-like event.",
&AriadneHandler::theResonanceFinder, true, false, true, false, false);
static Reference<AriadneHandler,QCDDipoleFinder> interfaceQCDFinder
("QCDFinder",
"The Object responsible for identifying QCD Diploles.",
&AriadneHandler::theQCDFinder, true, false, true, false, false);
static Reference<AriadneHandler,EMDipoleFinder> interfaceEMFinder
("EMFinder",
"The Object responsible for identifying electro-magnetic Diploles. "
"If null, all electro-magnetic radiation will be switched off.",
&AriadneHandler::theEMFinder, true, false, true, true, false);
static Reference<AriadneHandler,ColourResonanceModel>
interfaceColourResonanceModel
("ColourResonanceModel",
"The object responsible for radiating from decay products from "
"coloured resonances.",
&AriadneHandler::theColourResonanceModel, true, false, true, true, false);
static Reference<AriadneHandler,RemnantModel> interfaceRemnantModel
("RemnantModel",
"The object responsible for radiating from remnants.",
&AriadneHandler::theRemnantModel, true, false, true, false, false);
static Reference<AriadneHandler,ConsistencyChecker> interfaceConsistencyChecker
("ConsistencyChecker",
"The object responsible for checking the consistency of a DipoleState.",
&AriadneHandler::consistency, true, false, true, true, false);
interfacePTCut.rank(10);
interfaceLambdaQCD.rank(9);
interfaceNFlav.rank(8);
interfaceSoftMu.rank(7);
interfaceSoftAlpha.rank(6);
interfaceHardAlpha.rank(5);
interfaceBeta.rank(4);
}
diff --git a/Cascade/AriadneHandler.h b/Cascade/AriadneHandler.h
--- a/Cascade/AriadneHandler.h
+++ b/Cascade/AriadneHandler.h
@@ -1,621 +1,629 @@
// -*- C++ -*-
#ifndef Ariadne5_AriadneHandler_H
#define Ariadne5_AriadneHandler_H
//
// This is the declaration of the AriadneHandler class.
//
#include "Ariadne/Config/Ariadne5.h"
#include "ThePEG/Handlers/CascadeHandler.h"
#include "ThePEG/StandardModel/AlphaSBase.h"
#include "AriadneHandler.fh"
#include "History.fh"
#include "ReweightBase.fh"
#include "DipoleState.fh"
#include "EmitterBase.fh"
#include "BornCheckerBase.fh"
#include "ScaleSetter.fh"
#include "DISFinder.fh"
#include "QCDDipoleFinder.fh"
#include "EMDipoleFinder.fh"
#include "ResonanceFinder.fh"
#include "History.h"
#include "ConsistencyChecker.fh"
#include "QCDDipole.fh"
#include "EMDipole.fh"
#include "Ariadne/Cascade/Models/RemnantModel.fh"
#include "Ariadne/Cascade/Models/ColourResonanceModel.fh"
+#include "ThePEG/Analysis/FactoryBase.h"
namespace Ariadne5 {
using namespace ThePEG;
/**
* The AriadneHandler class inherits form the ThePEG::CascadeHandler
* base class and implements the Dipole Cascade Model for partonic
* cascades.
*
* @see \ref AriadneHandlerInterfaces "The interfaces"
* defined for AriadneHandler.
*/
class AriadneHandler: public CascadeHandler {
public:
/**
* Enum different options for running \f$\alpha_s\f$.
*/
enum RunningOption {
externalRunning = -1, /**< Use the \f$\alpha_s\f$ specified in the
current StandardModel object. */
noRunning = 0, /**< Use a fixed \f$\alpha_s\f$ specified by
alpha0(). */
simpleRunning = 1, /**< Use simple leading order running with
\f$\Lambda_{QCD}\f$ given by lambdaQCD() */
internalRunning = 2 /**< Use internal AlphaSBase object in
theInternalAlphaS for the running. */
};
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
AriadneHandler();
/**
* The destructor.
*/
virtual ~AriadneHandler();
//@}
public:
/**
* The main function to be overwritten by sub-classes. It is called
* by handle() after storing some information which is then
* available through simple access functions.
*/
virtual void cascade();
/**
* The AriadneHandler can be used inside the process generation to do
* so-called CKKW reweighting of the hard sub-process. In this case
* this function is called after information about the sub-process is
* made available through the LastXCombInfo base class. Only the
* function belonging to the primary AriadneHandler for the event to
* be generated is called. This default implementation of the function
* simply return one. The current sub-process is mixed together with
* other processes with a multiplicity of outgoing particles between
* \a minMult and \a maxMult.
* @throws CKKWMultiplicityException if the minMult > maxMult or if
* the number of outgoing particle is outside the specified range.
*/
virtual double reweightCKKW(int minMult, int maxMult);
/**
* Check if the particles in the dipole state passes the cuts from
* lastXComb.
*/
bool passCuts(tcDipoleStatePtr state);
public:
/** @name Functions relating to the running coupling. */
//@{
/**
* Strategy for \f$\alpha_S\f$. 0 means constant a constant
* \f$\alpha_S\f$ with a value given by alpha0(). 1 means a leading
* order running \f$\alpha_S\f$ with a \f$\Lambda_{QCD}\f$ given by
* lambdaQCD(). -1 means take whatever is specified in the current
* StandardModelBase object.
*/
inline RunningOption runningCoupling() const {
return theRunningCoupling;
}
/**
* The constant \f$\alpha_S\f$ to be used if runningCoupling() is 0.
*/
inline double alpha0() const {
return theAlpha0;
}
/**
* Return the \f$\alpha_S\f$ to be used for the given \a scale.
*/
double alphaS(Energy2 scale) const;
/**
* The \f$\Lambda_{QCD}\f$ to use in the one loop running
* \f$\alpha_S\f$ if runningCoupling is 1
*/
inline Energy lambdaQCD() const {
return theLambdaQCD;
}
/**
* An internal \f$\alpha_S\f$ object to be used if runningCoupling()
* is internalRunning.
*/
inline Ptr<AlphaSBase>::const_pointer internalAlphaS() const {
return theInternalAlphaS;
};
/**
* Return the flavour thresholds used in calculating \f$\alpha_S\f$.
*/
inline const set<Energy2> & flavourThresholds() const {
return theFlavourThresholds;
}
/**
* The constant \f$\alpha_{EM}\f$ to be used. If zero, use whatever
* is specified in the current StandardModelBase object.
*/
inline double alphaEM0() const {
return theAlphaEM0;
}
/**
* The number of different colour indices available to dipoles.
*/
inline int nCol() const {
return theNCol;
}
/**
* The number of possible flavours in a \f$g\to q\bar{q}\f$ splitting.
*/
inline int nFlav() const {
return theNFlav;
}
//@}
/**
* Check if photon emission are switched on or off.
*/
inline bool photonEmissions() const {
return thePhotonEmissions;
}
public:
/**
* Calculate the starting scale for the given DipoleState.
*/
Energy startingScale(const DipoleState & state) const;
/**
* If DIS-like scattered leptons are found in the given SubProcess,
* return properly setup corresponding hard remnants. The
* remnants will belong to the given DipoleState.
*/
pair<tRemParPtr,tRemParPtr>
findDISLeptons(SubProcess & sub, DipoleState & state) const;
/**
* If findDISLeptons() has found scattered leptons, find also the
* scattered quarks in the SubProcess and return them as hard
* remnants. The remnants will belong to the given DipoleState.
*/
pair<tRemParPtr,tRemParPtr>
findDISQuarks(pair<tRemParPtr,tRemParPtr> leptons,
SubProcess & sub, DipoleState & state) const;
/**
* Return a list of intermediate s-channel resonances in the given
* SubProcess. Resonances which have decayed into further resonances
* will come before their children in the list.
*/
tPVector resonances(SubProcess & sub) const;
/**
* Return the object responsible for radiating from decay products
* from coloured resonances.
*/
tColourResonanceModelPtr colourResonanceModel() const {
return theColourResonanceModel;
}
/**
* Return the object responsible for radiating from remnants.
*/
const RemnantModel & remnantModel() const {
return *theRemnantModel;
}
/**
* Find, create and return the QCD dipoles in the given
* DipoleState.
*/
vector<tQCDPtr> findQCDDipoles(DipoleState & state) const;
/**
* Find, create and return the electro-magnetic dipoles in the given
* DipoleState.
*/
vector<tEMDipPtr> findEMDipoles(DipoleState & state) const;
/**
* Check the given DipoleState for consistency.
*/
bool checkState(DipoleState &, tcEmPtr = tEmPtr());
/**
* The cutoff in invariant transverse momentum for QCD emissions.
*/
inline Energy pTCut() const {
return thePTCut;
}
/**
* The cutoff in invariant transverse momentum for QED emissions.
*/
inline Energy pTCutEM() const {
return thePTCutEM;
}
/**
* The inverse extension of a hadron remnant.
*/
inline Energy softMu() const {
return theSoftMu;
}
/**
* The dimension assumed for the extension of a hadron remnant.
*/
inline double softAlpha() const {
return theSoftAlpha;
}
/**
* The dimension assumed for the extension of a perturbative remnant.
*/
inline double hardAlpha() const {
return theHardAlpha >= 0? theHardAlpha: softAlpha();
}
/**
* The power in the suppression of radiation from extended dipoles.
*/
inline double beta() const {
return theBeta;
}
/**
* Return the vector of available reweighting objects.
*/
inline const vector<DipoleRWPtr> & reweighters() const {
return theReweighters;
}
/**
* The maximum number of allowed emission in the cascade.
*/
inline int maxEmissions() const {
return theMaxEmissions;
}
/**
* Return the vector of available emitter objects.
*/
inline const vector<EmitterPtr> & emitters() const {
return theEmitters;
}
/**
* Check if the given dipole state corresponds to a valid
* Born-level state in a CKKW-L merging and return the scale
* associated with that state. If not return the negative of the
* scale associated with the state.
*/
Energy checkBornState(const DipoleState &) const;
/**
* Check if the given dipole state corresponds to a valid state from
* a matrix element generator in a CKKW-L merging.
*/
bool checkTreeState(const DipoleState &) 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:
/**
* Exception class used if the outgoing multiplicity is out of range
* or the maximum multiplicity is larger than the minimum in
* reweightCKKW.
*/
struct CKKWMultiplicityException: Exception {};
/**
* Exception class used if no reasonable state was reconstructed in
* the CKKW-L algorithm.
*/
struct CKKWBornException: Exception {};
/**
* Exception class used if the dipole state is not self consistant at
* any stage of the cascade.
*/
struct IntegretyException: Exception {};
/**
* Exception class indicating no model could be found when requested.
*/
struct MissingModel: public Exception {};
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object 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();
/**
* Return true if this object needs to be initialized before all
* other objects because it needs to extract PDFs from the event
* file. This version will return true if runningCoupling() returns
* internalRunning and no internalAlphaS() object has been given.
*/
virtual bool preInitialize() const;
//@}
private:
/**
* Function used by the interface.
*/
Energy minPTCut() const;
/**
* Function used by the interface.
*/
Energy maxLambdaQCD() const;
private:
/**
* Strategy for \f$\alpha_S\f$. 0 means constant a constant
* \f$\alpha_S\f$ with a value given by alpha0(). 1 means a leading
* order running \f$\alpha_S\f$ with a \f$\Lambda_{QCD}\f$ given by
* lambdaQCD(). -1 means take whatever is specified in the current
* StandardModelBase object.
*/
RunningOption theRunningCoupling;
/**
* The constant \f$\alpha_S\f$ to be used if runningCoupling() is
* noRunning.
*/
double theAlpha0;
/**
* The \f$\Lambda_{QCD}\f$ to use in the one loop running
* \f$\alpha_S\f$ if runningCoupling is simpleRunning.
*/
Energy theLambdaQCD;
/**
* An internal \f$\alpha_S\f$ object to be used if runningCoupling()
* is internalRunning.
*/
Ptr<AlphaSBase>::pointer theInternalAlphaS;
/**
* Scale factor used to multiply the emission scales in the argument
* of \f$\alpha_S\f$.
*/
double scaleFactor;
/**
* The flavour thresholds used in calculating \f$\alpha_S\f$.
*/
set<Energy2> theFlavourThresholds;
/**
* The cutoff in invariant transverse momentum for QCD emissions.
*/
Energy thePTCut;
/**
* The constant \f$\alpha_{EM}\f$ to be used. If zero, use whatever
* is specified in the current StandardModelBase object.
*/
double theAlphaEM0;
/**
* The cutoff in invariant transverse momentum for QED emissions.
*/
Energy thePTCutEM;
/**
* The number of different colour indices available to dipoles.
*/
int theNCol;
/**
* The number of possible flavours in a \f$g\to q\bar{q}\f$ splitting.
*/
int theNFlav;
/**
* Switch to turn on and off photon emission in the cascade.
*/
bool thePhotonEmissions;
/**
* The inverse extension of a hadron remnant.
*/
Energy theSoftMu;
/**
* The dimension assumed for the extension of a hadron remnant.
*/
double theSoftAlpha;
/**
* The dimension assumed for the extension of a hard remnant.
*/
double theHardAlpha;
/**
* The power in the suppression of radiation from extended dipoles.
*/
double theBeta;
/**
* A list of reweighting objects which may be applied to dipole
* emissions.
*/
vector<DipoleRWPtr> theReweighters;
/**
* The maximum number of emissions allowed in the cascade.
*/
int theMaxEmissions;
/**
* The vector of EmittorBase objects responsible for generating and
* performing emissions according to the Dipole Cascade Model and
* its extentions.
*/
vector<EmitterPtr> theEmitters;
/**
* A vector of BornCheckerBase objects which are used in the CKKW-L
* algorithm.
*/
vector<BornCheckerBasePtr> theBornCheckers;
/**
* The object responsible for choosing the starting scale of the
* shower.
*/
ScaleSetterPtr theScaleSetter;
/**
* The object responsible for identifying DIS-like event.
*/
DISFinderPtr theDISFinder;
/**
* The object responsible for identifying DIS-like event.
*/
ResonanceFinderPtr theResonanceFinder;
/**
* The Object responsible for identifying QCD Diploles.
*/
QCDFinderPtr theQCDFinder;
/**
* The Object responsible for identifying electro-magnetic Diploles.
*/
EMFinderPtr theEMFinder;
/**
* The object responsible for radiating from decay products from
* coloured resonances.
*/
ColourResonanceModelPtr theColourResonanceModel;
/**
* The object responsible for radiating from remnants.
*/
RemnantModelPtr theRemnantModel;
/**
* The object responsible for checking the consistency of a DipoleState.
*/
CCheckPtr consistency;
/**
* If set true, the consistencyh checking of emissions in this
* dipole state is suspended as the initial state was not
* consistent.
*/
bool suspendConsistencyChecks;
/**
* A map containing the states generated by reweight CKKW. The map
* assiciates an XComb to a struct containing the dipole state, the
* the clustering History and a boolean which specifies if the state
* has the maximum parton multiplicity or not.
*/
struct CKKWState {
HistoryPtr history;
LorentzRotation rotation;
bool maxMult;
};
map<tXCPtr, CKKWState> theCKKWMap;
private:
/**
+ * Histograms for debugging purposes.
+ */
+ FactoryBase::tH1DPtr histswing, histall, histglue, histqq;
+
+private:
+
+ /**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
AriadneHandler & operator=(const AriadneHandler &);
};
}
#endif /* Ariadne5_AriadneHandler_H */
diff --git a/Cascade/DipoleState.cc b/Cascade/DipoleState.cc
--- a/Cascade/DipoleState.cc
+++ b/Cascade/DipoleState.cc
@@ -1,670 +1,680 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DipoleState class.
//
#include "DipoleState.h"
#include "Resonance.h"
#include "ResonanceParton.h"
#include "RemnantParton.h"
#include "AriadneHandler.h"
#include "EMDipole.h"
#include "QCDDipole.h"
#include "StateDipole.h"
#include "Junction.h"
#include "PartonTraits.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/EventRecord/Step.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/PDF/PartonBinInstance.h"
#include "ThePEG/Config/algorithm.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Utilities/MaxCmp.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Utilities/DebugItem.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Ariadne5;
DipoleState::DipoleState(tPPair inc): theIncoming(inc), theNe(0) {
dipindx(tcDBPtr());
parindx(tcParPtr());
}
DipoleState::~DipoleState() {}
ClonePtr DipoleState::clone() const {
return new_ptr(*this);
}
void DipoleState::setup(SubProcess & sub) {
set<tPPtr> left(sub.outgoing().begin(), sub.outgoing().end());
// First check if this is a DIS-like event.
pair<tRemParPtr,tRemParPtr> leptons = handler()->findDISLeptons(sub, *this);
pair<tRemParPtr,tRemParPtr> quarks =
handler()->findDISQuarks(leptons, sub, *this);
if ( quarks.first ) {
quarks.first->setupHard(1);
leptons.first->setupHard(1);
theRemnants.first.push_back(quarks.first);
theRemnants.first.push_back(leptons.first);
addHadronicFS(quarks.first);
addHardFS(leptons.first);
left.erase(quarks.first->orig());
left.erase(leptons.first->orig());
}
if ( quarks.second ) {
quarks.second->setupHard(-1);
leptons.second->setupHard(-1);
theRemnants.second.push_back(quarks.second);
theRemnants.second.push_back(leptons.second);
addHadronicFS(quarks.second);
addHardFS(leptons.second);
left.erase(quarks.second->orig());
left.erase(leptons.second->orig());
}
// Now setup the rest of the remnants.
tPBIPtr pb = handler()->lastXComb().partonBinInstance(sub.incoming().first);
while ( pb && pb->incoming() ) {
tRemParPtr rem = create<RemnantParton>();
rem->setup(*pb, 1);
addFS(rem);
theRemnants.first.push_back(rem);
pb = pb->incoming();
}
pb = handler()->lastXComb().partonBinInstance(sub.incoming().second);
while ( pb && pb->incoming() ) {
tRemParPtr rem = create<RemnantParton>();
rem->setup(*pb, -1);
addFS(rem);
theRemnants.second.push_back(rem);
pb = pb->incoming();
}
// After that we should identify all resonances.
tPVector resv = handler()->resonances(sub);
theResonances.resize(resv.size());
for ( int i = 0, N = resv.size(); i < N; ++i ) {
tResPtr res = create<Resonance>();
res->orig(resv[i]);
res->decaySystem(i + 1);
theResonances[i] = res;
}
for ( int i = 0, N = resv.size(); i < N; ++i ) {
tResPtr res = theResonances[i];
if ( resv[i]->parents().size() == 1 ) {
tPVector::iterator pit = find(resv, resv[i]->parents()[0]);
if ( pit != resv.end() )
res->parentResonance(theResonances[pit - resv.begin()]);
}
// Children in the final state should be of type ResonanceParton.
for ( int ic = 0, Nc = resv[i]->children().size(); ic < Nc; ++ic ) {
set<tPPtr>::iterator pit = find(left, resv[i]->children()[ic]);
if ( pit != left.end() ) {
tResParPtr rp = create<ResonanceParton>();
rp->orig(*pit);
rp->setResonance(res);
addHadronicFS(rp);
left.erase(pit);
}
}
}
// What is left is just ordinary partons.
setup(left);
// Save a pointer to the SubProcess
subprocess = &sub;
}
void DipoleState::setup(const set<tPPtr> & out) {
subprocess = tSubProPtr();
// Here we only have ordinary partons, life is simple.
for ( set<tPPtr>::iterator pit = out.begin(); pit != out.end(); ++pit ) {
tParPtr p = create<Parton>();
p->orig(*pit);
addHadronicFS(p);
}
// Now identify all QCD dipoles.
vector<tQCDPtr> dips = handler()->findQCDDipoles(*this);
active.insert(dips.begin(), dips.end());
// And the EM dipoles.
vector<tEMDipPtr> emdips = handler()->findEMDipoles(*this);
active.insert(emdips.begin(), emdips.end());
// Finally we add a StateDipole to be used by global models.
active.insert(create<StateDipole>());
sumTotalMomentum();
}
Energy DipoleState::select(Energy rhomin, Energy rhomax) {
theSelected = EmPtr();
// Energy rho = ZERO;
MaxCmp<Energy> rhosel;
for ( set<tDBPtr>::iterator it = active.begin(); it != active.end(); ++it ) {
pair<set<EmissionGenerator>::iterator, bool> ins =
generators.insert(EmissionGenerator(*it));
const EmissionGenerator & gen = *ins.first;
if ( ins.second ) gen.init();
else gen.reinit();
if ( rhosel(gen.generate(rhomin, rhomax)) && gen.emission ) {
if ( !gen.emission->geno ) gen.emission->geno = theNe + 1;
theSelected = gen.emission;
rhomin = gen.rho();
}
}
for ( BaseSet::iterator it = objects.begin(); it != objects.end(); ++it)
(**it).untouch();
if ( selected() ) emissions.push_back(selected());
untouch();
return rhosel;
}
bool DipoleState::perform() {
static DebugItem printsteps("Ariadne5::PrintSteps", 6);
touch();
if ( !selected() ) return false;
if ( printsteps )
cerr << "Ariadne5::PrintStep: Performing emission " << theNe + 1
<< " of type " << typeid(*selected()).name() << endl;
// Save four momentum before the emission.
LorentzMomentum ptot = totalMomentum();
if ( !selected()->perform() ) {
selected()->dipole->touch();
if ( printsteps )
cerr << "Ariadne5::PrintStep: Emission of type "
<< typeid(*selected()).name() << " failed" << endl;
return false;
}
// Inform other partons that an emission was performed.
for ( set<tParPtr>::const_iterator it = finalState().begin();
it != finalState().end(); ++it )
(**it).notify(*selected());
++theNe;
selected()->emno = theNe;
// Check conservation of four momentum.
Energy2 s = ptot.m2();
ptot -= sumTotalMomentum();
if( (ptot.vect().mag2() + sqr(ptot.e()) )/s > 1.0e-16)
throw MomentumException()
<< "Ariadne5::DipoleState::perform: Emission of type "
<< StringUtils::typeName(typeid(*selected()))
<< " did not conserve four momentum." << Exception::eventerror;
// Do additional consistency checks.
return selected()->failsafe ||
Current<AriadneHandler>()->checkState(*this, selected());
}
void DipoleState::fill(Step & step) const {
set<tParPtr> done;
// Determine the scale to be assigned to the produced particles.
Energy2 scale = sqr(handler()->pTCut());
if ( selected() ) scale = max(scale, sqr(selected()->rho));
// First create all the ThePEG::Particles
for ( set<tParPtr>::const_iterator it = finalState().begin();
it != finalState().end(); ++it )
if ( (**it).produceParticle() && (**it).particle() )
(**it).particle()->scale(scale);
// Then it's time to fix all colour lines.
vector<ColinePtr> colines;
for ( set<tDBPtr>::const_iterator it = active.begin();
it != active.end(); ++it )
if ( ColinePtr cl = (**it).colourLine() ) colines.push_back(cl);
// Now collect all relevant particles (ie. ell except soft
// remnants), andfind all original final state partons have actually
// radiated.
set<tPPtr> radiated;
vector<PPtr> final;
for ( set<tParPtr>::const_iterator it = hardFS().begin();
it != hardFS().end(); ++it ) {
if ( !(**it).orig() ) (**it).getOriginalParents(inserter(radiated));
if ( tPPtr p = (**it).particle() ) final.push_back(p);
}
// Now we can fix all mother-daughter relationships and add new
// particles to the step.
for ( set<tParPtr>::const_iterator it = hardFS().begin();
it != hardFS().end(); ++it )
if ( (**it).orig() && !member(radiated, (**it).orig()) )
step.setCopy((**it).orig(), (**it).particle());
else {
vector<tPPtr> parents;
(**it).getOriginalParents(inserter(parents));
step.addDecayProduct(parents.begin(), parents.end(),
(**it).particle(), false);
}
// Now add all intermediate resonances them to the step and fix up
// the parenthood.
for ( int i = 0, N = resonances().size(); i < N; ++i ) {
final.push_back(resonances()[i]->produceParticle());
step.setCopy(resonances()[i]->orig(), resonances()[i]->particle());
if ( tColinePtr cl = resonances()[i]->orig()->colourLine() )
cl->addColoured(resonances()[i]->particle());
if ( tColinePtr cl = resonances()[i]->orig()->antiColourLine() )
cl->addAntiColoured(resonances()[i]->particle());
}
for ( int i = 0, N = resonances().size(); i < N; ++i ) {
tPPtr orig = resonances()[i]->orig();
for ( int ic = 0, Nc = orig->children().size(); ic < Nc; ++ic ) {
tPPtr child = orig->children()[ic];
if ( child->next() )
step.addDecayProduct(orig->final(), child->final());
else if ( child->children().size() )
step.addDecayProduct(orig->final(), child->children().begin(),
child->children().end(), false);
}
}
if ( !subprocess ) return;
// Now start with the innermost soft remnants and go outwards and
// reextract the incoming particles.
// First skip all hard remnants.
vector<tRemParPtr>::const_iterator rit1 = remnants().first.begin();
while ( rit1 != remnants().first.end() && (**rit1).hard() ) ++rit1;
vector<tRemParPtr>::const_iterator rit2 = remnants().second.begin();
while ( rit2 != remnants().second.end() && (**rit2).hard() ) ++rit2;
PartonExtractor & pex = *(handler()->lastExtractor());
while ( rit1 != remnants().first.end() || rit2 != remnants().second.end() ) {
// If there are no soft remnants on one side oldinc=newinc is
// taken to be the colliding particle and nothing will happen on
// that side.
PPair oldinc = subprocess->incoming();
PPair newinc = subprocess->incoming();
if ( rit1 != remnants().first.end() ) {
oldinc.first = (**rit1).originalExtracted();
newinc.first = (**rit1).extracted();
}
if ( rit2 != remnants().second.end() ) {
oldinc.second = (**rit2).originalExtracted();
newinc.second = (**rit2).extracted();
}
Lorentz5Momentum p1 = newinc.first->momentum();
Lorentz5Momentum p2 = newinc.second->momentum();
// Re-extract the remnants and calculate the boost needed to put
// them on-shell.
PBIPair newbins = pex.newRemnants(oldinc, newinc, &step);
LorentzRotation tot =
pex.boostRemnants(newbins, p1, p2,
newbins.first && newbins.first->incoming(),
newbins.second && newbins.second->incoming());
// Perform the boost also for the final state.
Utilities::transform(final.begin(), final.end(), tot);
// Add the new remnants to the final state.
if ( rit1 != remnants().first.end() ) {
newinc.first->addChild(oldinc.first);
step.addIntermediate(newinc.first);
final.insert(final.end(), newbins.first->remnants().begin(),
newbins.first->remnants().end());
++rit1;
}
if ( rit2 != remnants().second.end() ) {
newinc.second->addChild(oldinc.second);
step.addIntermediate(newinc.second);
final.insert(final.end(), newbins.second->remnants().begin(),
newbins.second->remnants().end());
++rit2;
}
}
}
const Lorentz5Momentum & DipoleState::sumTotalMomentum() {
theTotalMomentum = Utilities::sumMomentum(finalState());
theHardMomentum = Utilities::sumMomentum(hardFS());
theHadronicMomentum = Utilities::sumMomentum(hadronicFS());
theTotalMomentum.rescaleMass();
theHardMomentum.rescaleMass();
theHadronicMomentum.rescaleMass();
return theTotalMomentum;
}
+double DipoleState::lambdaMeasure(Energy2 scale) const {
+ double lam = 0.0;
+ for ( set<tDBPtr>::const_iterator it = activeDipoles().begin();
+ it != activeDipoles().end(); ++it )
+ if ( tQCDPtr d = dynamic_ptr_cast<tQCDPtr>(*it) ) {
+ lam += log(d->sdip()/scale);
+ }
+ return lam;
+}
+
tParPtr DipoleState::create(tcPDPtr pd, tcParPtr parent, bool hfs) {
tParPtr p = create<Parton>();
p->data(pd);
if ( parent ) p->origSystem(parent->system());
if ( hfs ) addHadronicFS(p);
return p;
}
void DipoleState::forgetParton(tParPtr p) {
theHadronicFinalState.erase(p);
theHardFinalState.erase(p);
theFinalState.erase(p);
objects.erase(p);
}
void DipoleState::forgetDipole(tDBPtr d) {
active.erase(d);
objects.erase(d);
}
DipoleStatePtr DipoleState::fullclone() const {
TranslationMap trans;
DipoleStatePtr copy = preclone(trans);
copy->postclone(trans);
return copy;
}
DipoleStatePtr DipoleState::preclone(TranslationMap & trans) const {
CloneSet tocopy(objects.begin(), objects.end());
tocopy.insert(tcDipoleStatePtr(this));
CloneSet tocheck = tocopy;
while ( ! tocheck.empty() ) {
CloneSet additional;
for ( CloneSet::const_iterator it = tocheck.begin();
it != tocheck.end(); ++it )
if ( *it ) (**it).fillReferences(additional);
tocheck.clear();
for ( CloneSet::const_iterator it = additional.begin();
it != additional.end(); ++it )
if ( tocopy.insert(*it).second )
tocheck.insert(*it);
}
tocopy.erase(cDipoleStatePtr(this));
DipoleStatePtr copy = dynamic_ptr_cast<DipoleStatePtr>(clone());
trans[tcDipoleStatePtr(this)] = copy;
for ( CloneSet::const_iterator it = tocopy.begin();
it != tocopy.end(); ++it ) if ( *it ) trans[*it] = (**it).clone();
return copy;
}
void DipoleState::postclone(const TranslationMap & trans) const {
for ( TranslationMap::const_iterator it = trans.map().begin();
it != trans.map().end(); ++it ) it->second->rebind(trans);
}
void DipoleState::fillReferences(CloneSet & cs) const {
CascadeBase::fillReferences(cs);
cs.insert(theSelected);
for ( set<EmissionGenerator>::iterator it = generators.begin();
it != generators.end(); ++it )
if ( it->emission ) cs.insert(it->emission);
}
void DipoleState::rebind(const TranslationMap & trans) {
BaseSet old;
old.swap(objects);
objects.clear();
trans.translate(inserter(objects), old.begin(), old.end());
set<tParPtr> oldfs;
oldfs.swap(theFinalState);
theFinalState.clear();
trans.translate(inserter(theFinalState), oldfs.begin(), oldfs.end());
oldfs.swap(theHardFinalState);
theHardFinalState.clear();
trans.translate(inserter(theHardFinalState), oldfs.begin(), oldfs.end());
oldfs.swap(theHadronicFinalState);
theHadronicFinalState.clear();
trans.translate(inserter(theHadronicFinalState), oldfs.begin(), oldfs.end());
vector<tRemParPtr> orem;
orem.swap(theRemnants.first);
theRemnants.first.clear();
trans.translate(inserter(theRemnants.first), orem.begin(), orem.end());
orem.swap(theRemnants.second);
theRemnants.second.clear();
trans.translate(inserter(theRemnants.second), orem.begin(), orem.end());
vector<tResPtr> ores;
ores.swap(theResonances);
theResonances.clear();
trans.translate(inserter(theResonances), ores.begin(), ores.end());
set<tDBPtr> odip;
odip.swap(active);
active.clear();
trans.translate(inserter(active), odip.begin(), odip.end());
theSelected = trans.translate(theSelected);
set<EmissionGenerator> gold;
gold.swap(generators);
for ( set<EmissionGenerator>::iterator it = gold.begin();
it != gold.end(); ++ it ) {
EmissionGenerator gen(trans.translate(it->dipole));
gen.emission = trans.translate(it->emission);
generators.insert(gen);
}
dipindx.clear();
parindx.clear();
}
void DipoleState::persistentOutput(PersistentOStream & os) const {
os << subprocess << theIncoming << objects << theFinalState << theRemnants
<< theResonances
<< active << ounit(theTotalMomentum, GeV) << ounit(theHardMomentum, GeV)
<< ounit(theHadronicMomentum, GeV) << generators.size()
<< theNe << theSelected;
for ( set<EmissionGenerator>::const_iterator it = generators.begin();
it != generators.end(); ++it ) os << *it;
}
void DipoleState::persistentInput(PersistentIStream & is, int) {
int gsize;
is >> subprocess >> theIncoming >> objects >> theFinalState >> theRemnants
>> theResonances
>> active >> iunit(theTotalMomentum, GeV) >> iunit(theHardMomentum, GeV)
>> iunit(theHadronicMomentum, GeV) >> gsize >>
theNe >> theSelected;
generators.clear();
EmissionGenerator eg;
while ( gsize-- ) {
is >> eg;
generators.insert(eg);
}
}
void DipoleState::touchHadronicState() {
for ( int i = 0, N = theRemnants.first.size(); i < N; ++i )
if ( hadronicFS().find(theRemnants.first[i]) == hadronicFS().end() )
theRemnants.first[i]->touch();
for ( int i = 0, N = theRemnants.second.size(); i < N; ++i )
if ( hadronicFS().find(theRemnants.second[i]) == hadronicFS().end() )
theRemnants.second[i]->touch();
}
void DipoleState::untouchHadronicState() {
for ( int i = 0, N = theRemnants.first.size(); i < N; ++i )
if ( hadronicFS().find(theRemnants.first[i]) == hadronicFS().end() )
theRemnants.first[i]->untouch();
for ( int i = 0, N = theRemnants.second.size(); i < N; ++i )
if ( hadronicFS().find(theRemnants.second[i]) == hadronicFS().end() )
theRemnants.second[i]->untouch();
}
void DipoleState::transformHadronicState(const LorentzRotation & R) {
UtilityBase::transform(theHadronicFinalState, R);
touchHadronicState();
}
// The following static variable is needed for the type description
// system in ThePEG.
DescribeClass<DipoleState,CascadeBase>
describeAriadne5DipoleState("Ariadne5::DipoleState", "libAriadne5.so");
void DipoleState::Init() {}
pair<tQCDPtr,tQCDPtr> DipoleState::StringEnds(tQCDPtr d) {
pair<tQCDPtr,tQCDPtr> ret = make_pair(d, d);
while ( tQCDPtr dn = ret.second->next() ) {
if ( dn == d ) return ret;
ret.second = dn;
}
while ( tQCDPtr dp = ret.first->prev() ) {
if ( dp == d ) return ret;
ret.first = dp;
}
return ret;
}
void DipoleState::debugme() const {
CascadeBase::debugme();
cerr << "DipoleState:" << endl;
LorentzMomentum sum;
set<tParPtr> allpartons = finalState();
set<tDBPtr> dipoles = activeDipoles();
cerr << "> active dipoles:" << endl;
while ( !dipoles.empty() ) {
if ( tQCDPtr d = dynamic_ptr_cast<tQCDPtr>(*dipoles.begin()) ) {
dipoles.erase(d);
pair<tQCDPtr,tQCDPtr> range = StringEnds(d);
tParPtr p1 = range.first->iPart();
tParPtr p2 = range.second->oPart();
if ( p1 == p2 )
cerr << " string loop:" << endl;
else
cerr << " string: " << index(p1)
<< " ... " << index(p2) << endl;
while ( range.first ) {
range.first->iPart()->debug();
cerr << endl;
sum += range.first->iPart()->momentum();
allpartons.erase(range.first->iPart());
range.first->debug();
cerr << endl;
dipoles.erase(range.first);
if ( range.first == range.second )
range.first = tQCDPtr();
else
range.first = range.first->next();
}
if ( p1 != p2 ) {
p2->debug();
cerr << endl;
sum += p2->momentum();
allpartons.erase(p2);
}
} else {
cerr << " other dipole:" << endl;
(**dipoles.begin()).debug();
cerr << endl;
dipoles.erase(dipoles.begin());
}
}
cerr << "> rest of final state:" << endl;
while ( !allpartons.empty() ) {
(**allpartons.begin()).debug();
cerr << endl;
sum += (**allpartons.begin()).momentum();
allpartons.erase(allpartons.begin());
}
cerr << "> sum of momenta:" << setprecision(3)
<< setw(9) << ( abs(sum.x()) > MeV? sum.x(): 0.0*GeV )/GeV
<< setw(9) << ( abs(sum.y()) > MeV? sum.y(): 0.0*GeV )/GeV
<< setw(9) << ( abs(sum.z()) > MeV? sum.z(): 0.0*GeV )/GeV
<< setw(9) << sum.e()/GeV
<< setw(9) << sum.m()/GeV << endl;
if ( selected() ) {
cerr << "> selected emission:" << endl;
selected()->debug();
cerr << "selected dipole: " << index(selected()->cdipole);
cerr << endl;
}
}
void DipoleState::debugEmissions() const {
for ( unsigned i = 0; i < emissions.size(); ++i )
if ( emissions[i] ) emissions[i]->debug();
}
bool DipoleState::checkIntegrity() {
for ( set<tDBPtr>::iterator it = active.begin();
it != active.end(); ++it ) {
if( ! (*it)->checkIntegrity() ) {
return false;
}
}
return true;
}
int DipoleState::index(tcCascadeBasePtr o) const {
dipindx(tcDBPtr());
parindx(tcParPtr());
if ( dynamic_ptr_cast<tcDBPtr>(o) )
return dipindx(dynamic_ptr_cast<tcDBPtr>(o));
if ( dynamic_ptr_cast<tcParPtr>(o) )
return parindx(dynamic_ptr_cast<tcParPtr>(o));
return 0;
}
SaveDipoleState::SaveDipoleState(tDipoleStatePtr state, bool force)
: forced(force) {
if ( force ||
!( state->selected()->failsafe || state->selected()->reversible ) )
backup = state->preclone(trans);
else
backup = state;
}
tDipoleStatePtr SaveDipoleState::revert() {
static DebugItem printsteps("Ariadne5::PrintSteps", 6);
if ( printsteps )
cerr << "Ariadne5::PrintStep: Reverting emission " << backup->theNe + 1
<< " of type " << typeid(*(backup->selected())).name() << endl;
if ( forced ||
!( backup->selected()->failsafe || backup->selected()->reversible ) )
backup->postclone(trans);
else if ( backup->selected()->failsafe )
Throw<Exception>()
<< "The emitter model " << backup->selected()->model->fullName()
<< " has reported that an emission is failsafe, but the emission "
<< "failed anyway. Please contact the author to have this error "
<< "corrected." << Exception::runerror;
else if ( backup->selected()->reversible ) {
// Save four momentum before reverting.
LorentzMomentum ptot = backup->totalMomentum();
backup->selected()->revert();
Energy2 s = ptot.m2();
ptot -= backup->sumTotalMomentum();
if( (ptot.vect().mag2() + sqr(ptot.e()) )/s > 1.0e-16)
throw DipoleState::MomentumException()
<< "Ariadne5::DipoleState::perform: "
<< "Emission did not conserve four momentum."
<< Exception::eventerror;
}
backup->untouch();
return backup;
}
diff --git a/Cascade/DipoleState.h b/Cascade/DipoleState.h
--- a/Cascade/DipoleState.h
+++ b/Cascade/DipoleState.h
@@ -1,579 +1,585 @@
// -*- C++ -*-
#ifndef Ariadne5_DipoleState_H
#define Ariadne5_DipoleState_H
//
// This is the declaration of the DipoleState class.
//
#include "Ariadne/Cascade/CascadeBase.h"
#include "DipoleState.fh"
#include "DipoleBase.h"
#include "Parton.h"
#include "EmissionGenerator.h"
#include "Emission.h"
#include "RemnantParton.h"
#include "Resonance.fh"
#include "ThePEG/Utilities/ObjectIndexer.h"
namespace Ariadne5 {
using namespace ThePEG;
/**
* The DipoleState class describes a complete state of partons and
* dipoles in a Ariadne cascade at a given stage of the process.
*/
class DipoleState: public CascadeBase {
public:
/**
* A vector of CascadeBase pointers.
*/
typedef set<CascadeBasePtr> BaseSet;
/**
* A set of dipole pointers.
*/
typedef list<tDBPtr> DipoleSet;
/**
* The SaveDipoleState helper class is a friend.
*/
friend class SaveDipoleState;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor. Optional argument is the incoming
* particles to the collision.
*/
DipoleState(tPPair inc = tPPair());
/**
* The destructor.
*/
virtual ~DipoleState();
//@}
public:
/**
* Setup this DipoleState from a SubProcess.
*/
void setup(SubProcess &);
/**
* Setup this DipoleState from a simple set of final-state particles
*/
void setup(const set<tPPtr> &);
/**
* Select a dipole. Go through the list of active DipoleBase objects
* and tell them to calculate the scale of the next emission and
* select the Emitter with largest transverse momentum. Subsequent
* calls to selected() will return the selected dipole.
*
* @param rhomin the minimum allowed evolution scale.
*
* @param pt2max the maximum allowed evolution scale.
*
* @return the generated evolution scale. If less that rhomin, no
* emission was generated and subsequent calls to selected() will
* return null.
*/
Energy select(Energy rhomin, Energy rhomax);
/**
* Perform an emission.
*
* @return false if no emission was possible.
*/
bool perform();
/**
* After the cascade is done, fill the produced partons in the
* supplied Step.
*/
void fill(Step &) const;
/** @name Simple access functions. */
//@{
/**
* The next Emission
*/
inline tEmPtr selected() const {
return theSelected;
}
/**
* Set the next Emission.
*/
inline void selected(const tEmPtr & x) {
theSelected = x;
}
/**
* The number of emissions made from this state so far.
*/
inline int nEmissions() const {
return theNe;
}
/**
* Return the colliding particles if present.
*/
tPPair incoming() const {
return theIncoming;
}
/**
* Return the final-state partons. This includes all final-state
* particles, also remnants.
*/
const set<tParPtr> & finalState() const {
return theFinalState;
}
/**
* Add a parton to the final state.
*/
void addFS(tParPtr p) {
theFinalState.insert(p);
}
/**
* Add a parton to the hard final state.
*/
void addHardFS(tParPtr p) {
theHardFinalState.insert(p);
addFS(p);
}
/**
* Add a parton to the hadronic final state.
*/
void addHadronicFS(tParPtr p) {
theHadronicFinalState.insert(p);
addHardFS(p);
}
/**
* Remove the given parton from the final state and forget it ever
* existed.
*/
void forgetParton(tParPtr p);
/**
* Remove the given dipole from the active dipoles and forget it ever
* existed.
*/
void forgetDipole(tDBPtr d);
/**
* Return the partons in the hard final state. Returns finalState()
* excluding all soft remnants.
*/
const set<tParPtr> & hardFS() const {
return theHardFinalState;
}
/**
* Return the partons in the hadronic final state excluding all soft
* and non-coloured remnants.
*/
const set<tParPtr> & hadronicFS() const {
return theHadronicFinalState;
}
/**
* Return the remnants. Are a subset of theFinalState
*/
const pair< vector<tRemParPtr>, vector<tRemParPtr> > & remnants() const {
return theRemnants;
}
/**
* Return the set of active dipoles.
*/
const set<tDBPtr> & activeDipoles() const {
return active;
}
/**
* Return massive resonances.
*/
const vector<tResPtr> & resonances() const {
return theResonances;
}
//@}
public:
/** @name Functions relating to the book-keeping of included objects. */
//@{
/**
* Create a new object to be included in this DipoleState.
*/
template <typename Class>
inline typename Ptr<Class>::pointer create() {
typedef typename Ptr<Class>::pointer PTR;
PTR obj = ptr_new<PTR>();
objects.insert(objects.end(), obj);
obj->state(this);
if ( tDBPtr d = dynamic_ptr_cast<tDBPtr>(obj) ) {
active.insert(d);
dipindx(d);
}
if ( tParPtr p = dynamic_ptr_cast<tParPtr>(obj) ) parindx(p);
return obj;
}
/**
* Create a Parton of the given type \a pd. Optionally inherit the
* system properties of the \a parent parton and insert in the
* hadronic final state if \a hfs is true.
*/
tParPtr create(tcPDPtr pd, tcParPtr parent = tcParPtr(), bool hfs = true);
/**
* Remove a Parton from the final state. (The parton will still
* exist in the DipoleState).
*/
void remove(tParPtr p) {
theFinalState.erase(p);
theHardFinalState.erase(p);
theHadronicFinalState.erase(p);
}
/**
* Clone this DipoleState. Also cloning all included objects and fixing up
* their inter-dependence.
*/
DipoleStatePtr fullclone() const;
/**
* Clone this DipoleState. Also clone all included objects but do
* not fix up their inter-dependence, ie. the cloned objects in the
* new DipoleState may still point to objects in the this. Before
* the new DipoleState can be used, its postclone() function must
* be called with the \a trans object as argument. The \a trans
* object will contain the translation map between the objects in
* the old and new DipoleState.
*/
DipoleStatePtr preclone(TranslationMap & trans) const;
/**
* Fix up the inter-dependence among the included objects in this
* DipoleState which was produced by the preclone() function. The
* same \a trans object which was used in the preclone() call must
* be given here as argument.
*/
void postclone(const TranslationMap & trans) const;
//@}
public:
/**
* Calculate the total momentum of the dipole state.
*/
const Lorentz5Momentum & sumTotalMomentum();
/**
* Return the total momentum of the dipole state. This is
* recalculated after each emission.
*/
const Lorentz5Momentum & totalMomentum() const {
return theTotalMomentum;
}
/**
* Calculate the total momentum of the hardFinalState().
*/
const Lorentz5Momentum & sumHardMomentum();
/**
* Return the total momentum of the hardFS(). This is
* recalculated after each emission.
*/
const Lorentz5Momentum & hardMomentum() const {
return theHardMomentum;
}
/**
* Return the total momentum of the hadronicFS(). This is
* recalculated after each emission.
*/
const Lorentz5Momentum & hadronicMomentum() const {
return theHadronicMomentum;
}
/**
+ * Return the total sting length in terms of the lambda measure
+ * using the given scale.
+ */
+ double lambdaMeasure(Energy2 scale = 1.0*GeV2) const;
+
+ /**
* Transform the hadronicFS().
*/
void transformHadronicState(const LorentzRotation & R);
/**
* Touch all remnants to indicate that the hadronicFS() has changed
* collectively.
*/
void touchHadronicState();
/**
* Untouch all remnants to indicate that the hadronicFS() has
* reverted collectively.
*/
void untouchHadronicState();
/**
* Print out debugging information on std::cerr.
*/
virtual void debugme() const;
/**
* Print out more debugging information on std::cerr.
*/
void debugEmissions() const;
/**
* Check integrity of all dipoles. Return false if error is
* found.
*/
bool checkIntegrity();
/**
* Create and return an index for the given object. Only used in for
* debugging.
*/
int index(tcCascadeBasePtr o) const;
protected:
/**
* Trace a dipole to find (a piece of) string.
*/
static pair<tQCDPtr,tQCDPtr> StringEnds(tQCDPtr);
protected:
/** @name The virtual functions to be overridden in sub-classes. */
//@{
/**
* Return a simple clone of this object. Should be implemented as
* <code>return new_ptr(*this);</code> by a derived class.
*/
virtual ClonePtr clone() const;
/**
* Fill the provided set with all pointers to CloneBase objects used
* in this object.
*/
virtual void fillReferences(CloneSet &) const;
/**
* Rebind pointers to other CloneBase objects. Called after a number
* of interconnected CloneBase objects have been cloned, so that
* the cloned objects will refer to the cloned copies afterwards.
*
* @param trans a TranslationMap relating the original objects to
* their respective clones.
*/
virtual void rebind(const TranslationMap & trans);
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
private:
/**
* The SubProcess from which this state was created (may be null).
*/
tSubProPtr subprocess;
/**
* The colliding particles in the SubProcess or 0 if none was
* present.
*/
tPPair theIncoming;
/**
* The set of all diferent objects belonging to this state.
*/
BaseSet objects;
/**
* The final state partons.
*/
set<tParPtr> theFinalState;
/**
* The hard final state partons, excluding all soft remnants.
*/
set<tParPtr> theHardFinalState;
/**
* The hadronic final state, excluding all soft and non-coloured
* remnants.
*/
set<tParPtr> theHadronicFinalState;
/**
* The remnants. Are a subset of theFinalState
*/
pair< vector<tRemParPtr>, vector<tRemParPtr> > theRemnants;
/**
* The massive resonances.
*/
vector<tResPtr> theResonances;
/**
* The total momentum of this state.
*/
Lorentz5Momentum theTotalMomentum;
/**
* The total momentum of the hardFS().
*/
Lorentz5Momentum theHardMomentum;
/**
* The total momentum of the hadronicFS().
*/
Lorentz5Momentum theHadronicMomentum;
/**
* The active Dipoles.
*/
set<tDBPtr> active;
/**
* The set of EmissionGenerators corresponding to the dipoles
* included in this state.
*/
set<EmissionGenerator> generators;
/**
* The number of emissions made in this state.
*/
int theNe;
/**
* The Emitter selected to perform the next emission.
*/
EmPtr theSelected;
/**
* The list of emissions treated in this DipoleState.
*/
vector<EmPtr> emissions;
/**
* Keep track of indices for debugging purposes.
*/
mutable ObjectIndexer<int,const DipoleBase> dipindx;
/**
* Keep track of indices for debugging purposes.
*/
mutable ObjectIndexer<int,const Parton> parindx;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DipoleState & operator=(const DipoleState &);
public:
/** @cond EXCEPTIONCLASSES */
/** Exception class used by DipoleState if four momentum is
* not conserved.
*/
class MomentumException: public Exception {};
/** Exception class used by DipoleState if the SubProcess was
* querried when none was available.
*/
class SubProcessException: public Exception {};
/** @endcond */
};
/**
* Helper class to be able to revert changes made to a DipoleState.
*/
class SaveDipoleState {
public:
/**
* The constructor taking a pointer to the \a state to be saved. If
* necessary (or if \a force is true) the state will be cloned.
*/
SaveDipoleState(tDipoleStatePtr state, bool force = false);
/**
* Return the translation map relating the objects in the original
* state with the ones in the backup state.
*/
const DipoleState::TranslationMap & translationMap() const {
return trans;
}
/**
* Return a copy of the original state.
*/
tDipoleStatePtr revert();
private:
/**
* The backup state.
*/
DipoleStatePtr backup;
/**
* The translation map relating the objects in the original state
* with the ones in the backup state.
*/
DipoleState::TranslationMap trans;
/**
* Flag to indicate that precloning was forced.
*/
bool forced;
};
}
#endif /* Ariadne5_DipoleState_H */
diff --git a/Cascade/Models/DipoleSwinger.cc b/Cascade/Models/DipoleSwinger.cc
--- a/Cascade/Models/DipoleSwinger.cc
+++ b/Cascade/Models/DipoleSwinger.cc
@@ -1,300 +1,347 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DipoleSwinger class.
//
#include "FSGluonEmission.h"
#include "DipoleSwinger.h"
#include "Ariadne/Cascade/QCDDipole.h"
#include "Ariadne/Cascade/StateDipole.h"
#include "Ariadne/Cascade/DipoleState.h"
#include "Ariadne/Cascade/AriadneHandler.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/DebugItem.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
+#include "ThePEG/Interface/Command.h"
+#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Repository/CurrentGenerator.h"
+#include "ThePEG/Utilities/EnumIO.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "../../DIPSY/CPUTimer.h"
using namespace Ariadne5;
-DipoleSwinger::DipoleSwinger(): lambda(1.0) {}
+DipoleSwinger::DipoleSwinger(): lambda(1.0), Rmax(3.5*InvGeV*hbarc), linear(false) {}
DipoleSwinger::~DipoleSwinger() {}
IBPtr DipoleSwinger::clone() const {
return new_ptr(*this);
}
IBPtr DipoleSwinger::fullclone() const {
return new_ptr(*this);
}
void DipoleSwinger::TauRatio::plot(Time dt1, Time dt2, int N) const {
double rmax = maximum(dt1, dt2);
for ( int i = 0; i < N; ++i )
cerr << rat(dt1 + (i + 0.5)*(dt2 - dt1)/double(N))/rmax << endl;
}
-DipoleSwinger::TauRatio::TauRatio(const QCDDipole & d1, const QCDDipole & d2)
- : d1tau2(d1.iPart(), d1.oPart()),
- d2tau2(d2.iPart(), d2.oPart()),
- r1tau2(d1.iPart(), d2.oPart()),
- r2tau2(d2.iPart(), d1.oPart()) {}
+DipoleSwinger::TauRatio::TauRatio(const QCDDipole & d1, const QCDDipole & d2, Length Rmaxin)
+ : d1tau2(d1.iPart(), d1.oPart()), d2tau2(d2.iPart(), d2.oPart()),
+ r1tau2(d1.iPart(), d2.oPart()), r2tau2(d2.iPart(), d1.oPart()), Rmax(Rmaxin){}
double DipoleSwinger::TauRatio::rat(Time dt) const {
- return d1tau2(dt)*d2tau2(dt)/(r1tau2(dt)*r2tau2(dt));
+ if ( Rmax <= ZERO )
+ return d1tau2(dt)*d2tau2(dt)/(r1tau2(dt)*r2tau2(dt));
+ else
+ return sqr((exp(sqrt(d1tau2(dt))/Rmax) - 1.0)*(exp(sqrt(d2tau2(dt))/Rmax) - 1.0)/
+ ((exp(sqrt(r1tau2(dt))/Rmax) - 1.0)*(exp(sqrt(r2tau2(dt))/Rmax) - 1.0)));
}
double DipoleSwinger::TauRatio::maximum(Time dt1, Time dt2) const {
// return d1tau2.minmax(dt1, dt2).second*d2tau2.minmax(dt1, dt2).second/
// (r1tau2.minmax(dt1, dt2).first*r2tau2.minmax(dt1, dt2).first);
double rmax = max(rat(dt1), rat(dt2));
double rmax2 = max(rat(r1tau2.extreme(dt1, dt2)), rat(r2tau2.extreme(dt1, dt2)));
if ( rmax2 > rmax ) rmax = rmax2;
return 1.1*rmax;
}
bool DipoleSwinger::canHandle(const DipoleBase & e) const {
tcStateDipPtr d = dynamic_ptr_cast<tcStateDipPtr>(&e);
if ( !d ) return false;
return true;
}
bool DipoleSwinger::overrides(const EmitterBase &, DipoleBase &) const {
return false;
}
bool DipoleSwinger::touched(const DipoleBase & dipole) const {
return dipole.state()->touched() || dipole.touched() || dipole.state() != lastState;
}
EmPtr DipoleSwinger::
generate(const DipoleBase & dipole, Energy rhomin, Energy rhomax) const {
static DebugItem logme("Ariadne5::DipoleSwinger", 6);
static CPUClock cpuclock("Ariadne5::DipoleSwinger::generate");
static CPUClock cpuclock1("Ariadne5::DipoleSwinger::generate1");
static CPUClock cpuclock2("Ariadne5::DipoleSwinger::generate2");
static CPUClock cpuclock3("Ariadne5::DipoleSwinger::generate3");
CPUTimer timer(cpuclock);
const set<tDBPtr> & activeset = dipole.state()->activeDipoles();
map<ColourIndex, vector<tQCDPtr> > cactive;
map<ColourIndex, vector<tQCDPtr> > ctouched;
for ( set<tDBPtr>::const_iterator it = activeset.begin();
it != activeset.end(); ++it )
if ( tQCDPtr dp = dynamic_ptr_cast<tQCDPtr>(*it) ) {
cactive[dp->colourIndex()].push_back(dp);
if ( dp->touched() || dp->iPart()->touched() || dp->oPart()->touched() )
ctouched[dp->colourIndex()].push_back(dp);
}
CPUTimer timer1(cpuclock1);
Time dtmax = hbarc/rhomin;
Time dtmin = hbarc/rhomax;
+ Time t0 = hbarc/Current<AriadneHandler>()->pTCut();
bool redoFull = dipole.state() != lastState;
lastState = dipole.state();
if ( redoFull ) cache.clear();
if ( redoFull || dipole.touched() ) lastSelected = DipSwPtr();
if ( lastSelected && ( lastSelected->dipoles.first->touched() ||
lastSelected->dipoles.second->touched() ) )
lastSelected = DipSwPtr();
if ( logme && redoFull )
generator()->log() << "DipoleSwinger: Starting new event." << endl;
if ( logme ) generator()->log()
<< "DipoleSwinger: Starting new swing search." << endl;
for ( map<ColourIndex, vector<tQCDPtr> >::iterator ciit = cactive.begin();
ciit != cactive.end(); ++ciit ) {
const vector<tQCDPtr> & active = ciit->second;
if ( active.empty() ) continue;
const vector<tQCDPtr> & touched = ctouched[active[0]->colourIndex()];
int ti1 = touched.empty()? -1: 0;
for ( int i1 = 0, N = active.size(); i1 < N; ++i1 ) {
CPUTimer timer2(cpuclock2);
if ( ti1 >= 0 && touched[ti1] == active[i1] ) ++ti1;
QCDDipole & d1 = *active[i1];
bool redo = redoFull || d1.touched() ||
d1.iPart()->touched() || d1.oPart()->touched();
pair<CacheMap::iterator,bool> cit = cache.insert(make_pair(&d1, DipSwPtr()));
if ( redo ) cit.first->second = DipSwPtr();
DipSwPtr sel = cit.first->second;
if ( sel && ( sel->dipoles.second->touched() ||
sel->dipoles.second->iPart()->touched() ||
sel->dipoles.second->oPart()->touched() ) ) {
redo = true;
sel = cit.first->second = DipSwPtr();
}
const vector<tQCDPtr> * secondp = &active;
int i2 = i1 + 1;
if ( !redo ) {
if ( ti1 < 0 ) continue;
secondp = &touched;
i2 = ti1;
}
const vector<tQCDPtr> & second = *secondp;
CPUTimer timer3(cpuclock3);
for ( int N2 = second.size(); i2 < N2; ++i2 ) {
QCDDipole & d2 = *second[i2];
if ( d1.iPart() == d2.oPart() || d2.iPart() == d1.oPart() ) {
cerr << "Bullocks! Try to swing a gluon into colour singlet!" << endl;
}
- TauRatio tauRatio(d1, d2);
+ TauRatio tauRatio(d1, d2, Rmax);
double ratmax = tauRatio.maximum(dtmin, dtmax);
double Cinv = -1.0/(lambda*ratmax);
Time dt = dtmin;
Time dtcut = dtmax;
if ( sel ) dtcut = hbarc/sel->rho;
do {
- dt *= pow(UseRandom::rnd(), Cinv);
+ if ( linear )
+ dt += log(UseRandom::rnd())*Cinv*t0;
+ else
+ dt *= pow(UseRandom::rnd(), Cinv);
if ( dt >= dtcut ) break;
if ( tauRatio(dt) > ratmax ) {
cerr << "Ooops failed to overestimat swing probability: "
<< tauRatio(dt)/ratmax << endl;
}
} while ( dt < dtcut && tauRatio(dt) < UseRandom::rnd()*ratmax );
if ( dt < dtcut ) {
if ( sel ) sel->setup(d1, d2);
else sel = new_ptr(DipoleSwing(*this, dipole, d1, d2));
sel->rho = hbarc/dt;
}
}
cit.first->second = sel;
if ( sel && ( !lastSelected || sel->rho > lastSelected->rho ) ) lastSelected = sel;
}
}
return lastSelected;
}
void DipoleSwinger::swing(tQCDPtr d1, tQCDPtr d2) {
// Swing the o-partons.
tParPtr tmp = d1->oPart();
d1->oPart(d2->oPart());
d2->oPart(tmp);
// The new o-partons inhreits the previous colour line.
tColinePtr cltmp = d1->oPart()->origICol();
d1->oPart()->origICol(d2->oPart()->origICol());
d2->oPart()->origICol(cltmp);
// Fix pointers to neighboring dipoles.
tQCDPtr dtmp = d1->next();
d1->next(d2->next());
d2->next(dtmp);
if ( d1->next() ) d1->next()->prev(d1);
if ( d2->next() ) d2->next()->prev(d2);
}
bool DipoleSwinger::
perform(const Emission & emission) const {
static DebugItem stat("Ariadne5::DipoleSwinger::Statistics", 6);
static double sum = 0.0;
static CPUClock cpuclock("Ariadne5::DipoleSwinger::perform");
CPUTimer timer(cpuclock);
const DipoleSwing & e = dynamic_cast<const DipoleSwing &>(emission);
if ( e.dipoles.first->colourIndex() != e.dipoles.second->colourIndex() )
cerr << "Ooops, tried to swing dipoles of diferent colours!" << endl;
double lam0 = 0.0;
Energy2 s12 = ZERO;
Energy2 s34 = ZERO;
if ( stat ) {
lam0 = log(e.dipoles.first->sdip()/GeV2) + log(e.dipoles.second->sdip()/GeV2);
- TauRatio tauRatio(*e.dipoles.first, *e.dipoles.second);
+ TauRatio tauRatio(*e.dipoles.first, *e.dipoles.second, Rmax);
cerr << "ratio " << tauRatio(hbarc/e.rho);
s12 = e.dipoles.first->sdip();
s34 = e.dipoles.second->sdip();
}
swing(e.dipoles.first, e.dipoles.second);
if ( stat ) {
double lam = log(e.dipoles.first->sdip()/GeV2) + log(e.dipoles.second->sdip()/GeV2);
cerr << " (" << s12*s34/(e.dipoles.first->sdip()*e.dipoles.second->sdip())
<< "), swing diff " << lam - lam0
<< " (total " << (sum += lam - lam0) << ")" << endl;
}
e.dipoles.first->touch();
e.dipoles.second->touch();
e.dipole->touch();
return true;
}
void DipoleSwinger::revert(const Emission & emission) const {
static CPUClock cpuclock("Ariadne5::DipoleSwinger::revert");
CPUTimer timer(cpuclock);
const DipoleSwing & e = dynamic_cast<const DipoleSwing &>(emission);
swing(e.dipoles.first, e.dipoles.second);
e.dipoles.second->untouch();
e.dipole->untouch();
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void DipoleSwinger::persistentOutput(PersistentOStream & os) const {
- os << lambda << cache << lastState << lastSelected;
+ os << lambda << ounit(Rmax, femtometer) << oenum(linear)
+ << cache << lastState << lastSelected;
}
void DipoleSwinger::persistentInput(PersistentIStream & is, int) {
- is >> lambda >> cache >> lastState >> lastSelected;
+ is >> lambda >> iunit(Rmax, femtometer) >> ienum(linear)
+ >> cache >> lastState >> lastSelected;
+}
+
+string DipoleSwinger::setRmax(string cmd) {
+ istringstream is(cmd);
+ double x;
+ is >> x;
+ Rmax = x*InvGeV*hbarc;
+ return "";
}
DescribeClass<DipoleSwinger,EmitterBase>
describeAriadne5DipoleSwinger("Ariadne5::DipoleSwinger", "libAriadne5.so");
void DipoleSwinger::Init() {
static ClassDocumentation<DipoleSwinger> documentation
("The DipoleSwinger class implements the final-state swing method"
"for colour reconnections.");
static Parameter<DipoleSwinger,double> interfaceLambda
("Lambda",
"The frequency of the swings.",
&DipoleSwinger::lambda, 1.0, 1.0, 0.0, 0.0,
true, false, Interface::lowerlim);
+
+ static Parameter<DipoleSwinger,Length> interfaceRmax
+ ("Rmax",
+ "The typical hadronic size used to regularize large dipoles in the swing.",
+ &DipoleSwinger::Rmax, femtometer, 3.5*InvGeV*hbarc, 0.0*femtometer, 0*femtometer,
+ true, false, Interface::lowerlim);
+
+ static Command<DipoleSwinger> interfaceSetRmax
+ ("SetRmax",
+ "Set the value of <interface>Rmax</interface> with a value given in units of inverse GeV.",
+ &DipoleSwinger::setRmax, true);
+
+
+ static Switch<DipoleSwinger,bool> interfaceLinear
+ ("Linear",
+ "Normally we use logarithmic evolution in time. With this switch we can instead use linear evolution.",
+ &DipoleSwinger::linear, false, true, false);
+ static SwitchOption interfaceLinearLogarithmic
+ (interfaceLinear,
+ "Logarithmic",
+ "Use logarithmic evolution in time.",
+ false);
+ static SwitchOption interfaceLinearLinear
+ (interfaceLinear,
+ "Linear",
+ "Use linear evolution in time.",
+ true);
+
}
diff --git a/Cascade/Models/DipoleSwinger.h b/Cascade/Models/DipoleSwinger.h
--- a/Cascade/Models/DipoleSwinger.h
+++ b/Cascade/Models/DipoleSwinger.h
@@ -1,311 +1,331 @@
// -*- C++ -*-
#ifndef ARIADNE5_DipoleSwinger_H
#define ARIADNE5_DipoleSwinger_H
//
// This is the declaration of the DipoleSwinger class.
//
#include "Ariadne/Cascade/EmitterBase.h"
#include "DipoleSwing.h"
#include "Ariadne/Cascade/QCDDipole.h"
namespace Ariadne5 {
using namespace ThePEG;
/**
* The DipoleSwinger class implements the final-state swing method for
* colour reconections.
*
* @see \ref DipoleSwingerInterfaces "The interfaces"
* defined for DipoleSwinger.
*/
class DipoleSwinger: public EmitterBase {
public:
/**
* Convenient typedef.
*/
ThePEG_DECLARE_POINTERS(Ariadne5::DipoleSwing,DipSwPtr);
/**
* Another convenient typedef.
*/
typedef map<tcQCDPtr,DipSwPtr> CacheMap;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
DipoleSwinger();
/**
* The destructor.
*/
virtual ~DipoleSwinger();
//@}
/** @name Virtual functions to be overridden by sub-classes. */
//@{
/**
* Return true if and only if this model can handle the given
* Emitter.
*/
virtual bool canHandle(const DipoleBase &) const;
/**
* If the given EmissionModel overlaps with this model for the given
* Emitter, return true if this model should take precedence. Must
* only be called for Emitters for which canHandle() is true.
*/
virtual bool overrides(const EmitterBase &, DipoleBase &) const;
/**
* Check if objects related to the given \a dipole have been touched
* in a way such that emissions must be regenerated.
*/
virtual bool touched(const DipoleBase & dipole) const;
/**
* Generate the a phase space point for an emission corresponding to
* this model. Must only be called for Emitters for which
* canHandle() is true.
*/
virtual EmPtr generate(const DipoleBase &, Energy rhomin, Energy rhomax) const;
/**
* Perform an emission previously generated for this Emitter. Must
* only be called for Emitters for which canHandle() is true.
* @return true if the emission was successful
*/
virtual bool perform(const Emission &) const;
/**
* Reverse a previously performed emission. Sub-classes which has
* signalled that they can revert an emission but fails to do so,
* must throw a Exception::runerror.
*/
virtual void revert(const Emission & emission) const;
//@}
/**
* Swing the given dipoles.
*/
static void swing(tQCDPtr d1, tQCDPtr d2);
public:
/**
* Internal helper class to keep track of distnaces between partons.
*/
struct DeltaTau2 {
/**
* Constructor taked two partons and calculates the variables
* needed to get the distance as a function of time.
*/
DeltaTau2(tcParPtr p1, tcParPtr p2) {
LorentzDistance dx = p1->vertex() - p2->vertex();
dx2 = dx.m2();
LorentzVector<double> dp =
p1->momentum()/p1->momentum().e() - p2->momentum()/p2->momentum().e();
dp2 = dp.m2();
dpdotdx = dp*dx;
}
/**
* Return distance between partons after a given time \a dt.
*/
Area dTau2(Time dt) const {
return - dx2 - dt*2.0*dpdotdx - sqr(dt)*dp2;
}
/**
* Function version of dTau2
*/
Area operator()(Time dt) const {
return dTau2(dt);
}
/**
* If the distance between the partons may become zero in the
* given time interval, return the corresponding times. If there
* are no zeros, time=0 will be returned.
*/
pair<Length,Length> zeros(Time dt1, Time dt2) const {
pair<Length,Length> ret;
Area squarg = sqr(dpdotdx) - dx2*dp2;
if ( squarg < ZERO ) return ret;
Time t0 = - (dpdotdx + sqrt(squarg))/dp2;
if ( dt1 <= t0 && t0 <= dt2 ) ret.first = t0;
t0 = - (dpdotdx - sqrt(squarg))/dp2;
if ( dt1 <= t0 && t0 <= dt2 ) ret.second = t0;
return ret;
}
/**
* Return the minimum and maximum in the given interval.
*/
pair<Area,Area> minmax(Time dt1, Time dt2) const {
pair<Area,Area> ret(dTau2(dt1),dTau2(dt2));
if ( ret.first > ret.second ) swap(ret.first, ret.second);
Time dtext = -dpdotdx/dp2;
if ( dtext >= dt2 || dt1 >= dtext ) return ret;
Area dtau2ext = dTau2(dtext);
if ( dtau2ext > ret.second ) ret.second = dtau2ext;
if ( dtau2ext < ret.first ) ret.first = dtau2ext;
return ret;
}
/**
* Return the extreme point if in the given interval.
*/
Time extreme(Time dt1, Time dt2) const {
return min(max(dt1, -dpdotdx/dp2), dt2);
}
/**
* The initial distance between the partons.
*/
Area dx2;
/**
* The distance between the momentum vectors, scaled by their energies.
*/
double dp2;
/**
* The scalar product between the difference between momentum
* vectors (scaled with their energies) and the difference between
* the initial positions.
*/
Length dpdotdx;
};
/**
* Internal helper class to keep track of ratio of distances between
* original and reconected dipoles.
*/
struct TauRatio {
/**
* The constructor takes the original dipole pair as argument.
*/
- TauRatio(const QCDDipole & d1, const QCDDipole & d2);
+ TauRatio(const QCDDipole & d1, const QCDDipole & d2, Length Rmaxin);
/**
* Return the swing ratio at the given time.
*/
double rat(Time dt) const;
/**
* Return the swing ratio at the given time.
*/
double operator()(Time dt) const {
return rat(dt);
}
/**
* Estimate a maximum value of the swing ratio.
*/
double maximum(Time dt1, Time dt2) const;
void plot(Time dt1, Time dt2, int N = 20) const;
/**
- * The object for calculating distances between the original and
+ * The objects for calculating distances between the original and
* reconnected dipoles.
*/
DeltaTau2 d1tau2, d2tau2, r1tau2, r2tau2;
+ /**
+ * Hadronic size to regularize large dipoles.
+ */
+ Length Rmax;
+
};
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 only parameter giving the rate of swinging
+ * The parameter giving the rate of swinging
*/
double lambda;
/**
+ * The confinement scale used to prevent swinging to too large dipoles.
+ */
+ Length Rmax;
+
+ /**
+ * Should we use linear or logarithmic evolution in time?
+ */
+ bool linear;
+
+ /**
* Previously calculated cached swings.
*/
mutable CacheMap cache;
/**
* A pointer to the DipoleState last treated.
*/
mutable tDipoleStatePtr lastState;
/**
* A pointer to the last selected generated swing.
*/
mutable DipSwPtr lastSelected;
private:
/**
+ * Helper function used by the interface.
+ */
+ string setRmax(string);
+
+ /**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DipoleSwinger & operator=(const DipoleSwinger &);
};
}
#endif /* ARIADNE5_DipoleSwinger_H */
diff --git a/DIPSY/95.cc b/DIPSY/95.cc
new file mode 100644
--- /dev/null
+++ b/DIPSY/95.cc
@@ -0,0 +1,60 @@
+#include <iostream>
+#include <map>
+#include <cmath>
+
+using namespace std;
+
+int main() {
+
+ const double dummy = sqrt(dummy);
+
+ multimap<double,double> bdist;
+ multimap<double,double> fdist;
+ double w = 0.0;
+ double b = 0.0;
+ double n = 0.0;
+ double f = 0.0;
+ double m = 0.0;
+ double sumw = 0.0;
+ double sumnw = 0.0;
+ double sumnw2 = 0.0;
+
+ while ( cin >> w >> b >> n >> f) {
+ sumw += w;
+ sumnw += n*w;
+ sumnw2 += n*n*w;
+ bdist.insert(make_pair(b, w));
+ fdist.insert(make_pair(f, w));
+ }
+ cerr << "read " << bdist.size() << " lines" << endl;
+ cerr << "sumw " << sumw << endl;
+ cerr << "<n> " << sumnw/sumw
+ << " +- " << sqrt((sumnw2/sumw - sumnw*sumnw/(sumw*sumw))/bdist.size()) << endl;
+ double w5 = 0.05*sumw;
+
+ multimap<double,double>::iterator next = bdist.begin();
+ multimap<double,double>::iterator prev = next;
+ while ( ++next != bdist.end() ) {
+ next->second += prev->second;
+ if ( next->second > w5&& prev->second <= w5 ) {
+ double b5 = prev->first + (w5 - prev->second)*
+ (next->first - prev->first)/(next->second - prev->second);
+ cout << "5% centrality in b < " << b5 << endl;
+ break;
+ }
+ prev = next;
+ }
+
+ multimap<double,double>::reverse_iterator rnext = fdist.rbegin();
+ multimap<double,double>::reverse_iterator rprev = rnext;
+ while ( ++rnext != fdist.rend() ) {
+ rnext->second += rprev->second;
+ if ( rnext->second > w5 && rprev->second <= w5 ) {
+ double f5 = rprev->first + (w5 - rprev->second)*
+ (rnext->first - rprev->first)/(rnext->second - rprev->second);
+ cout << "5% centrality for f > " << f5 << endl;
+ break;
+ }
+ rprev = rnext;
+ }
+}
diff --git a/DIPSY/CMS_2012_PAS_FWD_11_003.cc b/DIPSY/CMS_2012_PAS_FWD_11_003.cc
new file mode 100644
--- /dev/null
+++ b/DIPSY/CMS_2012_PAS_FWD_11_003.cc
@@ -0,0 +1,173 @@
+// Samantha Dooling DESY
+// February 2012
+//
+// -*- C++ -*-
+// =============================
+//
+// Ratio of the energy deposited in the pseudorapditiy range
+// -6.6 < eta < -5.2 for events with a charged particle jet
+//
+// =============================
+#include "Rivet/Analysis.hh"
+#include "Rivet/RivetAIDA.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/Projections/FinalState.hh"
+#include "Rivet/Projections/ChargedFinalState.hh"
+#include "Rivet/Projections/FastJets.hh"
+#include "Rivet/Projections/Beam.hh"
+#include "Rivet/Projections/VetoedFinalState.hh"
+#include "LWH/Histogram1D.h"
+
+namespace Rivet {
+
+
+ class CMS_2012_PAS_FWD_11_003 : public Analysis {
+ public:
+
+ /// Constructor
+ CMS_2012_PAS_FWD_11_003()
+ : Analysis("CMS_2012_PAS_FWD_11_003")
+ { }
+
+ void init() {
+
+ // gives the range of eta and min pT for the final state from which I get the jets
+ FastJets jetpro (ChargedFinalState(-2.5, 2.5, 0.3*GeV), FastJets::ANTIKT, 0.5);
+ addProjection(jetpro, "Jets");
+
+ // skip Neutrinos and Muons
+ VetoedFinalState fsv(FinalState(-7.0, -4.0, 0.*GeV));
+ fsv.vetoNeutrinos();
+ fsv.addVetoPairId(MUON);
+ addProjection(fsv, "fsv");
+
+ // for the hadron level selection
+ VetoedFinalState sfsv(FinalState(-MAXRAPIDITY, MAXRAPIDITY, 0.*GeV));
+ sfsv.vetoNeutrinos();
+ sfsv.addVetoPairId(MUON);
+ addProjection(sfsv, "sfsv");
+
+ //counters
+ passedSumOfWeights = 0.;
+ inclEflow = 0.;
+
+ // Temporary histograms to fill the energy flow for leading jet events.
+ // Ratios are calculated in finalyze().
+ int id = 0;
+ if (fuzzyEquals(sqrtS()/GeV, 900, 1e-3)) id=1;
+ if (fuzzyEquals(sqrtS()/GeV, 2760, 1e-3)) id=2;
+ if (fuzzyEquals(sqrtS()/GeV, 7000, 1e-3)) id=3;
+ _tmp_jet = bookHistogram1D("eflow_jet", binEdges(id, 1, 1)); // Leading jet energy flow in pt
+ _tmp_njet = bookHistogram1D("number_jet", binEdges(id, 1, 1)); // Number of events in pt
+ }
+
+
+ /// Perform the per-event analysis
+ void analyze(const Event& event) {
+
+ const double weight = event.weight();
+
+ // Skip if the event is empty
+ const FinalState& fsv = applyProjection<FinalState>(event, "fsv");
+ if (fsv.empty()) vetoEvent;
+
+ // ====================== Minimum Bias selection
+
+ const FinalState& sfsv = applyProjection<FinalState>(event, "sfsv");
+ // ParticleVector parts = sfsv.particlesByRapidity();
+ ParticleVector parts = sfsv.particles(cmpParticleByAscRapidity);
+ if (parts.empty()) vetoEvent;
+
+ // find dymax
+ double dymax = 0;
+ int gap_pos = -1;
+ for (size_t i=0; i < parts.size()-1; ++i) {
+ double dy = parts[i+1].momentum().rapidity() - parts[i].momentum().rapidity();
+ if (dy > dymax) {
+ dymax = dy;
+ gap_pos = i;
+ }
+ }
+
+ // calculate mx2 and my2
+ FourMomentum xmom;
+ for (int i=0; i<=gap_pos; ++i) {
+ xmom += parts[i].momentum();
+ }
+ double mx2 = xmom.mass2();
+ if (mx2<0) vetoEvent;
+
+ FourMomentum ymom;
+ for (size_t i=gap_pos+1; i<parts.size(); ++i) {
+ ymom += parts[i].momentum();
+ }
+ double my2 = ymom.mass2();
+ if (my2<0) vetoEvent;
+
+ // calculate xix and xiy and xidd
+ double xix = mx2 / sqr(sqrtS());
+ double xiy = my2 / sqr(sqrtS());
+ double xidd = mx2*my2 / sqr(sqrtS()*0.938*GeV);
+
+ // combine the selection: xi cuts
+ bool passedHadronCuts = false;
+ if (fuzzyEquals(sqrtS()/GeV, 900, 1e-3) && (xix > 0.1 || xiy > 0.4 || xidd > 0.5)) passedHadronCuts = true;
+ if (fuzzyEquals(sqrtS()/GeV, 2760, 1e-3) && (xix > 0.07 || xiy > 0.2 || xidd > 0.5)) passedHadronCuts = true;
+ if (fuzzyEquals(sqrtS()/GeV, 7000, 1e-3) && (xix > 0.04 || xiy > 0.1 || xidd > 0.5)) passedHadronCuts = true;
+ if (!passedHadronCuts) vetoEvent;
+
+ // ============================== MINIMUM BIAS EVENTS
+
+ // loop over particles to calculate the energy
+ passedSumOfWeights += weight;
+
+ foreach (const Particle& p, fsv.particles()) {
+ if (-5.2 > p.momentum().eta() && p.momentum().eta() > -6.6) inclEflow += weight*p.momentum().E()/GeV;
+ }
+
+ // ============================== JET EVENTS
+
+ const FastJets& jetpro = applyProjection<FastJets>(event, "Jets");
+ const Jets& jets = jetpro.jetsByPt(1.0*GeV);
+ if (jets.size()<1) vetoEvent;
+
+ if (fabs(jets[0].momentum().eta()) < 2.0) {
+ _tmp_njet->fill(jets[0].momentum().pT()/GeV, weight);
+
+ // energy flow
+ foreach (const Particle& p, fsv.particles()) {
+ if (p.momentum().eta() > -6.6 && p.momentum().eta() < -5.2) { // ask for the CASTOR region
+ _tmp_jet->fill(jets[0].momentum().pT()/GeV, weight * p.momentum().E()/GeV);
+ }
+ }
+ }
+
+ }// analysis
+
+ void finalize() {
+ _tmp_jet->scale(passedSumOfWeights/inclEflow);
+
+ AIDA::IHistogramFactory& hf = histogramFactory();
+ if (fuzzyEquals(sqrtS()/GeV, 900, 1e-3)) hf.divide(histoDir() + "/d01-x01-y01", *_tmp_jet, *_tmp_njet);
+ if (fuzzyEquals(sqrtS()/GeV, 2760, 1e-3)) hf.divide(histoDir() + "/d02-x01-y01", *_tmp_jet, *_tmp_njet);
+ if (fuzzyEquals(sqrtS()/GeV, 7000, 1e-3)) hf.divide(histoDir() + "/d03-x01-y01", *_tmp_jet, *_tmp_njet);
+
+ hf.destroy(_tmp_jet);
+ hf.destroy(_tmp_njet);
+ }
+
+ private:
+ // counters
+ double passedSumOfWeights;
+ double inclEflow;
+
+ // histograms
+ AIDA::IHistogram1D* _tmp_jet;
+ AIDA::IHistogram1D* _tmp_njet;
+ };
+
+
+ // The hook for the plugin system
+ DECLARE_RIVET_PLUGIN(CMS_2012_PAS_FWD_11_003);
+
+}
diff --git a/DIPSY/CMS_FWD_11_003_PRIVATE.cc b/DIPSY/CMS_FWD_11_003_PRIVATE.cc
new file mode 100644
--- /dev/null
+++ b/DIPSY/CMS_FWD_11_003_PRIVATE.cc
@@ -0,0 +1,456 @@
+// Samantha Dooling DESY
+// February 2012
+//
+// -*- C++ -*-
+// =============================
+//
+// Ratio of the energy deposited in the pseudorapditiy range
+// -6.6 < eta < -5.2 for events with a charged particle jet
+//
+// =============================
+#include "Rivet/Analysis.hh"
+#include "Rivet/RivetAIDA.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/Projections/FinalState.hh"
+#include "Rivet/Projections/ChargedFinalState.hh"
+#include "Rivet/Projections/FastJets.hh"
+#include "Rivet/Projections/Beam.hh"
+#include "Rivet/Projections/VetoedFinalState.hh"
+#include "LWH/Histogram1D.h"
+#include <multimap>
+
+namespace Rivet {
+
+
+ class CMS_FWD_11_003_PRIVATE : public Analysis {
+ public:
+
+ /// Constructor
+ CMS_FWD_11_003_PRIVATE()
+ : Analysis("CMS_FWD_11_003_PRIVATE")
+ {
+ setBeams(PROTON,PROTON);
+ }
+
+ // counters
+ double evcounter_incl;
+ double evcounter_jet;
+
+ double norm_incl_eflow_09;
+ double norm_incl_eflow_276;
+ double norm_incl_eflow_7;
+
+ void init() {
+
+ //gives the range of eta and min pT for the final state from which I get the jets
+ const ChargedFinalState fsj(-2.5,2.5,0.3*GeV);
+ addProjection(fsj, "FSJ");
+
+ FastJets jetpro (fsj, FastJets::ANTIKT, 0.5);
+ jetpro.useInvisibles();
+ addProjection(jetpro, "Jets");
+
+ //gives the range of eta and min pT for the final state
+ const FinalState fs(-7.0,-4.0,0.0*GeV);
+ addProjection(fs, "FS");
+ VetoedFinalState fsv(fs);
+ // skip Neutrinos and Muons
+ fsv.vetoNeutrinos();
+ fsv.addVetoPairDetail(MUON, 0.0*GeV, 99999.9*GeV);
+ addProjection(fsv, "fsv");
+
+ // for the hadron level selection
+ const FinalState sfs(-1000.0,1000.0,0.0*GeV);
+ addProjection(sfs, "sfs");
+ VetoedFinalState sfsv(sfs);
+ sfsv.vetoNeutrinos();
+ sfsv.addVetoPairDetail(MUON, 0.0*GeV, 99999.9*GeV);
+ addProjection(sfsv, "sfsv");
+
+ //counters
+ evcounter_incl = 0;
+ evcounter_jet = 0;
+
+ // pt binning
+ // int NptBins = 7;
+ // double ptbinning[8] = {1.,2.,3.,5.,7.5,10,15,25};
+
+ if(fuzzyEquals(sqrtS()/GeV, 900, 1E-3)){
+
+ // temporary histograms to fill the energy flow for inclusive events and leading jet events
+ // in finalyze() I determine the ratios
+
+ _tmp_incl_09.reset(new LWH::Histogram1D(binEdges(1,1,1))); // inclusive energy flow in the eta range of CASTOR
+ _tmp_jet_09 = bookHistogram1D("eflow_jet_09",binEdges(3,1,1)); // Leading jet energy flow in pt
+ _tmp_njet_09 = bookHistogram1D("number_jet_09",binEdges(3,1,1)); // Number of events in pt
+
+ _hist_incl_09 = bookHistogram1D(1,1,1); // Histogram: inclusive energy flow in the eta range of CASTOR
+ _hist_jet_09 = bookHistogram1D(2,1,1); // Histogram: energy flow of events with a charged particle jet in the eta range of CASTOR
+ }
+
+ if(fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)){
+
+ _tmp_incl_276.reset(new LWH::Histogram1D(binEdges(4,1,1)));
+ _tmp_jet_276 = bookHistogram1D("eflow_jet_276",binEdges(6,1,1));
+ _tmp_njet_276 = bookHistogram1D("number_jet_276",binEdges(6,1,1));
+
+ _hist_incl_276 = bookHistogram1D(4,1,1);
+ _hist_jet_276 = bookHistogram1D(5,1,1);
+ }
+
+
+ if(fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)){
+
+ _tmp_incl_7.reset(new LWH::Histogram1D(binEdges(7,1,1)));
+ _tmp_jet_7 = bookHistogram1D("eflow_jet_7",binEdges(9,1,1));
+ _tmp_njet_7 = bookHistogram1D("number_jet_7",binEdges(9,1,1));
+
+ _hist_incl_7 = bookHistogram1D(7,1,1);
+ _hist_jet_7 = bookHistogram1D(8,1,1);
+ }
+
+ }
+
+
+ /// Perform the per-event analysis
+ void analyze(const Event& event) {
+
+ const double weight = event.weight();
+
+ // Skip if the event is empty
+ const FinalState& fsv = applyProjection<FinalState>(event, "fsv");
+ if (fsv.empty()) vetoEvent;
+
+ const FastJets& jetpro = applyProjection<FastJets>(event, "Jets");
+ const Jets& jets = jetpro.jetsByPt(1.0*GeV);
+
+ // ====================== Minimum Bias selection
+ // ============================== xi cuts
+
+ bool passedHadronCuts = false;
+
+ double xix = 10;
+ double xiy = 10;
+ double xidd = 10e10;
+ double Rapiditymax = -1;
+
+ // calculate xi of the event
+ // sort Particles in rapidity, from rapidity_min to rapidity_max
+
+ ParticleVector myTempParticles;
+ ParticleVector myRapiditySortedParticles;
+
+ // copy only final stable particles in tempvector
+ const FinalState& sfsv = applyProjection<FinalState>(event, "sfsv");
+ if (sfsv.empty()) vetoEvent;
+
+ foreach (const Particle& p, sfsv.particles()) {
+ myTempParticles.push_back(Particle(p));
+ }
+
+ while (myTempParticles.size() != 0){
+ double min_y = myTempParticles[0].momentum().rapidity();
+ int min_y_pos = 0;
+ for (unsigned int ipart = 1; ipart < myTempParticles.size(); ++ipart){
+ if (myTempParticles[ipart].momentum().rapidity() < min_y){
+ min_y = myTempParticles[ipart].momentum().rapidity();
+ min_y_pos = ipart;
+ }
+ }
+ myRapiditySortedParticles.push_back(myTempParticles[min_y_pos]);
+ myTempParticles.erase(myTempParticles.begin()+min_y_pos);
+ }
+
+
+ // find deltaymax
+ double deltaymax = 0;
+ int deltaymax_pos = -1;
+ for (unsigned int ipart=0; ipart < myRapiditySortedParticles.size()-1; ++ipart) {
+ double deltay = myRapiditySortedParticles[ipart+1].momentum().rapidity() - myRapiditySortedParticles[ipart].momentum().rapidity();
+ if (deltay > deltaymax) {
+ deltaymax = deltay;
+ deltaymax_pos = ipart;
+ }
+ }
+ Rapiditymax = deltaymax;
+
+ // calculate Mx2 and My2
+ FourMomentum Xfourmom;
+ FourMomentum Yfourmom;
+
+ for (int ipart=0; ipart <= deltaymax_pos; ++ipart) {
+ Xfourmom += myRapiditySortedParticles[ipart].momentum();
+ }
+ if(FourMomentum(Xfourmom).mass2() <0 )
+ vetoEvent;
+
+ long double Mx2 = FourMomentum(Xfourmom).mass()*FourMomentum(Xfourmom).mass();
+
+ for (unsigned int ipart = deltaymax_pos+1; ipart < myRapiditySortedParticles.size(); ++ipart) {
+ Yfourmom += myRapiditySortedParticles[ipart].momentum();
+ }
+ if(FourMomentum(Yfourmom).mass2() <0 )
+ vetoEvent;
+
+ long double My2 = FourMomentum(Yfourmom).mass()*FourMomentum(Yfourmom).mass() ;
+
+ // calculate xix and xiy and xidd
+ xix = Mx2/(sqrtS()/GeV*sqrtS()/GeV);
+ xiy = My2/(sqrtS()/GeV*sqrtS()/GeV);
+ xidd = (Mx2*My2)/(sqrtS()/GeV*sqrtS()/GeV*0.938*0.938);
+
+ // combine the selection
+ if(fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) {
+ if(xix > 0.1 || xiy > 0.4 || xidd > 0.5) passedHadronCuts = true;
+ }
+ if(fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)) {
+ if(xix > 0.07 || xiy > 0.2 || xidd > 0.5) passedHadronCuts = true;
+ }
+ if(fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) {
+ if(xix > 0.04 || xiy > 0.1 || xidd > 0.5) passedHadronCuts = true;
+ }
+
+ // skip the event if the hadron cut is not fullfilled
+ if(passedHadronCuts == false){
+ vetoEvent;
+ }
+
+ // ============================== MINIMUM BIAS EVENTS
+
+ // loop over particles to calculate the energy
+ evcounter_incl += weight;
+
+ foreach (const Particle& p, fsv.particles()) {
+
+ if(fuzzyEquals(sqrtS()/GeV, 900, 1E-3)) {
+
+ _tmp_incl_09 ->
+ fill(p.momentum().pseudorapidity(), weight * p.momentum().E()/GeV );
+ _hist_incl_09 ->
+ fill(p.momentum().pseudorapidity(), weight * p.momentum().E()/GeV );
+ }
+
+ if(fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)) {
+
+ _tmp_incl_276 ->
+ fill(p.momentum().pseudorapidity(), weight * p.momentum().E()/GeV );
+ _hist_incl_276 ->
+ fill(p.momentum().pseudorapidity(), weight * p.momentum().E()/GeV );
+ }
+
+ if(fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)) {
+
+ _tmp_incl_7 ->
+ fill(p.momentum().pseudorapidity(), weight * p.momentum().E()/GeV );
+ _hist_incl_7 ->
+ fill(p.momentum().pseudorapidity(), weight * p.momentum().E()/GeV );
+ }
+ }
+
+ // ============================== JET EVENTS
+
+ if(jets.size()>0){
+
+ signed int index_1 = -1; // for the jet with the 1.highest pT
+ double tempmax_1 = -100;
+
+ // jet with the 1.highest pt
+ for(signed int ijets = 0; ijets < (int)jets.size(); ++ijets){
+ if(tempmax_1 == -100 || tempmax_1 < jets[ijets].momentum().pT()/GeV){
+ tempmax_1 = jets[ijets].momentum().pT()/GeV;
+ index_1 = ijets;
+ }
+ }
+
+ if(index_1 != -1 ){
+
+ // *******
+ // 900 GeV
+ // *******
+
+ if(fuzzyEquals(sqrtS()/GeV, 900, 1E-3)){
+
+ //eta cut for the central jets
+ if(fabs(jets[index_1].momentum().pseudorapidity()) < 2.0 ){
+
+ // fill number of events
+ _tmp_njet_09->
+ fill(jets[index_1].momentum().pT()/GeV, weight );
+
+ // energy flow
+ foreach (const Particle& p, fsv.particles()){
+
+ // ask for the CASTOR region
+ if(-5.2 > p.momentum().pseudorapidity() && p.momentum().pseudorapidity() > -6.6 ){
+ _tmp_jet_09 ->
+ fill(jets[index_1].momentum().pT()/GeV, weight * p.momentum().E()/GeV );
+ }
+ }
+
+ // pt cut for the energy flow
+ if(jets[index_1].momentum().pT()/GeV > 10.0) {
+ evcounter_jet += weight;
+
+ // energy flow
+ foreach (const Particle& p, fsv.particles()){
+ _hist_jet_09 ->
+ fill(p.momentum().pseudorapidity(), weight * p.momentum().E()/GeV );
+ }
+ }// pt cut
+ }// eta
+ }// energy
+
+ // *******
+ // 2760 GeV
+ // *******
+
+ if(fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)){
+
+ // eta cut for the central jets
+ if(fabs(jets[index_1].momentum().pseudorapidity()) < 2.0 ){
+
+ // fill number of events
+ _tmp_njet_276->
+ fill(jets[index_1].momentum().pT()/GeV, weight );
+
+ // energy flow
+ foreach (const Particle& p, fsv.particles()){
+
+ // ask for the CASTOR region
+ if(-5.2 > p.momentum().pseudorapidity() && p.momentum().pseudorapidity() > -6.6 ){
+ _tmp_jet_276 ->
+ fill(jets[index_1].momentum().pT()/GeV, weight * p.momentum().E()/GeV );
+ }
+ }
+
+ // pt cut for the energy flow
+ if(jets[index_1].momentum().pT()/GeV > 10.0) {
+ evcounter_jet += weight;
+
+ // energy flow
+ foreach (const Particle& p, fsv.particles()){
+ _hist_jet_276 ->
+ fill(p.momentum().pseudorapidity(), weight * p.momentum().E()/GeV );
+ }
+ }// pt cut
+ }// eta
+ }// energy
+
+ // *******
+ // 7 TeV
+ // *******
+
+ if(fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)){
+
+ // eta cut for the central jets
+ if(fabs(jets[index_1].momentum().pseudorapidity()) < 2.0 ){
+
+ // fill number of events
+ _tmp_njet_7->
+ fill(jets[index_1].momentum().pT()/GeV, weight );
+
+ // energy flow
+ foreach (const Particle& p, fsv.particles()){
+
+ // ask for the CASTOR region
+ if(-5.2 > p.momentum().pseudorapidity() && p.momentum().pseudorapidity() > -6.6 ){
+ _tmp_jet_7 ->
+ fill(jets[index_1].momentum().pT()/GeV, weight * p.momentum().E()/GeV );
+ }
+ }
+
+ // pt cut for the energy flow
+ if(jets[index_1].momentum().pT()/GeV > 10.0) {
+ evcounter_jet += weight;
+
+ // energy flow
+ foreach (const Particle& p, fsv.particles()){
+ _hist_jet_7 ->
+ fill(p.momentum().pseudorapidity(), weight * p.momentum().E()/GeV );
+ }
+ }// pt cut
+ }// eta
+ }// energy
+ }// if jet index
+ }// if jets
+
+ }// analysis
+
+ void finalize() {
+
+ AIDA::IHistogramFactory& hf = histogramFactory();
+
+ const double norm_jet = evcounter_jet ;
+ const double norm_incl = evcounter_incl ;
+ const double binwidth = (6.6 - 5.2);
+
+ if(fuzzyEquals(sqrtS()/GeV, 900, 1E-3)){
+
+ norm_incl_eflow_09 = _tmp_incl_09->binHeight(0)/norm_incl; // no normalization to the binwidth
+ _tmp_jet_09->scale(1.0/norm_incl_eflow_09);
+
+ //the energy flow ratio
+ hf.divide(histoDir() + "/d03-x01-y01", *_tmp_jet_09, *_tmp_njet_09);
+ hf.destroy(_tmp_jet_09);
+ hf.destroy(_tmp_njet_09);
+
+ scale(_hist_incl_09,binwidth/norm_incl); //scale automatically normalizes to binwidth but I want to skip it here
+ scale(_hist_jet_09,binwidth/norm_jet);
+ }
+
+ if(fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)){
+
+ norm_incl_eflow_276 = _tmp_incl_276->binHeight(0)/norm_incl;
+ _tmp_jet_276->scale(1.0/norm_incl_eflow_276);
+
+ //the energy flow ratio
+ hf.divide(histoDir() + "/d06-x01-y01", *_tmp_jet_276, *_tmp_njet_276);
+ hf.destroy(_tmp_jet_276);
+ hf.destroy(_tmp_njet_276);
+
+ scale(_hist_incl_276,binwidth/norm_incl);
+ scale(_hist_jet_276,binwidth/norm_jet);
+ }
+
+ if(fuzzyEquals(sqrtS()/GeV, 7000, 1E-3)){
+
+ norm_incl_eflow_7 = _tmp_incl_7->binHeight(0)/norm_incl;
+ _tmp_jet_7->scale(1.0/norm_incl_eflow_7);
+
+ //the energy flow ratio
+ hf.divide(histoDir() + "/d09-x01-y01", *_tmp_jet_7, *_tmp_njet_7);
+ hf.destroy(_tmp_jet_7);
+ hf.destroy(_tmp_njet_7);
+
+ scale(_hist_incl_7,binwidth/norm_incl);
+ scale(_hist_jet_7,binwidth/norm_jet);
+ }
+
+ MSG_INFO(" " );
+ MSG_INFO("Number of inclusive events : " << norm_incl );
+ MSG_INFO("Number of jet events : " << norm_jet );
+
+ }
+
+ private:
+ shared_ptr<LWH::IHistogram1D> _tmp_incl_09;
+ AIDA::IHistogram1D *_tmp_jet_09,*_tmp_njet_09;
+ AIDA::IHistogram1D *_hist_incl_09, *_hist_jet_09;
+
+ shared_ptr<LWH::IHistogram1D> _tmp_incl_276;
+ AIDA::IHistogram1D *_tmp_jet_276, *_tmp_njet_276;
+ AIDA::IHistogram1D *_hist_incl_276, *_hist_jet_276;
+
+ shared_ptr<LWH::IHistogram1D> _tmp_incl_7;
+ AIDA::IHistogram1D *_tmp_jet_7, *_tmp_njet_7;
+ AIDA::IHistogram1D *_hist_incl_7, *_hist_jet_7;
+
+ };
+
+
+
+ // This global object acts as a hook for the plugin system
+ AnalysisBuilder<CMS_FWD_11_003_PRIVATE> plugin_CMS_FWD_11_003_PRIVATE;
+
+
+}
diff --git a/DIPSY/Makefile.am b/DIPSY/Makefile.am
--- a/DIPSY/Makefile.am
+++ b/DIPSY/Makefile.am
@@ -1,128 +1,128 @@
AUTOMAKE_OPTIONS = -Wno-portability
mySOURCES = DipoleEventHandler.cc WaveFunction.cc SimpleProton.cc \
VirtualPhoton.cc DipoleState.cc Dipole.cc Parton.cc DipoleXSec.cc \
ImpactParameterGenerator.cc WFInfo.cc PhotonWFInfo.cc \
VectorMesonBase.cc PhotonDipoleState.cc SimpleProtonState.cc \
Emitter.cc Swinger.cc DipoleAnalysisHandler.cc \
TotalXSecAnalysis.cc EventFiller.cc DipoleAbsorber.cc \
SmallDipoleAbsorber.cc EffectiveParton.cc RealParton.cc \
RealPartonState.cc DiffractiveEventFiller.cc
DOCFILES = DipoleEventHandler.h WaveFunction.h SimpleProton.h VirtualPhoton.h \
DipoleState.h Dipole.h Parton.h DipoleXSec.h ImpactParameters.h \
ImpactParameterGenerator.h WFInfo.h PhotonWFInfo.h VectorMesonBase.h \
PhotonDipoleState.h SimpleProtonState.h Emitter.h Swinger.h \
DipoleAnalysisHandler.h TotalXSecAnalysis.h EventFiller.h \
DipoleAbsorber.h SmallDipoleAbsorber.h EffectiveParton.h \
RealParton.h RealPartonState.h DiffractiveEventFiller.h
INCLUDEFILES = $(DOCFILES) DipoleEventHandler.icc DipoleEventHandler.fh \
WaveFunction.icc WaveFunction.fh SimpleProton.icc \
VirtualPhoton.icc DipoleState.fh DipoleState.icc Dipole.fh \
Dipole.icc Parton.fh Parton.icc DipoleXSec.fh DipoleXSec.icc \
ImpactParameters.icc ImpactParameterGenerator.fh \
WFInfo.fh WFInfo.icc \
PhotonWFInfo.icc VectorMesonBase.icc PhotonDipoleState.icc \
SimpleProtonState.icc Emitter.fh Emitter.icc Swinger.fh \
Swinger.icc DipoleAnalysisHandler.fh EventFiller.fh \
DipoleAbsorber.fh EffectiveParton.fh \
RealParton.fh RealPartonState.fh DiffractiveEventFiller.fh
pkglib_LTLIBRARIES = libDIPSY.la OldStyleEmitter.la PT1DEmitter.la \
ElasticXSecAnalysis.la LargePTDipoleAbsorber.la \
SimpleNucleus.la RecoilSwinger.la FSDipoleOrdering.la \
FSDipole5Ordering.la FSAnalysis.la GapAnalysis.la HIAnalysis.la \
FixedImpactGenerator.la
INPUTFILES = TestXSecs.in RivetAnalyses.in DIPSYRemove.in DIPSYDefaults.in CurrentTune.in Tune27.in
CLEANFILES = done-all-links
dist_pkgdata_DATA = $(INPUTFILES)
# Version info should be updated if any interface or persistent I/O
# function is changed
libDIPSY_la_LDFLAGS = -module -version-info 1:0:0
libDIPSY_la_SOURCES = $(mySOURCES) $(INCLUDEFILES)
OldStyleEmitter_la_LDFLAGS = -module -version-info 1:0:0
OldStyleEmitter_la_SOURCES = OldStyleEmitter.cc OldStyleEmitter.h
FSDipoleOrdering_la_LDFLAGS = -module -version-info 1:0:0
FSDipoleOrdering_la_SOURCES = FSDipoleOrdering.cc FSDipoleOrdering.h
FSDipole5Ordering_la_LDFLAGS = -module -version-info 1:0:0
FSDipole5Ordering_la_SOURCES = FSDipole5Ordering.cc FSDipole5Ordering.h
FSAnalysis_la_LDFLAGS = -module -version-info 1:0:0
FSAnalysis_la_SOURCES = FSAnalysis.cc FSAnalysis.h
GapAnalysis_la_LDFLAGS = -module -version-info 1:0:0
GapAnalysis_la_SOURCES = GapAnalysis.cc GapAnalysis.h
HIAnalysis_la_LDFLAGS = -module -version-info 1:0:0
HIAnalysis_la_SOURCES = HIAnalysis.cc HIAnalysis.h
PT1DEmitter_la_LDFLAGS = -module -version-info 1:0:0
PT1DEmitter_la_SOURCES = PT1DEmitter.cc PT1DEmitter.h
ElasticXSecAnalysis_la_LDFLAGS = -module -version-info 1:0:0
ElasticXSecAnalysis_la_SOURCES = ElasticXSecAnalysis.cc ElasticXSecAnalysis.h
LargePTDipoleAbsorber_la_LDFLAGS = -module -version-info 1:0:0
LargePTDipoleAbsorber_la_SOURCES = LargePTDipoleAbsorber.cc LargePTDipoleAbsorber.h
SimpleNucleus_la_LDFLAGS = -module -version-info 1:0:0
SimpleNucleus_la_SOURCES = SimpleNucleus.cc SimpleNucleus.h \
SimpleNucleusState.cc SimpleNucleusState.h
RecoilSwinger_la_LDFLAGS = -module -version-info 1:0:0
RecoilSwinger_la_SOURCES = RecoilSwinger.cc RecoilSwinger.h
FixedImpactGenerator_la_LDFLAGS = -module -version-info 1:0:0
FixedImpactGenerator_la_SOURCES = FixedImpactGenerator.cc FixedImpactGenerator.h
EXTRACLASSES = OldStyleEmitter.la PT1DEmitter.la ElasticXSecAnalysis.la \
LargePTDipoleAbsorber.la RecoilSwinger.la SimpleNucleus.la \
FSDipoleOrdering.la FSAnalysis.la GapAnalysis.la HIAnalysis.la \
FixedImpactGenerator.la FSDipole5Ordering.la
done-all-links:
@EMPTY@ifdef SHOWCOMMAND
for file in $(INPUTFILES); do \
if test ! -f $$file; then $(LN_S) $(srcdir)/$$file $$file; fi; done
echo "stamp" > done-all-links
@EMPTY@else
@echo "sym-linking input files files..."
@for file in $(INPUTFILES); do \
if test ! -f $$file; then $(LN_S) $(srcdir)/$$file $$file; fi; done
@echo "stamp" > done-all-links
@EMPTY@endif
DIPSYDefaults.rpo: done-all-links DIPSYDefaults.in libDIPSY.la $(EXTRACLASSES) $(THEPEGLIB)/ThePEGDefaults.rpo
$(SETUPTHEPEG) -L../lib -L$(THEPEGLIB) -L ../../TheP8I/lib -L .libs -L ../src/.libs --exitonerror --init -r ../lib/Ariadne5Defaults.rpo -o DIPSYDefaults.rpo DIPSYDefaults.in
Test5%.run: Test5%.in DIPSYDefaults.rpo CurrentTune.in RivetAnalyses.in
$(SETUPTHEPEG) --exitonerror -L .libs -L ../../Pythia7/src/.libs -r DIPSYDefaults.rpo $<
%.run: %.in DIPSYDefaults.rpo CurrentTune.in RivetAnalyses.in
- $(SETUPTHEPEG) --exitonerror -L .libs -r DIPSYDefaults.rpo $<
+ $(SETUPTHEPEG) -L../lib --exitonerror -L .libs -r DIPSYDefaults.rpo $<
%.out: %.run
time $(RUNTHEPEG) --tics -d 0 $<
valgrind:
valgrind --leak-check=full --num-callers=25 --track-fds=yes --freelist-vol=100000000 --leak-resolution=med --trace-children=yes $(top_builddir)/../ThePEG/src/setupThePEG -L ../../ThePEG/lib -L../../Pythia7/lib --exitonerror -L../../TheP8I/lib -L .libs -L ../lib -r ../lib/Ariadne5Defaults.rpo TestFull.in &> /tmp/valgrind.out
valgrind --leak-check=full --num-callers=25 --track-fds=yes --freelist-vol=100000000 --leak-resolution=med --trace-children=yes $(top_builddir)/../ThePEG/src/runThePEG -N 10 TestFull.run >> /tmp/valgrind.out 2>&1
install-data-local:
LD_LIBRARY_PATH=$(DESTDIR)$(pkglibdir):$$LD_LIBRARY_PATH $(DESTDIR)$(bindir)/setupThePEG --exitonerror --init -i $(DESTDIR)$(pkgdatadir) -r $(DESTDIR)$(pkglibdir)/ThePEGDefaults.rpo -o $(DESTDIR)$(pkglibdir)/ThePEGDefaults.rpo $(srcdir)/DIPSYDefaults.in
unregister-from-repo:
LD_LIBRARY_PATH=$(DESTDIR)$(pkglibdir):$$LD_LIBRARY_PATH $(DESTDIR)$(bindir)/setupThePEG --exitonerror --init -r $(DESTDIR)$(pkglibdir)/ThePEGDefaults.rpo -o $(DESTDIR)$(pkglibdir)/ThePEGDefaults.rpo $(srcdir)/DIPSYRemove.in
include $(top_srcdir)/Config/Makefile.aminclude
diff --git a/DIPSY/RHICAuAu.in b/DIPSY/RHICAuAu.in
new file mode 100644
--- /dev/null
+++ b/DIPSY/RHICAuAu.in
@@ -0,0 +1,140 @@
+cd /DIPSY
+
+read Tune30.in
+set stdEmitter:MinusOrderingMode EffectiveParton
+set stdEmitter:PMinusOrdering 1.5
+set stdProton:R0 2.7
+set stdAntiProton:R0 2.7
+mset . DIPSY::DipoleEventHandler LambdaQCD 0.22
+
+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
+
+erase RHICGenerator:AnalysisHandlers[0]
+
+library HIAnalysis.so
+create DIPSY::HIAnalysis HIAna
+create ThePEG::ProgressLog Logger ProgressLog.so
+set Logger:Interval 600
+
+set RHICGenerator:NumberOfEvents 100
+set RHICEventHandler:PreSamples 0
+set RHICEventHandler:ConsistencyLevel 0
+set RHICEventHandler:DecayHandler:MaxLifeTime 10
+set RHICEventHandler:BGen:Width 25
+set RHICEventHandler:XSecFn:CheckOffShell false
+insert RHICGenerator:AnalysisHandlers[0] Logger
+insert RHICGenerator:AnalysisHandlers[0] HIAna
+set RHICEventHandler:LuminosityFunction:Energy 39400
+set RHICGenerator:EventHandler:EventFiller:PTCut 0.6
+set RHICEventHandler:EventFiller:FSSwingTime 1.0
+set RHICEventHandler:EventFiller:FSSwingTimeStep 0.1
+
+saverun RHICAuAu RHICGenerator
+
+set RHICEventHandler:EventFiller:FSSwingTime 0.0
+saverun RHICAuAu0 RHICGenerator
+set RHICEventHandler:EventFiller:FSSwingTime 1.0
+
+
+create DIPSY::FixedImpactGenerator FixedB FixedImpactGenerator.so
+set RHICEventHandler:BGen FixedB
+set FixedB:MaxB 3.2
+set FixedB:MinB 0.0
+saverun RHICAuAuC RHICGenerator
+set RHICEventHandler:EventFiller:FSSwingTime 0.0
+saverun RHICAuAuC0 RHICGenerator
+set stdSwinger:Lambda 2.0
+saverun RHICAuAuC02 RHICGenerator
+set stdSwinger:Lambda 5.0
+saverun RHICAuAuC05 RHICGenerator
+set stdSwinger:Lambda 20.0
+saverun RHICAuAuC020 RHICGenerator
+set stdSwinger:Lambda 100.0
+saverun RHICAuAuC0100 RHICGenerator
+set stdSwinger:Lambda 1.0
+
+# New Final state swing
+# With maximum energy scale
+create Ariadne5::DipoleSwinger Swinger
+set Swinger:Rmax 0.0
+insert RHICGenerator:EventHandler:CascadeHandler:Emitters[0] Swinger
+saverun RHICAuAuCNM1 RHICGenerator
+set Swinger:Lambda 2
+saverun RHICAuAuCNM2 RHICGenerator
+set Swinger:Lambda 10
+saverun RHICAuAuCNM10 RHICGenerator
+set Swinger:Lambda 50
+saverun RHICAuAuCNM50 RHICGenerator
+set Swinger:Lambda 100
+saverun RHICAuAuCNM100 RHICGenerator
+set Swinger:Lambda 1
+set Swinger:Rmax 2.7
+saverun RHICAuAuCNM1R RHICGenerator
+set Swinger:Lambda 10
+saverun RHICAuAuCNM10R RHICGenerator
+set Swinger:Lambda 1
+set Swinger:Rmax 0.0
+set Swinger:Linear Linear
+saverun RHICAuAuCNM1L RHICGenerator
+set Swinger:Lambda 10
+saverun RHICAuAuCNM10L RHICGenerator
+set Swinger:Lambda 1
+set Swinger:Linear Logarithmic
+
+# With maximum pt scale
+create Ariadne5::MaxPTScale MaxPTScale MaxPTScale.so
+set MaxPTScale:Strategy ActualPT
+set RHICGenerator:EventHandler:CascadeHandler:ScaleSetter MaxPTScale
+saverun RHICAuAuCN1 RHICGenerator
+set Swinger:Lambda 2
+saverun RHICAuAuCN2 RHICGenerator
+set Swinger:Lambda 0.1
+saverun RHICAuAuCN01 RHICGenerator
+set Swinger:Lambda 5
+saverun RHICAuAuCN5 RHICGenerator
+set Swinger:Lambda 10
+saverun RHICAuAuCN10 RHICGenerator
+set Swinger:Lambda 20
+saverun RHICAuAuCN20 RHICGenerator
+set Swinger:Lambda 50
+saverun RHICAuAuCN50 RHICGenerator
+set Swinger:Lambda 100
+saverun RHICAuAuCN100 RHICGenerator
+set Swinger:Lambda 200
+saverun RHICAuAuCN200 RHICGenerator
+set Swinger:Rmax 2.7
+set Swinger:Lambda 1
+saverun RHICAuAuCN1R RHICGenerator
+set Swinger:Lambda 10
+saverun RHICAuAuCN10R RHICGenerator
+set Swinger:Lambda 20
+saverun RHICAuAuCN20R RHICGenerator
+set Swinger:Linear Linear
+set Swinger:Lambda 1
+saverun RHICAuAuCN1L RHICGenerator
+set Swinger:Lambda 10
+saverun RHICAuAuCN10L RHICGenerator
+set Swinger:Lambda 100
+saverun RHICAuAuCN100L RHICGenerator
+set Swinger:Linear Logarithmic
+set Swinger:Rmax 0.0
+set Swinger:Lambda 1
+
+# With scalar sum pt scale
+set MaxPTScale:Strategy ScalarSum
+saverun RHICAuAuCNS1 RHICGenerator
+set Swinger:Lambda 2
+saverun RHICAuAuCNS2 RHICGenerator
+set Swinger:Lambda 10
+saverun RHICAuAuCNS10 RHICGenerator
+set Swinger:Rmax 2.7
+set Swinger:Linear Linear
+saverun RHICAuAuCNS10L RHICGenerator
+set Swinger:Lambda 100
+saverun RHICAuAuCNS100L RHICGenerator
+
+
diff --git a/DIPSY/Test5ppFS.in b/DIPSY/Test5ppFS.in
--- a/DIPSY/Test5ppFS.in
+++ b/DIPSY/Test5ppFS.in
@@ -1,237 +1,270 @@
cd /DIPSY
+# read Tune30.in
+# set stdEmitter:MinusOrderingMode EffectiveParton
+# set stdEmitter:PMinusOrdering 1.5
+# set stdProton:R0 2.7
+# set stdAntiProton:R0 2.7
+# mset . DIPSY::DipoleEventHandler LambdaQCD 0.22
+
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 LHCGenerator:NumberOfEvents 10000
set LHCEventHandler:PreSamples 0
set LHCEventHandler:ConsistencyLevel Never
set LHCEventHandler:DecayHandler:MaxLifeTime 10
insert LHCGenerator:AnalysisHandlers[0] Logger
set LHCGenerator:MaxErrors 10000
set TevatronGenerator:NumberOfEvents 1000000
set TevatronEventHandler:PreSamples 0
set TevatronEventHandler:DecayHandler:MaxLifeTime 10
insert TevatronGenerator:AnalysisHandlers[0] Logger
set TevatronGenerator:MaxErrors 10000
set LHCGenerator:EventHandler:EventFiller:PTCut 0.6
set LHCLumi:Energy 70.0
#saverun RHIC LHCGenerator
#saverun RHICCuAu LHCGenerator
#saverun LHCPbPb LHCGenerator
set LHCGenerator:RandomNumberGenerator:Seed 1
set LHCLumi:Energy 7000.0
saverun Test5ppFS70debug LHCGenerator
-insert LHCGenerator:AnalysisHandlers[0] RivetLHC
-saverun LHCpp30d LHCGenerator
+create ThePEG::RivetAnalysis RivetCMS
+insert RivetCMS:Analyses[0] CMS_2012_PAS_FWD_11_003
+insert RivetCMS:Analyses[0] CMS_FWD_11_003_PRIVATE
+insert LHCGenerator:AnalysisHandlers[0] RivetCMS
+set LHCGenerator:NumberOfEvents 1000000
+saverun Test5ppFS70CMS LHCGenerator
+set LHCLumi:Energy 2760.0
+saverun Test5ppFS27CMS LHCGenerator
+set LHCLumi:Energy 900.0
+saverun Test5ppFS09CMS LHCGenerator
+set LHCLumi:Energy 7000.0
+set LHCGenerator:EventHandler:YFrametest 0.3
+saverun Test5ppFS70CMS03 LHCGenerator
+set LHCLumi:Energy 2760.0
+saverun Test5ppFS27CMS03 LHCGenerator
+set LHCLumi:Energy 900.0
+saverun Test5ppFS09CMS03 LHCGenerator
+set LHCLumi:Energy 7000.0
+set LHCGenerator:EventHandler:YFrametest 0.5
+set LHCGenerator:EventHandler:YFrametest 0.7
+saverun Test5ppFS70CMS07 LHCGenerator
+set LHCLumi:Energy 2760.0
+saverun Test5ppFS27CMS07 LHCGenerator
+set LHCLumi:Energy 900.0
+saverun Test5ppFS09CMS07 LHCGenerator
+set LHCLumi:Energy 7000.0
+set LHCGenerator:EventHandler:YFrametest 0.5
+set LHCGenerator:NumberOfEvents 10000
+# saverun LHCpp30d LHCGenerator
set LHCGenerator:EventHandler:YFrametest 0.3
saverun LHCpp30cFrame03 LHCGenerator
set LHCGenerator:EventHandler:YFrametest 0.5
set LHCGenerator:RandomNumberGenerator:Seed 1
set LHCGenerator:NumberOfEvents 0
set LHCEventHandler:PreSamples 4000000
set LHCLumi:Energy 100.0
saverun pp100Inc LHCGenerator
set LHCLumi:Energy 200.0
#saverun pp200Inc LHCGenerator
set LHCLumi:Energy 400.0
#saverun pp400Inc LHCGenerator
set LHCLumi:Energy 900.0
#saverun pp900Inc LHCGenerator
set LHCLumi:Energy 1800.0
#saverun pp1800Inc LHCGenerator
set LHCLumi:Energy 2360.0
#saverun pp2360Inc LHCGenerator
set LHCLumi:Energy 7000.0
saverun pp7000Inc4MSeed1 LHCGenerator
set LHCLumi:Energy 14000.0
#saverun pp14000Inc LHCGenerator
create DIPSY::FixedImpactGenerator FixedB FixedImpactGenerator.so
set FixedB:MaxB 1.28
set FixedB:MinB 1.28
set FixedB:BAngle 0
set LHCEventHandler:BGen FixedB
set LHCLumi:Energy 200.0
#saverun ppW200B128temp LHCGenerator
set LHCLumi:Energy 7000.0
set LHCEventHandler:BGen stdBGenerator
#saverun Test5ppFS23 LHCGenerator
set LHCLumi:Energy 900.0
#insert LHCGenerator:AnalysisHandlers[0] RivetLHC
set LHCLumi:Energy 900.0
#saverun Test5ppFS09 LHCGenerator
set LHCGenerator:EventHandler:EventFiller:PTCut 0.6
#saverun Test5ppFS09cut LHCGenerator
set /DIPSY/FSOrdering:OnlyOriginal OnlyOriginalEmissions
#saverun Test5ppFS09orig LHCGenerator
#set /DIPSY/stdXSec:Smoothing 0.1
#saverun Test5ppFS09smooth LHCGenerator
#set /DIPSY/stdXSec:Smoothing 0.0
set /DIPSY/FSOrdering:Fudge 1.2
#saverun Test5ppFS09fudge LHCGenerator
set /DIPSY/FSOrdering:Fudge 1.0
set /DIPSY/FSOrdering:PTMin 2.0
set LHCGenerator:EventHandler:EventFiller:PTCut 2.0
#saverun Test5ppFS09ptlim LHCGenerator
set LHCGenerator:EventHandler:EventFiller:PTCut 0.6
set /DIPSY/FSOrdering:PTMin 0.0
set /DIPSY/FSOrdering:Fudge 1.0
set /DIPSY/FSOrdering:OnlyOriginal AllEmissions
set LHCLumi:Energy 7000.0
#set LHCGenerator:RandomNumberGenerator:Seed 0
#saverun LHCpp0 LHCGenerator
#...
#set LHCGenerator:RandomNumberGenerator:Seed 20
#saverun LHCpp20 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 21
#saverun LHCpp21 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 22
#saverun LHCpp22 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 23
#saverun LHCpp23 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 24
#saverun LHCpp24 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 25
#saverun LHCpp25 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 26
#saverun LHCpp26 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 27
#saverun LHCpp27 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 28
#saverun LHCpp28 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 29
#saverun LHCpp29 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 30
#saverun LHCpp30 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 31
#saverun LHCpp31 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 32
#saverun LHCpp32 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 33
#saverun LHCpp33 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 34
#saverun LHCpp34 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 35
#saverun LHCpp35 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 36
#saverun LHCpp36 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 37
#saverun LHCpp37 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 38
#saverun LHCpp38 LHCGenerator
#set LHCGenerator:RandomNumberGenerator:Seed 39
#saverun LHCpp39 LHCGenerator
set LHCGenerator:EventHandler:YFrametest 0.3
set LHCLumi:Energy 900.0
#saverun Test5ppFS09frame03 LHCGenerator
set LHCLumi:Energy 7000.0
set LHCGenerator:NumberOfEvents 200000
#saverun Test5ppFS70frame03 LHCGenerator
set LHCGenerator:EventHandler:YFrametest 0.5
set TevatronLumi:Energy 1800.0
insert TevatronGenerator:AnalysisHandlers[0] RivetCDF
#saverun Test5ppFS18 TevatronGenerator
set TevatronGenerator:EventHandler:YFrametest 0.3
#saverun Test5ppFS18frame03 TevatronGenerator
set TevatronGenerator:EventHandler:YFrametest 0.5
cp LHCGenerator LHCcentral
create DIPSY::FixedImpactGenerator FixedB FixedImpactGenerator.so
set FixedB:MaxB = 0.0
set FixedB:MinB = 0.0
set LHCcentral:EventHandler:BGen FixedB
erase LHCcentral:AnalysisHandlers[0]
#erase LHCcentral:AnalysisHandlers[0]
erase LHCcentral:AnalysisHandlers[0]
insert LHCcentral:AnalysisHandlers[0] Logger
create DIPSY::HIAnalysis HIAna HIAnalysis.so
insert LHCcentral:AnalysisHandlers[0] HIAna
cp LHCcentral LHCppC
do LHCcentral:AddInterface /DIPSY/stdFiller:FSSwingTime 0, 0.5, 1.0, 1.5, 2.0
do LHCcentral:AddInterface /DIPSY/stdFiller:FSSwingTimeStep 0.05, 0.1
set LHCcentral:NumberOfEvents 1000
#saverun LHCppc LHCcentral
do LHCppC:AddInterface /DIPSY/stdFiller:FSSwingTime 0, 0.25, 0.5, 0.75, 1.0
do LHCppC:AddInterface /DIPSY/stdFiller:FSSwingTimeStep 0.05, 0.1, 0.2
set LHCppC:NumberOfEvents 1000
#saverun LHCppC LHCppC
cp LHCGenerator LHCridge0
set LHCridge0:NumberOfEvents 4000000
cp FixedB SmallB
set SmallB:MinB 0.0
set SmallB:MaxB 0.1
set LHCridge0:EventHandler:BGen SmallB
#erase LHCridge0:AnalysisHandlers[0]
erase LHCridge0:AnalysisHandlers[0]
erase LHCridge0:AnalysisHandlers[0]
insert LHCridge0:AnalysisHandlers[0] Logger
create ThePEG::CMSlrnsac CMSlrnsac CMSlrnsac.so
insert LHCridge0:AnalysisHandlers[0] CMSlrnsac
#saverun LHCridge0 LHCridge0
set SmallB:MaxB 0.2
#saverun LHCridge1 LHCridge0
set SmallB:MaxB 0.3
#saverun LHCridge2 LHCridge0
set stdEmitter:BothOrderedEvo False
set stdEmitter:BothOrderedInt False
set stdEmitter:BothOrderedFS False
set stdFiller:FSSwingTime 0.6
set stdFiller:FSSwingTimeStep 0.2
#saverun LHCridge3 LHCridge0
set SmallB:MaxB 0.2
#saverun LHCridge4 LHCridge0
set SmallB:MaxB 0.3
set SmallB:MinB 0.2
saverun LHCridge5 LHCridge0
set SmallB:MaxB 0.2
set SmallB:MinB 0.1
saverun LHCridge6 LHCridge0
set LHCcentral:EventHandler:BGen stdBGenerator
#saverun Test5ppFS70swing LHCGenerator
set LHCLumi:Energy 900.0
#saverun Test5ppFS09swing LHCGenerator
set LHCLumi:Energy 7000.0
set stdEmitter:BothOrderedEvo True
set stdEmitter:BothOrderedInt True
set stdEmitter:BothOrderedFS True
#saverun Test5ppFS70sw2 LHCGenerator
set LHCLumi:Energy 900.0
#saverun Test5ppFS09sw2 LHCGenerator
set LHCLumi:Energy 7000.0
set stdFiller:FSSwingTime 0.0
set stdFiller:FSSwingTimeStep 0.1

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:02 PM (23 h, 7 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242373
Default Alt Text
(142 KB)

Event Timeline