Page MenuHomeHEPForge
Files F10881469

No OneTemporary
Actions

View Options

diff --git a/DIPSY/GlauberAnalysis.cc b/DIPSY/GlauberAnalysis.cc
--- a/DIPSY/GlauberAnalysis.cc
+++ b/DIPSY/GlauberAnalysis.cc
@@ -1,428 +1,455 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the GlauberAnalysis class.
//
#include "GlauberAnalysis.h"
#include "DipoleXSec.h"
#include "DipoleState.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "gsl/gsl_sf_bessel.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace DIPSY;
template <typename T>
valarray<T> sqr(valarray<T> v) { return v*v; }
GlauberAnalysis::GlauberAnalysis()
: nnSigTot(50.0*millibarn), nnSigEl(9.2*millibarn), nnSigInND(35.9*millibarn), Nt(1000), tMax(10.0*GeV2) {}
GlauberAnalysis::~GlauberAnalysis() {}
void GlauberAnalysis::initialize() {
generator()->histogramFactory()->initrun();
generator()->histogramFactory()->registerClient(this);
ntot = 0;
- sigtot = signd = sigel = sigdt = sigdr = sigdl = sigdd
+ sigtot = signd = sigel = sigdt = sigdr = sigdl = sigdd = sigql = sigqr
= valarray<double>(ZERO, nstr);
- sigtot2 = signd2 = sigel2 = sigdt2 = sigdr2 = sigdl2 = sigdd2
+ sigtot2 = signd2 = sigel2 = sigdt2 = sigdr2 = sigdl2 = sigdd2 = sigql2 = sigqr2
= valarray<double>(ZERO, nstr);
sTwLR = sTw2LR = sTwLR2 = sTwL2R = sTwR2L = valarray<double>(ZERO, nstr);
sT2wLR = sT2w2LR = sT2wLR2 = sT2wL2R = sT2wR2L = valarray<double>(ZERO, nstr);
swLR = sw2LR = swLR2 = swL2R = swR2L = 0.0;
- hists = vector<FactoryBase::tH1DPtr>(3*nstr);
+ hists = vector<FactoryBase::tH1DPtr>(5*nstr);
bookHistos();
}
void GlauberAnalysis::
setupProb(int Nr, const DipoleState & dr,
int Nl, const DipoleState & dl) const {
double ar = sqrt(2.0*nnSigEl/nnSigTot);
if ( probl.find(&dl) == probl.end() )
for ( int il = 0; il < Nl; ++il )
probl[&dl].push_back(UseRandom::rndbool(ar));
if ( probr.find(&dr) == probr.end() )
for ( int ir = 0; ir < Nr; ++ir )
probr[&dr].push_back(UseRandom::rndbool(ar));
}
valarray<double>
GlauberAnalysis::getT(const DipoleState & dr, const DipoleState & dl,
const ImpactParameters & b, const DipoleXSec & xsec,
double fsum) const {
valarray<double> T(0.0, nstr);
T[0] = xsec.unitarize(fsum);
int Ndl = dl.initialDipoles().size();
if ( Ndl%3 && Ndl != 1 ) return T;
int Ndr = dr.initialDipoles().size();
if ( Ndr%3 && Ndr != 1 ) return T;
int Nl = dl.initialDipoles().size()/3;
if ( Ndl == 1 ) Nl = 1;
int Nr = dr.initialDipoles().size()/3;
if ( Ndr == 1 ) Nr = 1;
setupProb(Nr, dr, Nl, dl);
double alpha = 2.0*(1.0 - nnSigInND/nnSigTot);
int nalpha = 0;
bool strikprob = UseRandom::rndbool(2.0*nnSigEl/nnSigTot);
for ( int il = 0; il < Nl; ++il ) {
Parton::Point pl =
(dl.initialDipoles()[0]->partons().first->position() +
dl.initialDipoles()[0]->partons().second->position())/2.0;
if ( Ndl > 1 )
pl = (dl.initialDipoles()[il*3]->partons().first->position() +
dl.initialDipoles()[il*3 + 1]->partons().first->position() +
dl.initialDipoles()[il*3 + 2]->partons().first->position())/3.0;
for ( int ir = 0; ir < Nr; ++ir ) {
Parton::Point pr =
(dr.initialDipoles()[0]->partons().first->position() +
dr.initialDipoles()[0]->partons().second->position())/2.0;
if ( Ndr > 1 )
pr = (dr.initialDipoles()[ir*3]->partons().first->position() +
dr.initialDipoles()[ir*3 + 1]->partons().first->position() +
dr.initialDipoles()[ir*3 + 2]->partons().first->position())/3.0;
Length r = b.difference(pl, pr).pt()*hbarc;
// Black disc
Length R = sqrt(nnSigTot/(2.0*Constants::pi));
if ( r <= R ) T[1] = 1.0;
// Gray disc
R = nnSigTot/sqrt(4.0*Constants::pi*nnSigEl);
double a = 2.0*nnSigEl/nnSigTot;
if ( r <= R && UseRandom::rndbool(a) ) T[2] = 1.0;
// Gray3 disc
R = nnSigTot/sqrt(4.0*Constants::pi*nnSigEl);
a = nnSigEl/(nnSigTot - nnSigInND);
if ( r <= R && UseRandom::rndbool(a) ) ++nalpha;
// Gussian
R = nnSigTot/sqrt(8.0*Constants::pi*nnSigEl);
a = 4.0*nnSigEl/nnSigTot;
if ( UseRandom::rndbool(a*exp(-sqr(r/R))) ) T[4] = 1.0;
// Old black Disk
R = sqrt((nnSigTot - nnSigEl)/Constants::pi);
if ( r <= R ) T[5] = 1.0;
// Old gray disk
R = sqrt(nnSigTot/Constants::pi);
a = (nnSigTot - nnSigEl)/nnSigTot;
if ( r <= R && UseRandom::rndbool(a) ) T[6] = 1.0;
// ND black Disk
R = sqrt(nnSigInND/Constants::pi);
if ( r <= R ) T[7] = 1.0;
// Gray disc 2
R = nnSigTot/sqrt(4.0*Constants::pi*nnSigEl);
if ( r <= R && probr[&dr][ir] && probl[&dl][il] ) T[8] = 1.0;
// Gray disc strikman
R = nnSigTot/sqrt(4.0*Constants::pi*nnSigEl);
if ( r <= R && strikprob ) T[9] = 1.0;
}
}
T[3] = 1.0 - pow(1.0 - alpha, nalpha);
return T;
}
void GlauberAnalysis::
analyze(const vector<DipoleStatePtr> & vr, const vector<DipoleStatePtr> & vl,
const vector<ImpactParameters> & vb, const DipoleXSec & xsec,
const Vec3D & probs, double jac) {
int Nr = vr.size();
int Nl = vl.size();
if ( Nr*Nl*vb.size() == 0 ) return;
for ( int ib = 0, Nb = vb.size(); ib < Nb; ++ib ) {
CrossSection bweight = sqr(hbarc)*vb[ib].weight()*jac;
double bw = bweight/nanobarn;
double bw2 = sqr(bw);
valarray<double> sumLR(0.0, nstr); // T averaged over L and R
int nLR = 0;
valarray<double> sum2LR(0.0, nstr); // T^2 average over L and R
int n2LR = 0;
valarray<double> sumLR2(0.0, nstr); // (T averaged over L and R )^2
int nLR2 = 0;
valarray<double> sumL2R(0.0, nstr); // (T averaged over L)^2 averaged over R
int nL2R = 0;
valarray<double> sumR2L(0.0, nstr); // (T averaged over R)^2 averaged over L
int nR2L = 0;
vector< vector < valarray<double> > >
UT(Nl, vector < valarray<double> >(Nr));
probl.clear();
probr.clear();
for ( int ir1 = 0; ir1 < Nr; ++ir1 )
for ( int il1 = 0; il1 < Nl; ++il1 )
UT[il1][ir1] = getT(*vl[il1], *vr[ir1], vb[ib], xsec, probs[il1][ir1][ib]);
for ( int ir1 = 0; ir1 < Nr; ++ir1 ) for ( int il1 = 0; il1 < Nl; ++il1 ) {
valarray<double> UT1 = UT[il1][ir1];
double w1 = vr[ir1]->weight()*vl[il1]->weight();
sumLR += UT1*w1;
++nLR;
sTwLR += UT1*w1*bw;
sT2wLR += sqr(UT1)*sqr(w1*bw);
swLR += 1.0;
// sum2LR += sqr(UT1)*sqr(w1);
sum2LR += sqr(UT1)*w1;
++n2LR;
sTw2LR += sqr(UT1)*w1*bw;
sT2w2LR += sqr(sqr(UT1))*sqr(w1*bw);
sw2LR += 1.0;
for ( int ir2 = 0; ir2 < Nr; ++ir2 ) for ( int il2 = 0; il2 < Nl; ++il2 ) {
valarray<double> UT2 = UT[il2][ir2];
double w2 = vr[ir2]->weight()*vl[il2]->weight();
// valarray<double> UT12 = UT1*w1*UT2*w2;
valarray<double> UT12 = UT1*UT2;
if ( il1 != il2 && ir1 != ir2 ){
sumLR2 += UT12*w1*w2;
++nLR2;
sTwLR2 += UT12*w1*w2*bw;
sT2wLR2 += sqr(UT12)*sqr(w1*w2*bw);
swLR2 += 1.0;
}
if ( il1 != il2 && ir1 == ir2 ) {
sumL2R += UT12*w1*w2;
++nL2R;
sTwL2R += UT12*w1*w2*bw;
sTwL2R += sqr(UT12)*sqr(w1*w2*bw);
swL2R += 1.0;
}
if ( il1 == il2 && ir1 != ir2 ) {
sumR2L += UT12*w1*w2;
++nR2L;
sTwR2L += UT12*w1*w2*bw;
sTwR2L += sqr(UT12)*sqr(w1*w2*bw);
swR2L += 1.0;
}
}
}
sumLR /=double(nLR);
sum2LR /=double(n2LR);
sumL2R /=double(nL2R);
sumR2L /=double(nR2L);
sumLR2 /=double(nLR2);
sigtot += 2.0*sumLR*bw;
sigtot2 += 4.0*sqr(sumLR)*bw2;
signd += (2.0*sumLR - sum2LR)*bw;
signd2 += (2.0*sumLR - sum2LR)*(2.0*sumLR - sum2LR)*bw2;
sigel += sumLR2*bw;
sigel2 += sumLR2*sumLR2*bw2;
sigdl += (sumL2R - sumLR2)*bw;
sigdl2 += (sumL2R - sumLR2)*(sumL2R - sumLR2)*bw2;
+ sigql += sumL2R*bw;
+ sigql2 += sumL2R*sumL2R*bw2;
sigdr += (sumR2L - sumLR2)*bw;
sigdr2 += (sumR2L - sumLR2)*(sumR2L - sumLR2)*bw2;
+ sigqr += sumR2L*bw;
+ sigqr2 += sumR2L*sumR2L*bw2;
sigdd += (sum2LR - sumR2L - sumL2R + sumLR2)*bw;
sigdd2 += (sum2LR - sumR2L - sumL2R + sumLR2)*(sum2LR - sumR2L - sumL2R + sumLR2)*bw2;
sigdt += (sum2LR - sumLR2)*bw;
sigdt2 += (sum2LR - sumLR2)*(sum2LR - sumLR2)*bw2;
InvEnergy b = vb[ib].bVec().pt();
double tw = vb[ib].weight()*jac*GeV2;
for ( int it = 0; it < Nt; ++it ) {
Energy q = sqrt((double(it) + 0.5)*tMax/double(Nt));
double J0 = gsl_sf_bessel_J0(b*q);
double sq = sqrt(q/GeV);
for ( int i = 0; i < nstr; ++i ) {
hists[i]->fill(sqr(q)/GeV2, sq*J0*sumLR[i]*bw);
if ( J0 > 0.0 )
hists[i + nstr]->fill(sqr(q)/GeV2, J0*sumLR2[i]*tw);
else
hists[i + 2*nstr]->fill(sqr(q)/GeV2, J0*sumLR2[i]*tw);
+ hists[i + 3*nstr]->fill(sqr(q)/GeV2, sq*J0*sumL2R[i]*bw);
+ hists[i + 4*nstr]->fill(sqr(q)/GeV2, sq*J0*sumR2L[i]*bw);
}
}
++ntot;
}
}
string GlauberAnalysis::getStrat(int i) const {
string strat = "DIPSY";
if ( i == 1 ) strat = "black disc";
if ( i == 2 ) strat = "grey disc";
if ( i == 3 ) strat = "grey3 disc";
if ( i == 4 ) strat = "Gaussian";
if ( i == 5 ) strat = "old black disc";
if ( i == 6 ) strat = "old grey disc";
if ( i == 7 ) strat = "blakc disc ND";
if ( i == 8 ) strat = "grey2 disc";
if ( i == 9 ) strat = "grey strick";
return strat;
}
void GlauberAnalysis::bookHistos() {
generator()->histogramFactory()->mkdirs("/Glauber");
generator()->histogramFactory()->mkdirs("/tmp");
for ( int i = 0; i < nstr; ++i ) {
string I(1, '0' + i);
hists[i] = generator()->histogramFactory()->createHistogram1D
("/tmp/dSigmadt-a-" + I, Nt, 0.0, tMax/GeV2);
hists[nstr + i] = generator()->histogramFactory()->createHistogram1D
("/tmp/dSigmadt-b-" + I, Nt, 0.0, tMax/GeV2);
hists[2*nstr + i] = generator()->histogramFactory()->createHistogram1D
("/tmp/-dSigmadt-b-" + I, Nt, 0.0, tMax/GeV2);
+ hists[3*nstr + i] = generator()->histogramFactory()->createHistogram1D
+ ("/tmp/dSigmadt-ql-" + I, Nt, 0.0, tMax/GeV2);
+ hists[4*nstr + i] = generator()->histogramFactory()->createHistogram1D
+ ("/tmp/dSigmadt-qr-" + I, Nt, 0.0, tMax/GeV2);
}
}
void GlauberAnalysis::print(valarray<double> sig, valarray<double> sig2,
int ntot, string xstype) const {
for ( int i = 0; i < nstr; ++i ) {
string strat = getStrat(i);
sig[i] /= double(ntot);
sig2[i] /= double(ntot);
double err = sqrt((sig2[i] - sqr(sig[i]))/double(ntot));
stub(": " + xstype + ": (" + strat + "):")
<< ouniterr(sig[i], err, 1.0) << " nb." << endl;
}
}
void GlauberAnalysis::finalize(long neve) {
if ( neve <= 0 ) return;
stub(": Total: (pp):")
<< ouniterr(double(nnSigTot/nanobarn), 10000.0, 1.0)
<< " nb." << endl;
print(sigtot, sigtot2, ntot, "Total");
stub(": Inelastic(ND): (pp):")
<< ouniterr(double(nnSigInND/nanobarn), 10000.0, 1.0)
<< " nb." << endl;
print(signd, signd2, ntot, "Inelastic(ND)");
stub(": Inelastic(tot): (pp):")
<< ouniterr(double((nnSigTot - nnSigEl)/nanobarn), 10000.0, 1.0)
<< " nb." << endl;
print(sigtot - sigel, sigtot2 - sigel2, ntot, "Inelastic(tot)");
stub(": Elastic: (pp):")
<< ouniterr(double(nnSigEl/nanobarn), 10000.0, 1.0) << " nb." << endl;
print(sigel, sigel2, ntot, "Elastic");
stub(": Diff. total: (pp):")
<< ouniterr(double((nnSigTot - nnSigEl - nnSigInND)/nanobarn), 10000.0, 1.0)
<< " nb." << endl;
print(sigdt, sigdt2, ntot, "Diff. total");
print(sigdr, sigdr2, ntot, "Diff. exc. (R)");
print(sigdl, sigdl2, ntot, "Diff. exc. (L)");
+ print(sigqr, sigqr2, ntot, "Quasi El. (R)");
+ print(sigql, sigql2, ntot, "Quasi El. (L)");
print(sigdd, sigdd2, ntot, "Double diff. exc.");
for (int i = 0; i < nstr; ++i ) {
string I(1, '0' + i);
FactoryBase::tH1DPtr tmp =
generator()->histogramFactory()->histogramFactory().multiply
("/Glauber/dSigmadt-a-" + I, *hists[i], *hists[i]);
tmp->setTitle("dSigma/dt a (" + getStrat(i) + ")");
tmp->scale(nanobarn*GeV2/(sqr(hbarc)*double(ntot*ntot)*2.0*Constants::pi));
generator()->histogramFactory()->histogramFactory().destroy(hists[i]);
tmp = generator()->histogramFactory()->histogramFactory().subtract
("/Glauber/dSigmadt-b-" + I, *hists[i + nstr], *hists[i + 2*nstr]);
tmp->setTitle("dSigma/dt b (" + getStrat(i) + ")");
tmp->scale(1.0/(double(ntot)*2.0*Constants::pi));
+ generator()->histogramFactory()->histogramFactory().destroy
+ (hists[i + nstr]);
+ generator()->histogramFactory()->histogramFactory().destroy
+ (hists[i + 2*nstr]);
+ tmp = generator()->histogramFactory()->histogramFactory().multiply
+ ("/Glauber/dSigmadt-ql-" + I, *hists[3*nstr + i], *hists[3*nstr+ 1]);
+ tmp->setTitle("dSigma/dt ql (" + getStrat(i) + ")");
+ tmp->scale(nanobarn*GeV2/(sqr(hbarc)*double(ntot*ntot)*2.0*Constants::pi));
+ tmp = generator()->histogramFactory()->histogramFactory().multiply
+ ("/Glauber/dSigmadt-qr-" + I, *hists[4*nstr + i], *hists[4*nstr+ 1]);
+ tmp->setTitle("dSigma/dt qr (" + getStrat(i) + ")");
+ tmp->scale(nanobarn*GeV2/(sqr(hbarc)*double(ntot*ntot)*2.0*Constants::pi));
+ generator()->histogramFactory()->histogramFactory().destroy(hists[i]);
}
}
IBPtr GlauberAnalysis::clone() const {
return new_ptr(*this);
}
IBPtr GlauberAnalysis::fullclone() const {
return new_ptr(*this);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void GlauberAnalysis::persistentOutput(PersistentOStream & os) const {
const CrossSection nb = nanobarn;
os << ounit(nnSigTot, nb) << ounit(nnSigEl, nb) << ounit(nnSigInND, nb)
<< ntot << sigtot << signd << sigel
- << sigdt << sigdl << sigdr << sigdd
+ << sigdt << sigdl << sigdr << sigdd << sigql << sigqr
<< sigtot2 << signd2 << sigel2
- << sigdt2 << sigdl2 << sigdr2 << sigdd2 << Nt << ounit(tMax,GeV2);
+ << sigdt2 << sigdl2 << sigdr2 << sigdd2 << sigql2 << sigqr2
+ << Nt << ounit(tMax,GeV2);
}
void GlauberAnalysis::persistentInput(PersistentIStream & is, int) {
const CrossSection nb = nanobarn;
is >> iunit(nnSigTot, nb) >> iunit(nnSigEl, nb) >> iunit(nnSigInND, nb)
>> ntot >> sigtot >> signd >> sigel
- >> sigdt >> sigdl >> sigdr >> sigdd
+ >> sigdt >> sigdl >> sigdr >> sigdd >> sigql >> sigqr
>> sigtot2 >> signd2 >> sigel2
- >> sigdt2 >> sigdl2 >> sigdr2 >> sigdd2 >> Nt >> iunit(tMax, GeV2);
+ >> sigdt2 >> sigdl2 >> sigdr2 >> sigdd2 >> sigql2 >> sigqr2
+ >> Nt >> iunit(tMax, GeV2);
hists.resize(3*nstr);
// bookHistos();
}
// Static variable needed for the type description system in ThePEG.
#include "ThePEG/Utilities/DescribeClass.h"
DescribeClass<GlauberAnalysis,DIPSY::DipoleAnalysisHandler>
describeDIPSYGlauberAnalysis("DIPSY::GlauberAnalysis",
"GlauberAnalysis.so");
void GlauberAnalysis::Init() {
static ClassDocumentation<GlauberAnalysis> documentation
("There is no documentation for the GlauberAnalysis class");
static Parameter<GlauberAnalysis,CrossSection> interfaceTotalnnXSec
("TotalnnXSec",
"The assumed total nucleon-nucleon cross section assumed in the calculation (in mb).",
&GlauberAnalysis::nnSigTot, millibarn, 96.0*millibarn, 0.0*millibarn, 0*millibarn,
true, false, Interface::lowerlim);
static Parameter<GlauberAnalysis,CrossSection> interfaceElasticnnXSec
("ElasticnnXSec",
"The assumed elastic nucleon-nucleon cross section assumed in the calculation (in mb).",
&GlauberAnalysis::nnSigEl, millibarn, 24.0*millibarn, 0.0*millibarn, 0*millibarn,
true, false, Interface::lowerlim);
static Parameter<GlauberAnalysis,CrossSection> interfaceElasticnnXSecInND
("InElasticnnXSec",
"The assumed inelastic, non-diffractive nucleon-nucleon cross section "
"assumed in the calculation (in mb).",
&GlauberAnalysis::nnSigInND, millibarn, 30.0*millibarn, 0.0*millibarn, 0*millibarn,
true, false, Interface::lowerlim);
static Parameter<GlauberAnalysis,int> interfaceNt
("Nt",
"Number of bins in t-plots.",
&GlauberAnalysis::Nt, 1000, 0, 0,
true, false, Interface::lowerlim);
static Parameter<GlauberAnalysis,Energy2> interfacetMax
("tMax",
"Maximum t-value in t-plot.",
&GlauberAnalysis::tMax, GeV2, 10.0*GeV2, 0.0*GeV2, 0*GeV2,
true, false, Interface::lowerlim);
}
diff --git a/DIPSY/GlauberAnalysis.h b/DIPSY/GlauberAnalysis.h
--- a/DIPSY/GlauberAnalysis.h
+++ b/DIPSY/GlauberAnalysis.h
@@ -1,221 +1,223 @@
// -*- C++ -*-
#ifndef DIPSY_GlauberAnalysis_H
#define DIPSY_GlauberAnalysis_H
//
// This is the declaration of the GlauberAnalysis class.
//
#include "DipoleAnalysisHandler.h"
#include "ThePEG/Analysis/FactoryBase.h"
#include "ThePEG/Analysis/LWH/AnalysisFactory.h"
#include <valarray>
namespace DIPSY {
using namespace ThePEG;
using std::valarray;
/**
* Here is the documentation of the GlauberAnalysis class.
*
* @see \ref GlauberAnalysisInterfaces "The interfaces"
* defined for GlauberAnalysis.
*/
class GlauberAnalysis: public DIPSY::DipoleAnalysisHandler {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
GlauberAnalysis();
/**
* The destructor.
*/
virtual ~GlauberAnalysis();
//@}
public:
/** @name Standard virtual functions inherited from the classes. */
//@{
/**
* Initialize the analysis object.
*/
virtual void initialize();
/**
* Analyze a given set of collision given a set of left- and
* right-moving dipole states in \a vl and \a vr, an a set of
* ImpactParameters object in \a vb, and a cross section object, \a
* xsec. Also a jacobian, \a jac, is supplied in case of varying
* total energy.
*/
virtual void analyze(const vector<DipoleStatePtr> & vr, const vector<DipoleStatePtr> & vl,
const vector<ImpactParameters> & vb, const DipoleXSec & xsec,
const Vec3D & probs, double jac);
/**
* Finalize the analysis, (compute statistics etc.). \a neve is the
* number of times analyze() has been called since last
* initialize().
*/
virtual void finalize(long neve);
//@}
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;
//@}
private:
/**
* Convenient typedef.
*/
typedef QTY<4,0,0>::Type CrossSection2;
/**
* The assumed total nucleon - nucleon cross section
*/
CrossSection nnSigTot;
/**
* The assumed elastic nucleon - nucleon cross section
*/
CrossSection nnSigEl;
/**
* The assumed inelastic, non-diffracitve nucleon - nucleon cross section
*/
CrossSection nnSigInND;
/**
* Number of strategies.
*/
static const int nstr = 10;
/**
* The sums.
*/
- valarray<double> sigtot, signd, sigel, sigdt, sigdr, sigdl, sigdd;
+ valarray<double> sigtot, signd, sigel,
+ sigdt, sigdr, sigdl, sigdd, sigql, sigqr;
valarray<double> sumlr, sum2lr, sumlr2, suml2r, sumr2l;
/**
* The sum of squares.
*/
- valarray<double> sigtot2, signd2, sigel2, sigdt2, sigdr2, sigdl2, sigdd2;
+ valarray<double> sigtot2, signd2, sigel2,
+ sigdt2, sigdr2, sigdl2, sigdd2,sigql2, sigqr2;
valarray<double> sTwLR, sTw2LR, sTwLR2, sTwL2R, sTwR2L;
valarray<double> sT2wLR, sT2w2LR, sT2wLR2, sT2wL2R, sT2wR2L;
double swLR, sw2LR, swLR2, swL2R, swR2L;
/**
* The number of b-values.
*/
long ntot;
/**
* The number of t-bins for the elastic cross section
*/
int Nt;
/**
* The maximum t-value for the elastic cross section.
*/
Energy2 tMax;
/**
* The t-values, first and second strategy
*/
vector<FactoryBase::tH1DPtr> hists;
/**
* The grey2 probabilities for the nucleons,
*/
mutable map<const DipoleState *, vector<bool> > probl, probr;
/**
* Calulate all transition probabilities
*/
valarray<double> getT(const DipoleState & dr, const DipoleState & dl,
const ImpactParameters & b, const DipoleXSec & xsec,
double fsum) const;
/**
* Setup probabilities for grey2 nucleons.
*/
void setupProb(int Nr, const DipoleState & dr,
int Nl, const DipoleState & dl) const;
/**
* helper function for printing
*/
void print(valarray<double> sig, valarray<double> sig2,
int ntot, string xstype) const;
/**
* helper function for histogram booking;
*/
void bookHistos();
/**
* helper function for printing
*/
string getStrat(int i) const;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
GlauberAnalysis & operator=(const GlauberAnalysis &);
};
}
#endif /* DIPSY_GlauberAnalysis_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:25 AM (1 d, 11 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983018
Default Alt Text
(21 KB)

Event Timeline