Page MenuHomeHEPForge
Files F10664348

No OneTemporary
Actions

View Options

diff --git a/DIPSY/PAXSec.in b/DIPSY/PAXSec.in
--- a/DIPSY/PAXSec.in
+++ b/DIPSY/PAXSec.in
@@ -1,457 +1,469 @@
cd /DIPSY
## First we setup some tuned parameters
#read CurrentTune.in
read Tune31.in
## Now we set up an event generator. We start with running pp rahter than pA.
cp EventHandler PAEventHandler
set PAEventHandler:WFL stdProton
set PAEventHandler:WFR stdProton
set PAEventHandler:ConsistencyLevel 0
set PAEventHandler:XSecFn:CheckOffShell false
set PAEventHandler:CascadeHandler NULL
set PAEventHandler:HadronizationHandler NULL
set PAEventHandler:DecayHandler NULL
create ThePEG::LuminosityFunction PALumi
set PAEventHandler:LuminosityFunction PALumi
cp Generator PAGenerator
set PAGenerator:EventHandler PAEventHandler
set PAGenerator:NumberOfEvents 0
set PAGenerator:EventHandler:EventFiller:PTCut 0.6
set PAEventHandler:BGen:Width 10
set PAEventHandler:EffectivePartonMode Colours
## These are the analysess we will run
## First the Glauber analyses
create DIPSY::GlauberAnalysis Glauber GlauberAnalysis.so
insert PAEventHandler:AnalysisHandlers[0] Glauber
## Some semi-inclusive cross section for DIPSY which need at least
## four combinations of left- and right-moving cascades.
create DIPSY::SemiInclusiveXSecAnalysis SemiIncl SemiInclusiveXSecAnalysis.so
insert PAEventHandler:AnalysisHandlers[0] SemiIncl
set PAEventHandler:PreSampleL 2
set PAEventHandler:PreSampleR 2
## This is just to keep track of the progress of a run
create DIPSY::AnalysisProgress AnaLog AnalysisProgress.so
set AnaLog:Interval 600
insert PAEventHandler:AnalysisHandlers[0] AnaLog
## Set the interaction frame
set PAEventHandler:YFrametest 0.5
# set PAEventHandler:YFrametest 0.9
## The sample rates need to be adjusted so that we get a reasonable
## statistics in a reasonable time. It is typically efficient to
## sample a number of impact parameter values for each pair of DIPSY
## cascades.
set PAEventHandler:PreSampleB 10
## But we need a good sample of cascades
set PAEventHandler:PreSamples 100000
## Set the pp energy we want to run with
set PALumi:BeamEMaxA 100
set PALumi:BeamEMaxB 100
## We run pp to get the nucleon--nucleon cross sections
saverun PA01pp0 PAGenerator
set PAEventHandler:BGen:Width 5
saverun PA01pp0t PAGenerator
set PAEventHandler:BGen:Width 10
saverun PA01pp0o PAGenerator
saverun PA01pp0c PAGenerator
# set PAEventHandler:EffectivePartonMode Relatives
set PALumi:BeamEMaxA 200
set PALumi:BeamEMaxB 200
saverun PA02pp0 PAGenerator
set PALumi:BeamEMaxA 500
set PALumi:BeamEMaxB 500
saverun PA05pp0 PAGenerator
set PALumi:BeamEMaxA 1000
set PALumi:BeamEMaxB 1000
saverun PA10pp0 PAGenerator
set PALumi:BeamEMaxA 2500
set PALumi:BeamEMaxB 2500
saverun PA25pp0 PAGenerator
set PALumi:BeamEMaxA 3500
set PALumi:BeamEMaxB 3500
saverun PA35pp0 PAGenerator
## Now we take these cross sections and feed them into the Glauber analysis
## My numbers for 100 GeV (in millibarns)
set Glauber:TotalnnXSec 49.71
set Glauber:ElasticnnXSec 9.10
set Glauber:InElasticnnXSec 35.80
set PALumi:BeamEMaxA 100
set PALumi:BeamEMaxB 100
## Run again to make sure that the cross sections are reproduced, at least by "grey3 disc"
saverun PA01pp1 PAGenerator
set PAEventHandler:PreSamples 1000000
saverun PA01pp2 PAGenerator
saverun PA01pp3 PAGenerator
## My numbers for 200 GeV
set Glauber:TotalnnXSec 58.52
set Glauber:ElasticnnXSec 11.64
set Glauber:InElasticnnXSec 41.05
set PALumi:BeamEMaxA 200
set PALumi:BeamEMaxB 200
set PAEventHandler:PreSamples 100000
saverun PA02pp1 PAGenerator
set PAEventHandler:PreSamples 1000000
saverun PA02pp2 PAGenerator
saverun PA02pp3 PAGenerator
## My numbers for 500 GeV
set Glauber:TotalnnXSec 70.53
set Glauber:ElasticnnXSec 15.36
set Glauber:InElasticnnXSec 47.88
set PALumi:BeamEMaxA 500
set PALumi:BeamEMaxB 500
set PAEventHandler:PreSamples 100000
saverun PA05pp1 PAGenerator
set PAEventHandler:PreSamples 1000000
saverun PA05pp2 PAGenerator
saverun PA05pp3 PAGenerator
## My numbers for 1000 GeV
set Glauber:TotalnnXSec 79.63
set Glauber:ElasticnnXSec 18.35
set Glauber:InElasticnnXSec 52.85
set PALumi:BeamEMaxA 1000
set PALumi:BeamEMaxB 1000
set PAEventHandler:PreSamples 100000
saverun PA10pp1 PAGenerator
set PAEventHandler:PreSamples 1000000
saverun PA10pp2 PAGenerator
saverun PA10pp3 PAGenerator
## My numbers for 2500 GeV
set Glauber:TotalnnXSec 91.43
set Glauber:ElasticnnXSec 22.35
set Glauber:InElasticnnXSec 59.12
set PALumi:BeamEMaxA 2500
set PALumi:BeamEMaxB 2500
set PAEventHandler:PreSamples 100000
saverun PA25pp1 PAGenerator
saverun PA25pp2 PAGenerator
set PAEventHandler:PreSamples 1000000
saverun PA25pp3 PAGenerator
set Glauber:TotalnnXSec 95.86
set Glauber:ElasticnnXSec 23.88
set Glauber:InElasticnnXSec 61.41
set PALumi:BeamEMaxA 3500
set PALumi:BeamEMaxB 3500
set PAEventHandler:PreSamples 100000
saverun PA35pp1 PAGenerator
saverun PA35pp2 PAGenerator
## Finally it's time to run with a heavy ion.
create DIPSY::GaussianImpactGenerator BGauss GaussianImpactGenerator.so
set PAEventHandler:BGen BGauss
set BGauss:Width 15
set PAEventHandler:WFR Oxygen
set PALumi:BeamEMaxA 100
set PALumi:BeamEMaxB 1600
set Glauber:TotalnnXSec 49.71
set Glauber:ElasticnnXSec 9.10
set Glauber:InElasticnnXSec 35.80
set Copper:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 1000
saverun PA01pO1 PAGenerator
set Copper:InterNucleonSwing Off
saverun PA01pO2 PAGenerator
set Oxygen:Rn 0
saverun PA01pO2a PAGenerator
set Oxygen:Rn -1.3
saverun PA01pO2b PAGenerator
set Oxygen:Rn 1.3
saverun PA01pO2c PAGenerator
set Oxygen:w 0.051
saverun PA01pO2d PAGenerator
set Oxygen:Rn 0
saverun PA01pO2e PAGenerator
set PAEventHandler:WFR Copper
# set PAEventHandler:WFR Lead
set PALumi:BeamEMaxA 100
set PALumi:BeamEMaxB 6300
set Glauber:TotalnnXSec 49.71
set Glauber:ElasticnnXSec 9.10
set Glauber:InElasticnnXSec 35.80
set Copper:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 1000
saverun PA01pCu1 PAGenerator
set Copper:InterNucleonSwing Off
saverun PA01pCu2 PAGenerator
set PALumi:BeamEMaxA 200
set PALumi:BeamEMaxB 12600
set Glauber:TotalnnXSec 58.52
set Glauber:ElasticnnXSec 11.64
set Glauber:InElasticnnXSec 41.05
set Copper:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 1000
saverun PA02pCu1 PAGenerator
set Copper:InterNucleonSwing Off
saverun PA02pCu2 PAGenerator
set PALumi:BeamEMaxA 500
set PALumi:BeamEMaxB 31500
set Glauber:TotalnnXSec 70.53
set Glauber:ElasticnnXSec 15.36
set Glauber:InElasticnnXSec 47.88
set Copper:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 1000
saverun PA05pCu1 PAGenerator
set Copper:InterNucleonSwing Off
saverun PA05pCu2 PAGenerator
set PALumi:BeamEMaxA 1000
set PALumi:BeamEMaxB 63000
set Glauber:TotalnnXSec 79.63
set Glauber:ElasticnnXSec 18.35
set Glauber:InElasticnnXSec 52.85
set Copper:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 500
saverun PA10pCu1 PAGenerator
set Copper:InterNucleonSwing Off
saverun PA10pCu2 PAGenerator
set PALumi:BeamEMaxA 2500
set PALumi:BeamEMaxB 157500
set Glauber:TotalnnXSec 91.43
set Glauber:ElasticnnXSec 22.35
set Glauber:InElasticnnXSec 59.12
set Copper:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 500
saverun PA25pCu1 PAGenerator
set Copper:InterNucleonSwing Off
saverun PA25pCu2 PAGenerator
# Let's do lead as well
set PAEventHandler:WFR Lead
# Just testing Glauber first
set Lead:Rn -0.9
set Lead:R 6.407
set Lead:a 0.459
set PALumi:BeamEMaxA 10
set PALumi:BeamEMaxB 2080
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 500
set Glauber:TotalnnXSec 49.71
set Glauber:ElasticnnXSec 9.10
set Glauber:InElasticnnXSec 35.80
saverun PA01pPbG PAGenerator
set Glauber:TotalnnXSec 58.52
set Glauber:ElasticnnXSec 11.64
set Glauber:InElasticnnXSec 41.05
saverun PA02pPbG PAGenerator
set Glauber:TotalnnXSec 70.53
set Glauber:ElasticnnXSec 15.36
set Glauber:InElasticnnXSec 47.88
saverun PA05pPbG PAGenerator
set Glauber:TotalnnXSec 79.63
set Glauber:ElasticnnXSec 18.35
set Glauber:InElasticnnXSec 52.85
saverun PA10pPbG PAGenerator
set Glauber:TotalnnXSec 91.43
set Glauber:ElasticnnXSec 22.35
set Glauber:InElasticnnXSec 59.12
saverun PA25pPbG PAGenerator
set Glauber:TotalnnXSec 95.6
set Glauber:ElasticnnXSec 23.82
set Glauber:InElasticnnXSec 61.22
saverun PA35pPbG PAGenerator
set Glauber:TotalnnXSec 95.86
set Glauber:ElasticnnXSec 23.88
set Glauber:InElasticnnXSec 61.41
saverun PA35pPbG2 PAGenerator
set Glauber:TotalnnXSec 95.79
set Glauber:ElasticnnXSec 23.86
set Glauber:InElasticnnXSec 61.36
saverun PA35pPbG3 PAGenerator
+set PAEventHandler:PreSamples 10
+set Lead:R 6.407
+set Lead:a 0.459
+set Lead:Rn 0.9
+set Lead:w 0
+set Lead:NTry 0
+saverun PA35pPb9GL0 PAGenerator
+set Lead:NTry 1
+saverun PA35pPb9GL1 PAGenerator
+set Lead:NTry 100
+saverun PA35pPb9GL100 PAGenerator
+
do Lead:SetNucleus Pb
set Lead:Rn 1.3
# Now include also DIPSY
set PALumi:BeamEMaxA 100
set PALumi:BeamEMaxB 20800
set Glauber:TotalnnXSec 49.71
set Glauber:ElasticnnXSec 9.10
set Glauber:InElasticnnXSec 35.80
set Lead:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 500
saverun PA01pPb1 PAGenerator
saverun PA01pPb1o PAGenerator
set Lead:Rn -1.3
saverun PA01pPb1c PAGenerator
set Lead:Rn 1.3
saverun PA01pPb1C PAGenerator
set Lead:Rn 0.0
saverun PA01pPb1h PAGenerator
set Lead:Rn -0.9
set Lead:R 6.407
set Lead:a 0.459
saverun PA01pPb1g PAGenerator
set PAEventHandler:PreSamples 50
saverun PA01pPb1g2 PAGenerator
set PAEventHandler:PreSamples 500
do Lead:SetNucleus Pb
set Lead:Rn 1.3
set Lead:InterNucleonSwing Off
saverun PA01pPb2 PAGenerator
set PALumi:BeamEMaxA 200
set PALumi:BeamEMaxB 41600
set Glauber:TotalnnXSec 58.52
set Glauber:ElasticnnXSec 11.64
set Glauber:InElasticnnXSec 41.05
set Lead:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 250
saverun PA02pPb1 PAGenerator
set Lead:InterNucleonSwing Off
saverun PA02pPb2 PAGenerator
set Lead:InterNucleonSwing On
set Lead:Rn -0.9
set Lead:R 6.407
set Lead:a 0.459
saverun PA02pPb3 PAGenerator
set Lead:InterNucleonSwing Off
saverun PA02pPb4 PAGenerator
set Lead:InterNucleonSwing On
do Lead:SetNucleus Pb
set Lead:Rn 1.3
set PALumi:BeamEMaxA 500
set PALumi:BeamEMaxB 10400
set Glauber:TotalnnXSec 70.53
set Glauber:ElasticnnXSec 15.36
set Glauber:InElasticnnXSec 47.88
set Lead:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 250
saverun PA05pPb1 PAGenerator
set Lead:InterNucleonSwing Off
saverun PA05pPb2 PAGenerator
set PALumi:BeamEMaxA 1000
set PALumi:BeamEMaxB 208000
set Glauber:TotalnnXSec 79.63
set Glauber:ElasticnnXSec 18.35
set Glauber:InElasticnnXSec 52.85
set Lead:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 250
saverun PA10pPb1 PAGenerator
set Lead:InterNucleonSwing Off
saverun PA10pPb2 PAGenerator
set PALumi:BeamEMaxA 2500
set PALumi:BeamEMaxB 416000
set Glauber:TotalnnXSec 91.43
set Glauber:ElasticnnXSec 22.35
set Glauber:InElasticnnXSec 59.12
set Lead:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 100
saverun PA25pPb1 PAGenerator
set Lead:InterNucleonSwing Off
saverun PA25pPb2 PAGenerator
set Lead:InterNucleonSwing On
set Lead:Rn -0.9
set Lead:R 6.407
set Lead:a 0.459
saverun PA25pPb3 PAGenerator
set Lead:InterNucleonSwing Off
saverun PA25pPb4 PAGenerator
set Lead:InterNucleonSwing On
do Lead:SetNucleus Pb
set Lead:Rn 1.3
set PALumi:BeamEMaxA 3500
set PALumi:BeamEMaxB 728000
set Glauber:TotalnnXSec 95.79
set Glauber:ElasticnnXSec 23.86
set Glauber:InElasticnnXSec 61.36
set Lead:InterNucleonSwing On
set PAEventHandler:PreSampleB 40
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 100
saverun PA35pPb1 PAGenerator
set Lead:InterNucleonSwing Off
saverun PA35pPb2 PAGenerator
set Lead:InterNucleonSwing On
set Lead:Rn -0.9
set Lead:R 6.407
set Lead:a 0.459
saverun PA35pPb3 PAGenerator
set Lead:InterNucleonSwing Off
saverun PA35pPb4 PAGenerator
set Lead:InterNucleonSwing On
do Lead:SetNucleus Pb
set Lead:Rn 1.3
# set PALumi:BeamEMaxB 12600
# set PALumi:BeamEMaxB 31500
# set PALumi:BeamEMaxB 63000
# set PALumi:BeamEMaxB 157500
## We may need to adjust the sampling. Especially since the nucleus
## cascade is time consuming, it may be worth running several protons
## per nucleus.
# set PAEventHandler:PreSampleB 10
set PAEventHandler:PreSampleL 10
set PAEventHandler:PreSamples 2000
## Let's go
# saverun PA01pCu2 PAGenerator
## And let's run one more time with the swing between nucleons turned
## off
set Copper:InterNucleonSwing Off
# saverun PA01pCu3 PAGenerator
set PAEventHandler:BGen:Width 10
# saverun PA01pCu4 PAGenerator
set PAEventHandler:BGen:Width 15
# saverun PA01pCu5 PAGenerator
set PAEventHandler:PreSamples 200
create DIPSY::GaussianImpactGenerator BGauss GaussianImpactGenerator.so
# saverun PA01pCu6 PAGenerator
set BGauss:Width 21
set PAEventHandler:BGen BGauss
# saverun PA01pCu7 PAGenerator
set BGauss:Width 15
set Copper:InterNucleonSwing On
set PAEventHandler:PreSamples 1000
# saverun PA01pCu80 PAGenerator
set PAEventHandler:PreSampleL 20
set PAEventHandler:PreSamples 500
# saverun PA01pCu81 PAGenerator
set PAEventHandler:PreSampleR 4
set PAEventHandler:PreSamples 250
# saverun PA01pCu82 PAGenerator
diff --git a/DIPSY/SimpleNucleus.cc b/DIPSY/SimpleNucleus.cc
--- a/DIPSY/SimpleNucleus.cc
+++ b/DIPSY/SimpleNucleus.cc
@@ -1,316 +1,351 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SimpleNucleus class.
//
#include "SimpleNucleus.h"
#include "SimpleNucleusState.h"
#include "DipoleEventHandler.h"
#include "NucleusData.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "ThePEG/PDT/StandardMatchers.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "gsl/gsl_sf_erf.h"
using namespace DIPSY;
SimpleNucleus::~SimpleNucleus() {}
void SimpleNucleus::initialize(const DipoleEventHandler & eh) {
if ( proton() ) proton()->initialize(eh);
if ( neutron() ) neutron()->initialize(eh);
if ( nucleon() ) nucleon()->initialize(eh);
}
Energy2 SimpleNucleus::m2() const {
if ( abs(particle()->id()) > 1000000000 ) return sqr(particle()->mass());
return sqr(A*particle()->mass());
}
DipoleStatePtr SimpleNucleus::
generate(const DipoleEventHandler & eh, Energy plus) {
vector<Point> positions;
- Length Rc = min(R - Rn/2.0, R);
while ( int(positions.size()) < A ) {
Point pos(2.0*UseRandom::rnd(-R, R),
2.0*UseRandom::rnd(-R, R),
2.0*UseRandom::rnd(-R, R));
Length r = pos.mag();
- if ( !UseRandom::rndbool((1.0 + w*sqr(r)/sqr(Rc))/
- (1.0 + exp((r - Rc)/a))) ) continue;
- bool overlap = false;
- // for ( int ir = 0; ir < 100; ir++) {
+ if ( !UseRandom::rndbool(ws(r)) ) continue;
+ int ntry = 0;
+ do {
+ bool overlap = false;
for ( int i = 0, N = positions.size(); i < N && !overlap; ++i )
- if ( (positions[i] - pos).mag() < abs(Rn) ) {
- overlap = true;
- break;
- }
- // if ( overlap ) {
- // // double theta = UseRandom::rnd( 0.,M_PI);
- // double theta = acos(UseRandom::rnd(-1.0,1.0));
- // double phi = UseRandom::rnd( 0.,2.*M_PI);
- // pos = Point(r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta));
- // overlap = false;
- // } else break;
- // }
+ if ( (positions[i] - pos).mag() < abs(Rn) ) overlap = true;
- if ( overlap ) continue;
+ if ( !overlap ) {
+ positions.push_back(pos);
+ break;
+ }
- positions.push_back(pos);
+ if ( nTry == 0 ) positions.clear();
+ else if ( ++ntry < nTry ) {
+ double theta = acos(UseRandom::rnd(-1.0,1.0));
+ double phi = UseRandom::rnd( 0.,2.*M_PI);
+ pos = Point(r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta));
+ }
+ } while ( ++ntry < nTry );
}
+ if ( rdist ) {
+ Point cms;
+ for ( int i = 0, N = positions.size(); i < N; ++i )
+ cms += positions[i];
+ cms /= double(positions.size());
+ for ( int i = 0, N = positions.size(); i < N; ++i ) {
+ double r = (positions[i] - cms).mag()/femtometer;
+ rdist->fill(r, sqr(1.0/r));
+ }
+ }
+
+
Energy minus = m2()/plus;
return new_ptr(SimpleNucleusState(eh, *this, plus, minus, positions,
new_ptr(WFInfo(this, R/Constants::hbarc))));
}
void SimpleNucleus::persistentOutput(PersistentOStream & os) const {
os << A << Z << ounit(R, femtometer) << ounit(a, femtometer)
- << w << wf << wfp << wfn << ounit(Rn, femtometer) << oenum(useInterNucleonSwing);
+ << w << wf << wfp << wfn << ounit(Rn, femtometer)
+ << oenum(useInterNucleonSwing) << nTry;
}
void SimpleNucleus::persistentInput(PersistentIStream & is, int) {
is >> A >> Z >> iunit(R, femtometer) >> iunit(a, femtometer)
- >> w >> wf >> wfp >> wfn >> iunit(Rn, femtometer) >> ienum(useInterNucleonSwing);
+ >> w >> wf >> wfp >> wfn >> iunit(Rn, femtometer)
+ >> ienum(useInterNucleonSwing) >> nTry;
}
string SimpleNucleus::setParameters(string cmd) {
cmd = StringUtils::stripws(cmd);
if ( Ptr<NucleusData>::pointer nucl = Repository::GetPtr<Ptr<NucleusData>::pointer>(cmd) ) {
A = nucl->A();
Z = nucl->Z();
setParticle(nucl);
} else {
if ( cmd == "He" ) {
A = 4;
Z = 2;
} else if ( cmd == "Li" ) {
A = 6;
Z = 3;
} else if ( cmd == "C" ) {
A = 12;
Z = 6;
} else if ( cmd == "O" ) {
A = 16;
Z = 8;
} else if ( cmd == "Al" ) {
A = 27;
Z = 13;
} else if ( cmd == "S" ) {
A = 32;
Z = 16;
} else if ( cmd == "Ca" ) {
A = 40;
Z = 20;
} else if ( cmd == "Ni" ) {
A = 58;
Z = 28;
} else if ( cmd == "Cu" ) {
A = 63;
Z = 29;
} else if ( cmd == "W" ) {
A = 186;
Z = 74;
} else if ( cmd == "Au" ) {
A = 197;
Z = 79;
} else if ( cmd == "Pb" ) {
A = 208;
Z = 82;
} else if ( cmd == "U" ) {
A = 238;
Z = 92;
} else {
return string("Nucleus \"") + cmd +
string("\" not available Please set individual parameters explicitly.");
}
}
switch ( abs(Z)*1000+abs(A) ) {
case 2004:
R = 1.71*femtometer;
a = 0.5*femtometer;
w = 0.0;
break;
case 3006:
R = 1.96*femtometer;
a = 0.5*femtometer;
w = 0.0;
break;
case 6012:
R = 2.47*femtometer;
a = 0.00001*femtometer;
w = 0.0;
break;
case 8016:
R = 2.608*femtometer;
a = 0.513*femtometer;
w = -0.051;
break;
case 13027:
R = 3.07*femtometer;
a = 0.519*femtometer;
w = 0.0;
break;
case 16032:
R = 3.458*femtometer;
a = 0.61*femtometer;
w = 0.0;
break;
case 20040:
R = 3.76*femtometer;
a = 0.586*femtometer;
w = -0.161;
break;
case 28058:
R = 4.309*femtometer;
a = 0.516*femtometer;
w = -0.1308;
break;
case 29063:
R = 4.2*femtometer;
a = 0.596*femtometer;
w = 0.0;
break;
case 74186:
R = 6.51*femtometer;
a = 0.535*femtometer;
w = 0.0;
break;
case 79197:
R = 6.38*femtometer;
a = 0.535*femtometer;
w = 0.0;
break;
case 82208:
R = 6.68*femtometer;
a = 0.546*femtometer;
w = 0.0;
break;
case 92238:
R = 6.68*femtometer;
a = 0.6*femtometer;
w = 0.0;
break;
default:
return string("Nucleus \"") + cmd +
string("\" not available. Please set individual parameters explicitly.");
}
return "";
}
+void SimpleNucleus::dofinish() {
+ WaveFunction::dofinish();
+}
+
+void SimpleNucleus::doinitrun() {
+ WaveFunction::doinitrun();
+ if ( HistFacPtr fac = generator()->histogramFactory() ) {
+ fac->initrun();
+ fac->registerClient(this);
+ fac->mkdirs("/DIPSY/SimpleNucleus");
+ rdist = fac->createHistogram1D("/DIPSY/SimpleNucleus/RDist",100,0.0,10.0);
+ }
+}
+
+
// Static variable needed for the type description system in ThePEG.
#include "ThePEG/Utilities/DescribeClass.h"
DescribeClass<SimpleNucleus,DIPSY::WaveFunction>
describeDIPSYSimpleNucleus("DIPSY::SimpleNucleus", "SimpleNucleus.so");
void SimpleNucleus::Init() {
static ClassDocumentation<SimpleNucleus> documentation
("The SimpleNucleus class represents the unevolved necleus wave function "
"described in terms of particles distributed according to a Wood-Saxon "
"potential.");
static Parameter<SimpleNucleus,int> interfaceA
("A",
"The mass number of the nucleus.",
&SimpleNucleus::A, 208, 0, 0,
true, false, Interface::lowerlim);
static Parameter<SimpleNucleus,int> interfaceZ
("Z",
"The charge of the nucleus.",
&SimpleNucleus::Z, 82, 0, 0,
true, false, Interface::lowerlim);
static Parameter<SimpleNucleus,Length> interfaceR
("R",
"The size of the nucleus in units of fermi.",
&SimpleNucleus::R, femtometer, 6.68*femtometer, 0.0*femtometer,
0.0*femtometer, true, false, Interface::lowerlim);
static Parameter<SimpleNucleus,Length> interfacea
("a",
"The suppression at the edge of the nucleus in units of fermi.",
&SimpleNucleus::a, femtometer, 0.546*femtometer, 0.0*femtometer,
0.0*femtometer, true, false, Interface::lowerlim);
static Parameter<SimpleNucleus,double> interfacew
("w",
"Optional quadratic term in Wood-Saxon potential",
&SimpleNucleus::w, 0.0, 0, 0,
true, false, Interface::nolimits);
static Reference<SimpleNucleus,WaveFunction> interfaceWF
("WF",
"The wave function describing individual nucleons.",
&SimpleNucleus::wf, true, false, true, true, false);
static Reference<SimpleNucleus,WaveFunction> interfaceWFP
("WFP",
"The wave function describing individual protons.",
&SimpleNucleus::wfp, true, false, true, true, false);
static Reference<SimpleNucleus,WaveFunction> interfaceWFN
("WFN",
"The wave function describing individual neutrons.",
&SimpleNucleus::wfn, true, false, true, true, false);
static Parameter<SimpleNucleus,Length> interfaceRn
("Rn",
"Size in units of fermi of an individual nucleon for the exclusion "
"mechanism.",
&SimpleNucleus::Rn, femtometer, 1.3*femtometer, 0.0*femtometer,
0.0*femtometer, true, false, Interface::nolimits);
static Command<SimpleNucleus> interfaceSetNucleus
("SetNucleus",
"Set parameters of the given nucleus according to measured values "
"if available.",
&SimpleNucleus::setParameters, true);
static Switch<SimpleNucleus,bool> interfaceInterNucleonSwing
("InterNucleonSwing",
"Allow for swings between the nucleons in the evolution.",
&SimpleNucleus::useInterNucleonSwing, true, true, false);
static SwitchOption interfaceInterNucleonSwingYes
(interfaceInterNucleonSwing,
"On",
"Allow for inter-nucleon swings.",
true);
static SwitchOption interfaceInterNucleonSwingNo
(interfaceInterNucleonSwing,
"Off",
"Disallow inter-nucleon swings.",
false);
+
+ static Parameter<SimpleNucleus,int> interfaceNTry
+ ("NTry",
+ "Number of times to try new angles for the same r-value if "
+ "overlapping nucleons are found. If zero, the whole set of "
+ "generated positions is discarded if anoverlap is found.",
+ &SimpleNucleus::nTry, 1, 0, 0,
+ true, false, Interface::lowerlim);
+
interfaceSetNucleus.rank(10);
interfaceA.rank(9);
interfaceZ.rank(8);
interfaceR.rank(7);
interfacea.rank(6);
interfaceWF.rank(5);
}
diff --git a/DIPSY/SimpleNucleus.h b/DIPSY/SimpleNucleus.h
--- a/DIPSY/SimpleNucleus.h
+++ b/DIPSY/SimpleNucleus.h
@@ -1,234 +1,277 @@
// -*- C++ -*-
#ifndef DIPSY_SimpleNucleus_H
#define DIPSY_SimpleNucleus_H
//
// This is the declaration of the SimpleNucleus class.
//
#include "Ariadne/DIPSY/WaveFunction.h"
+#include "ThePEG/Analysis/FactoryBase.h"
namespace DIPSY {
using namespace ThePEG;
/**
* The SimpleNucleus class represents the unevolved proton wave
* function described in terms of an equilatteral triangle of dipoles
* with a size distributed as a Gaussian.
*
* @see \ref SimpleNucleusInterfaces "The interfaces"
* defined for SimpleNucleus.
*/
class SimpleNucleus: public WaveFunction {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
SimpleNucleus()
: A(208), Z(82), R(6.68*femtometer), a(0.546*femtometer), w(0.0),
wf(WaveFunctionPtr()), wfp(WaveFunctionPtr()), wfn(WaveFunctionPtr()),
- Rn(1.3*femtometer), useInterNucleonSwing(true) {}
+ Rn(1.3*femtometer), useInterNucleonSwing(true), nTry(1), rdist(0) {}
/**
* The destructor.
*/
virtual ~SimpleNucleus();
//@}
public:
/** @name Main virtual functions to be overridden in sub-classes. */
//@{
/**
* Initialize the wavefunction for the given DipoleEventHandler.
*/
virtual void initialize(const DipoleEventHandler &);
/**
* Generate a dipole state according to this wavefunction, given an
* event handler and a positive light-cone momentum component.
*/
virtual DipoleStatePtr generate(const DipoleEventHandler & eh, Energy plus);
/**
* Return the invariant mass squared of the particle.
*/
virtual Energy2 m2() const;
//@}
- /** @name Simple access functions. */
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object. Called in the run phase just before
+ * a run begins.
+ */
+ virtual void doinitrun();
+
+ /**
+ * Finalize this object. Called in the run phase just after a
+ * run has ended. Used eg. to write out statistics.
+ */
+ virtual void dofinish();
+
+ //@}
+
+public:
+
+ /** @Name Simple access functions. */
//@{
/**
* Return the nucleon wave function.
*/
tWaveFunctionPtr nucleon() const {
return wf;
}
/**
* Return the proton wave function.
*/
tWaveFunctionPtr proton() const {
return wfp;
}
/**
* Return the number of protons.
*/
int nProtons() const {
return abs(Z);
}
/**
* Return the neutron wave function.
*/
tWaveFunctionPtr neutron() const {
return wfn;
}
/**
* Return the number of protons.
*/
int nNeutrons() const {
return abs(A) - nProtons();
}
/**
* Return the mass number.
*/
int massNumber() const {
return A;
}
+ /**
+ * Return true if we want inter-nucleon swings.
+ */
bool interNucleonSwing() const {
return useInterNucleonSwing;
}
+ /**
+ * Return the Wood-Saxon potential.
+ */
+ double ws(Length r) const {
+ Length Rc = min(R - Rn/2.0, R);
+ return (1.0 + w*sqr(r)/sqr(Rc))/(1.0 + exp((r - Rc)/a));
+ }
+
//@}
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 {
return new_ptr(*this);
}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const {
return new_ptr(*this);
}
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The mass number of the nucleaus.
*/
int A;
/**
* The charge of the nucleaus.
*/
int Z;
/**
* The size of the nucleus.
*/
Length R;
/**
* The suppression at the edge of the nucleus.
*/
Length a;
/**
* Optional quadratic term in Wood-Saxon potential
*/
double w;
/**
* Wave function of the individual nucleons.
*/
WaveFunctionPtr wf;
/**
* Wave function of the individual protons.
*/
WaveFunctionPtr wfp;
/**
* Wave function of the individual neutrons.
*/
WaveFunctionPtr wfn;
/**
* Size of an individual nucleon for the exclusion mechanism.
*/
Length Rn;
/**
* Allow swings between nucleons in the nuclei.
*/
bool useInterNucleonSwing;
/**
+ * Number of times to try new angles for the same r-value if
+ * overlapping nucleons are found. If zero, the whole set of
+ * generated positions is discarded if anoverlap is found.
+ */
+ int nTry;
+
+ /** The distribution of nucleon centers. */
+ FactoryBase::tH1DPtr rdist;
+
+
+ /**
* Helper function to set measured parametrs of a nucleus given by
* the name.
*/
string setParameters(string);
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SimpleNucleus & operator=(const SimpleNucleus &);
};
}
#endif /* DIPSY_SimpleNucleus_H */

File Metadata

Mime Type
text/x-diff
Expires
Thu, Apr 24, 6:37 AM (1 d, 19 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4878203
Default Alt Text
(30 KB)

Event Timeline