Page MenuHomeHEPForge
Files F10881835

No OneTemporary
Actions

View Options

diff --git a/DIPSY/GlauberAnalysis.cc b/DIPSY/GlauberAnalysis.cc
--- a/DIPSY/GlauberAnalysis.cc
+++ b/DIPSY/GlauberAnalysis.cc
@@ -1,1011 +1,1015 @@
// -*- 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 "Bin.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; }
valarray<double> vsqr(const valarray<double> & v) { return v*v; }
GlauberAnalysis::GlauberAnalysis()
: nnSigTot(50.0*millibarn), nnSigEl(9.2*millibarn), nnSigInND(35.9*millibarn)
, nnSigSD(2.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 = sigql = sigqr
= valarray<double>(ZERO, nstr);
sigtot2 = signd2 = sigel2 = sigdt2 = sigdr2 = sigdl2 = sigdd2 = sigql2 = sigqr2
= valarray<double>(ZERO, nstr);
hists = vector<FactoryBase::tH1DPtr>(8*nstr);
hdyl = vector<FactoryBase::tH1DPtr>(20);
hdyr = vector<FactoryBase::tH1DPtr>(20);
sumlr = sumnd = sum2lr = sumlr2 = suml2r = sumr2l =
Bin< valarray<double> >(valarray<double>(0.0, nstr));
bookHistos();
npartdiff = vector< vector<double> >(50, vector<double>(50, 0.0));
grey2R.sigTot = nnSigTot;
grey2R.sigEl = nnSigEl;
grey2R.sigD = nnSigTot - nnSigEl - nnSigInND;
grey2R.sigSD = nnSigSD;
grey2R.init();
GlauberBase::fit(grey2R);
stub(" grey2R Fitted sigTot:") << grey2R.stot() << " (" << nnSigTot/millibarn << ") mb" << endl;
stub(" grey2R Fitted sigEl:") << grey2R.sel() << " (" << nnSigEl/millibarn << ") mb" << endl;
stub(" grey2R Fitted sigD:") << grey2R.sD() << " (" << grey2R.sigD/millibarn << ") mb" << endl;
stub(" grey2R Fitted sigSD:") << grey2R.sSD() << " (" << nnSigSD/millibarn << ") mb" << endl;
stub(" grey2R Fitted c:") << grey2R.c << endl;
stub(" grey2R Fitted alpha:") << grey2R.alpha << endl;
stub(" grey2R Fitted R:") << grey2R.R*sqrt(millibarn)/femtometer << " fm" << endl;
stub(" grey2R Fitted r:") << grey2R.r*sqrt(millibarn)/femtometer << " fm" << endl;
}
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, 3*nstr);
T[0] = xsec.unitarize(fsum);
T[10] = 1.0;
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);
vector< vector<bool> > lpart(Nl, vector<bool>(nstr, false));
vector< vector<bool> > rpart(Nr, vector<bool>(nstr, false));
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;
tPartonPtr ppl = dl.initialDipoles()[il*3]->partons().first;
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;
tPartonPtr ppr = dr.initialDipoles()[ir*3]->partons().first;
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;
lpart[il][1] = rpart[ir][1] = true;
}
// 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;
lpart[il][2] = rpart[ir][2] = true;
}
// Gray3 disc
R = nnSigTot/sqrt(4.0*Constants::pi*nnSigEl);
a = nnSigEl/(nnSigTot - nnSigInND);
if ( r <= R && UseRandom::rndbool(a) ) {
++nalpha;
if ( UseRandom::rndbool(alpha*(2.0 - alpha)) )
lpart[il][3] = rpart[ir][3] = true;
}
// 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;
lpart[il][4] = rpart[ir][4] = true;
}
// Old black Disk
R = sqrt((nnSigTot - nnSigEl)/Constants::pi);
if ( r <= R ) {
T[5] = 1.0;
lpart[il][5] = rpart[ir][5] = true;
}
// Old gray disk
R = sqrt(nnSigTot/Constants::pi);
a = (nnSigTot - nnSigEl)/nnSigTot;
if ( r <= R && UseRandom::rndbool(a) ) {
T[6] = 1.0;
lpart[il][6] = rpart[ir][6] = true;
}
// ND black Disk
R = sqrt(nnSigInND/Constants::pi);
if ( r <= R ) {
T[7] = 1.0;
lpart[il][7] = rpart[ir][7] = true;
}
// Gray disc 2
R = nnSigTot/sqrt(4.0*Constants::pi*nnSigEl);
if ( r <= R && probr[&dr][ir] && probl[&dl][il] ) {
T[8] = 1.0;
lpart[il][8] = rpart[ir][8] = true;
}
// Gray disc strikman
R = nnSigTot/sqrt(4.0*Constants::pi*nnSigEl);
if ( r <= R && strikprob ) {
T[9] = 1.0;
lpart[il][9] = rpart[ir][9] = true;
}
// Gray disk with two radii
double TT = grey2R.T(ppr, ppl, r);
T[10] *= 1.0 - TT;
if ( UseRandom::rndbool(TT*(2.0 - TT)) )
lpart[il][10] = rpart[ir][10] = true;
// Gray disk with two radii but let right side fluctuate for each subcollision
TT = grey2R.T(ppl, r);
T[11] *= 1.0 - TT;
if ( UseRandom::rndbool(TT*(2.0 - TT)) )
lpart[il][11] = rpart[ir][11] = true;
// Gray disk with two radii but count also diffractively excited participants
TT = grey2R.T(ppr, ppl, r);
T[12] *= 1.0 - TT;
if ( UseRandom::rndbool(grey2R.Pr(ppl, r)) ) {
rpart[ir][12] = true;
if ( UseRandom::rndbool(grey2R.Par(ppl, r)/grey2R.Pr(ppl, r)) )
rpart[ir][13] = true;
}
if ( UseRandom::rndbool(grey2R.Pl(ppr, r)) ) {
lpart[il][12] = true;
if ( UseRandom::rndbool(grey2R.Pal(ppr, r)/grey2R.Pl(ppr, r)) )
lpart[il][13] = true;
}
}
}
T[3] = 1.0 - pow(1.0 - alpha, nalpha);
T[10] = 1.0 - T[10];
T[11] = 1.0 - T[11];
T[13] = T[12] = 1.0 - T[12];
for ( int istr = 1; istr < nstr; ++istr ) {
for ( int il = 0; il < Nl; ++il )
if ( lpart[il][istr] ) T[nstr + istr] += 1.0;
for ( int ir = 0; ir < Nr; ++ir )
if ( rpart[ir][istr] ) T[2*nstr + istr] += 1.0;
}
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) {
if ( a2 ) {
analyze2(vr, vl, vb, xsec, probs, jac);
return;
}
int Nr = vr.size();
int Nl = vl.size();
int Nb = vb.size();
if ( Nr*Nl*vb.size() == 0 ) return;
grey2R.clearStates();
vector < vector< vector < valarray<double> > > >
UTB(Nb, vector< vector < valarray<double> > >
(Nl, vector < valarray<double> >
(Nr, valarray<double>(nstr))));
for ( int ib = 0; 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 )
for ( int il1 = 0; il1 < Nl; ++il1 ) {
valarray<double> ret = getT(*vl[il1], *vr[ir1], vb[ib], xsec, probs[il1][ir1][ib]);
for ( int istr = 0; istr < nstr; ++istr )
UTB[ib][il1][ir1][istr] = UT[il1][ir1][istr] = ret[istr];
}
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;
sum2LR += sqr(UT1)*w1;
++n2LR;
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;
}
if ( il1 != il2 && ir1 == ir2 ) {
sumL2R += UT12*w1*w2;
++nL2R;
}
if ( il1 == il2 && ir1 != ir2 ) {
sumR2L += UT12*w1*w2;
++nR2L;
}
}
}
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();
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);
for ( int i = 0; i < nstr; ++i ) {
hists[i]->fill(sqr(q)/GeV2, J0*sumLR[i]*bw);
// hists[i + nstr]->fill(sqr(q)/GeV2, J0*sumLR[i]*bw);
// hists[i + 2*nstr]->fill(sqr(q)/GeV2, 2.0*Constants::pi*J0*sumLR[i]*bw);
hists[i + 3*nstr]->fill(sqr(q)/GeV2, J0*sqrt(sumL2R[i])*bw);
hists[i + 4*nstr]->fill(sqr(q)/GeV2, J0*sqrt(sumR2L[i])*bw);
}
}
++ntot;
}
double xBLR = double((Nb-1)*Nl*(Nl-1)*Nr*(Nr-1))/8.0;
vector< vector<double> > J0(Nb, vector<double>(Nt));
vector< vector< vector<double> > >
J02(Nb, vector< vector<double> >(Nb, vector<double>(Nt)));
for ( int it = 0; it < Nt; ++it ) {
Energy q = sqrt((double(it) + 0.5)*tMax/double(Nt));
for ( int ib1 = 0; ib1 < Nb; ++ib1 ) {
J0[ib1][it] = gsl_sf_bessel_J0(vb[ib1].bVec().pt()*q);
for ( int ib2 = ib1 + 1; ib2 < Nb; ++ib2 )
J02[ib1][ib2][it] =
gsl_sf_bessel_J0((vb[ib1].bVec() - vb[ib2].bVec()).pt()*q);
}
}
vector<double> bw(Nb);
for ( int ib1 = 0; ib1 < Nb; ++ib1 )
bw[ib1] = sqr(hbarc)*vb[ib1].weight()*jac/nanobarn;
for ( int ib1 = 0; ib1 < Nb; ++ib1 )
for ( int ir1 = 0; ir1 < Nr; ++ir1 )
for ( int il1 = 0; il1 < Nl; ++il1 ) {
valarray<double> UT1 = UTB[ib1][il1][ir1];
double w1 = vr[ir1]->weight()*vl[il1]->weight();
for ( int ib2 = ib1; ib2 < Nb; ++ib2 )
for ( int ir2 = ir1; ir2 < Nr; ++ir2 )
for ( int il2= il1; il2 < Nl; ++il2 ) {
valarray<double> UT2 = UTB[ib2][il2][ir2];
double w2 = vr[ir2]->weight()*vl[il2]->weight();
if ( ib1 != ib2 && il1 != il2 && ir1 != ir2 ) {
for ( int it = 0; it < Nt; ++it ) {
double q2 = (double(it) + 0.5)*tMax/double(Nt)/GeV2;
for ( int i = 0; i < nstr; ++i ) {
hists[i + nstr]->fill(q2,
J0[ib1][it]*UT1[i]*w1*bw[ib1]*
J0[ib2][it]*UT2[i]*w2*bw[ib2]/xBLR);
hists[i + 2*nstr]->fill(q2,
J02[ib1][ib2][it]*
UT1[i]*w1*bw[ib1]*
UT2[i]*w2*bw[ib2]/xBLR);
}
}
}
}
}
}
void GlauberAnalysis::
analyze2(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();
int Nb = vb.size();
if ( Nr*Nl*Nb == 0 ) return;
grey2R.clearStates();
vector < vector< vector < valarray<double> > > >
UTB(Nb, vector< vector < valarray<double> > >
(Nl, vector < valarray<double> >
(Nr, valarray<double>(nstr))));
for ( int ib1 = 0; ib1 < Nb; ++ib1 ) {
probl.clear();
probr.clear();
for ( int ir1 = 0; ir1 < Nr; ++ir1 )
for ( int il1 = 0; il1 < Nl; ++il1 ) {
valarray<double> ret = getT(*vl[il1], *vr[ir1], vb[ib1], xsec, probs[il1][ir1][ib1]);
double weight =
sqr(hbarc)*vb[ib1].weight()*vr[ir1]->weight()*vl[il1]->weight()
*(jac/nanobarn)/double(Nb*Nl*Nr);
int ia = int(ret[nstr + 13] + 0.5);
int id = int(ret[nstr + 12] + 0.5);
if ( ia >= 0 && id >= 0 && ia < 50 && id < 50 )
npartdiff[id][ia] += weight;
for ( int istr = 1; istr < nstr; ++istr ) {
UTB[ib1][il1][ir1][istr] = ret[istr];
hists[5*nstr + istr]->fill(ret[nstr + istr], weight);
hists[6*nstr + istr]->fill(ret[2*nstr + istr], weight);
hists[7*nstr + istr]->fill((ret[nstr + istr] + ret[2*nstr + istr])/2.0, weight);
}
pair< vector<GlauberAnalysis::Nucleon>, vector<GlauberAnalysis::Nucleon> >
wounded = dipsyParticipants(*vl[il1], *vr[ir1], vb[ib1], xsec, probs[il1][ir1][ib1]);
int nwl = 0;
for ( int iln = 0, Nln = wounded.first.size(); iln < Nln; ++iln )
if ( wounded.first[iln].participating ) ++nwl;
int nwr = 0;
for ( int irn = 0, Nrn = wounded.second.size(); irn < Nrn; ++irn )
if ( wounded.second[irn].participating ) ++nwr;
hists[5*nstr]->fill(nwl, weight);
hists[6*nstr]->fill(nwr, weight);
hists[7*nstr]->fill((nwl + nwr)/2.0, weight);
if ( nwl > 1 && nwl < 20 )
for ( int iln = 0, Nln = wounded.first.size(); iln < Nln; ++iln )
if ( wounded.first[iln].participating && wounded.first[iln].maxdeltay < 100 )
hdyl[nwl - 2]->fill(wounded.first[iln].maxdeltay, weight/(nwl - 1.0));
if ( nwr > 1 && nwr < 20 )
for ( int irn = 0, Nrn = wounded.second.size(); irn < Nrn; ++irn )
if ( wounded.second[irn].participating && wounded.second[irn].maxdeltay < 100 )
hdyr[nwr - 2]->fill(wounded.first[irn].maxdeltay, weight/(nwr - 1.0));
}
}
vector< vector<double> > J0(Nb, vector<double>(Nt));
vector< vector< vector<double> > >
J02(Nb, vector< vector<double> >(Nb, vector<double>(Nt)));
for ( int it = 0; it < Nt; ++it ) {
Energy q = sqrt((double(it) + 0.5)*tMax/double(Nt));
for ( int ib1 = 0; ib1 < Nb; ++ib1 ) {
J0[ib1][it] = gsl_sf_bessel_J0(vb[ib1].bVec().pt()*q);
for ( int ib2 = ib1 + 1; ib2 < Nb; ++ib2 )
J02[ib1][ib2][it] =
gsl_sf_bessel_J0((vb[ib1].bVec() - vb[ib2].bVec()).pt()*q);
}
}
vector<double> bw(Nb);
for ( int ib1 = 0; ib1 < Nb; ++ib1 )
bw[ib1] = sqr(hbarc)*vb[ib1].weight()*jac/nanobarn;
for ( int ib1 = 0; ib1 < Nb; ++ib1 )
for ( int ir1 = 0; ir1 < Nr; ++ir1 )
for ( int il1 = 0; il1 < Nl; ++il1 ) {
valarray<double> UT1 = UTB[ib1][il1][ir1];
double w1 = vr[ir1]->weight()*vl[il1]->weight();
sumlr += 2.0*UT1*w1*bw[ib1];
sumnd += (2.0*UT1 - vsqr(UT1))*w1*bw[ib1];
sum2lr += vsqr(UT1)*w1*bw[ib1];
for ( int it = 0; it < Nt; ++it )
tel0[it] += J0[ib1][it]*UT1*w1*bw[ib1];
for ( int ib2 = ib1; ib2 < Nb; ++ib2 )
for ( int ir2 = ir1; ir2 < Nr; ++ir2 )
for ( int il2 = il1; il2 < Nl; ++il2 ) {
valarray<double> UT2 = UTB[ib2][il2][ir2];
double w2 = vr[ir2]->weight()*vl[il2]->weight();
valarray<double> UT12 = UT1*UT2;
if ( ib1 == ib2 && il1 != il2 && ir1 != ir2 )
sumlr2 += UT12*w1*w2*bw[ib1];
if ( ib1 == ib2 && il1 != il2 && ir1 == ir2 )
suml2r += UT12*w1*w2*bw[ib1];
if ( ib1 == ib2 && il1 == il2 && ir1 != ir2 )
sumr2l += UT12*w1*w2*bw[ib1];
if ( ib1 != ib2 && il1 != il2 && ir1 != ir2 )
for ( int it = 0; it < Nt; ++it )
tel[it] += J02[ib1][ib2][it]*UT12*w1*w2*bw[ib1]*bw[ib2];
if ( ib1 != ib2 && il1 != il2 && ir1 == ir2 )
for ( int it = 0; it < Nt; ++it )
telql[it] += J02[ib1][ib2][it]*UT12*w1*w2*bw[ib1]*bw[ib2];
if ( ib1 != ib2 && il1 == il2 && ir1 != ir2 )
for ( int it = 0; it < Nt; ++it )
telqr[it] += J02[ib1][ib2][it]*UT12*w1*w2*bw[ib1]*bw[ib2];
}
}
}
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";
if ( i == 10 ) strat = "grey2R";
if ( i == 11 ) strat = "grey2R(fluctp)";
if ( i == 12 ) strat = "grey2R(diff)";
if ( i == 13 ) strat = "grey2R(diff/a)";
return strat;
}
void GlauberAnalysis::bookHistos() {
generator()->histogramFactory()->mkdirs("/Glauber");
generator()->histogramFactory()->mkdirs("/tmp");
for ( int i = 0; i < nstr; ++i ) {
string I(1, i < 10? '0' + i: 'a' + i - 10);
if ( a2 ) {
tel = tel0 = telql= telqr = vector< Bin< valarray<double> > >
(Nt, Bin< valarray<double> >(valarray<double>(0.0, nstr)));
hists[5*nstr + i] = generator()->histogramFactory()->createHistogram1D
("/tmp/dSigmadNl-" + I, Nt, -0.5, 199.5);
hists[6*nstr + i] = generator()->histogramFactory()->createHistogram1D
("/tmp/dSigmadNr-" + I, Nt, -0.5, 199.5);
hists[7*nstr + i] = generator()->histogramFactory()->createHistogram1D
("/tmp/dSigmadNlr-" + I, Nt, -0.5, 199.5);
} else {
hists[i] = generator()->histogramFactory()->createHistogram1D
("/tmp/dSigmadt-a-" + I, Nt, 0.0, tMax/GeV2);
hists[nstr + i] = generator()->histogramFactory()->createHistogram1D
("/Glauber/dSigmadt-b-" + I, Nt, 0.0, tMax/GeV2);
hists[2*nstr + i] = generator()->histogramFactory()->createHistogram1D
("/Glauber/dSigmadt-c-" + 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);
}
}
for ( int i = 0; i < 20; ++i ) {
string I(1, i + 2 < 10? '0' + i + 2: 'a' + i - 8);
hdyl[i] = generator()->histogramFactory()->createHistogram1D
("/tmp/deltayL" + I, 20, 0.0, 20.0);
hdyr[i] = generator()->histogramFactory()->createHistogram1D
("/tmp/deltayR" + I, 20, 0.0, 20.0);
}
}
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::print(valarray<double> sig, valarray<double> err,
string xstype) const {
for ( int i = 0; i < nstr; ++i ) {
string strat = getStrat(i);
stub(": " + xstype + ": (" + strat + "):")
<< ouniterr(sig[i], err[i], 1.0) << " nb." << endl;
}
}
void GlauberAnalysis::print(const Bin< valarray<double> > & sig,
string xstype) const {
for ( int i = 0; i < nstr; ++i ) {
string strat = getStrat(i);
stub(": " + xstype + ": (" + strat + "):")
<< ouniterr(sig.average()[i], sig.err()[i], 1.0) << " nb." << endl;
}
}
void GlauberAnalysis::finalize(long neve) {
if ( neve <= 0 ) return;
if ( a2 ) {
finalize2(neve);
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");
stub(": Diff. sing. exc. L+R: (pp):")
<< ouniterr(double(nnSigSD/nanobarn), 10000.0, 1.0)
<< " nb." << endl;
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, i < 10? '0' + i: 'a' + i - 10);
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)*4.0*Constants::pi));
generator()->histogramFactory()->histogramFactory().destroy(hists[i]);
stub(": Elastic cross check a: (" + getStrat(i) + ")")
<< tmp->sumAllBinHeights()*(tMax/GeV2)/double(Nt) << endl;
tmp = hists[i + nstr];
tmp->setTitle("dSigma/dt b (" + getStrat(i) + ")");
tmp->scale(nanobarn*GeV2/(sqr(hbarc)*double(ntot)*4.0*Constants::pi));
stub(": Elastic cross check b: (" + getStrat(i) + ")")
<< tmp->sumAllBinHeights()*(tMax/GeV2)/double(Nt) << endl;
tmp = hists[i + 2*nstr];
tmp->setTitle("dSigma/dt c (" + getStrat(i) + ")");
tmp->scale(nanobarn*GeV2/(sqr(hbarc)*double(ntot)*4.0*Constants::pi));
stub(": Elastic cross check c: (" + getStrat(i) + ")")
<< tmp->sumAllBinHeights()*(tMax/GeV2)/double(Nt) << endl;
tmp = generator()->histogramFactory()->histogramFactory().multiply
("/Glauber/dSigmadt-ql-" + I, *hists[3*nstr + i], *hists[3*nstr+ i]);
tmp->setTitle("dSigma/dt ql (" + getStrat(i) + ")");
tmp->scale(nanobarn*GeV2/(sqr(hbarc)*double(ntot*ntot)*2.0*Constants::pi));
generator()->histogramFactory()->histogramFactory().destroy(hists[3*nstr + i]);
tmp = generator()->histogramFactory()->histogramFactory().multiply
("/Glauber/dSigmadt-qr-" + I, *hists[4*nstr + i], *hists[4*nstr+ i]);
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[4*nstr + i]);
}
}
void GlauberAnalysis::finalize2(long neve) {
if ( neve <= 0 ) return;
stub(": Total: (pp) [New procedure]:")
<< ouniterr(double(nnSigTot/nanobarn), 10000.0, 1.0)
<< " nb." << endl;
print(sumlr, "Total");
stub(": Inelastic(ND): (pp):")
<< ouniterr(double(nnSigInND/nanobarn), 10000.0, 1.0)
<< " nb." << endl;
print(sumnd, "Inelastic(ND)");
stub(": Inelastic(tot): (pp):")
<< ouniterr(double((nnSigTot - nnSigEl)/nanobarn), 10000.0, 1.0)
<< " nb." << endl;
print(sumlr - sumlr2, "Inelastic(tot)");
stub(": Elastic: (pp):")
<< ouniterr(double(nnSigEl/nanobarn), 10000.0, 1.0) << " nb." << endl;
print(sumlr2, "Elastic");
stub(": Diff. total: (pp):")
<< ouniterr(double((nnSigTot - nnSigEl - nnSigInND)/nanobarn), 10000.0, 1.0)
<< " nb." << endl;
print(sum2lr - sumlr2, "Diff. total");
stub(": Diff. sing. exc. L+R: (pp):")
<< ouniterr(double(nnSigSD/nanobarn), 10000.0, 1.0)
<< " nb." << endl;
print(sumr2l - sumlr2, "Diff. exc. (R)");
print(suml2r - sumlr2, "Diff. exc. (L)");
print(sumr2l, "Quasi El. (R)");
print(suml2r, "Quasi El. (L)");
print(sum2lr - sumr2l - suml2r + sumlr2, "Double diff. exc.");
for (int i = 0; i < nstr; ++i )
stub(": <npart>L90 (" + getStrat(i) + "):")
<< avpart(hists[5*nstr + i], 0.9) << endl;
for (int i = 0; i < nstr; ++i )
stub(": <npart>R90 (" + getStrat(i) + "):")
<< avpart(hists[6*nstr + i], 0.9) << endl;
for (int i = 0; i < nstr; ++i )
stub(": <npart>L+R90 (" + getStrat(i) + "):")
<< avpart(hists[7*nstr + i], 0.9) << endl;
for (int i = 0; i < nstr; ++i )
stub(": <npart>L95 (" + getStrat(i) + "):")
<< avpart(hists[5*nstr + i], 0.95) << endl;
for (int i = 0; i < nstr; ++i )
stub(": <npart>R95 (" + getStrat(i) + "):")
<< avpart(hists[6*nstr + i], 0.95) << endl;
for (int i = 0; i < nstr; ++i )
stub(": <npart>L+R95 (" + getStrat(i) + "):")
<< avpart(hists[7*nstr + i], 0.95) << endl;
ofstream os((generator()->path() + "/" + generator()->runName() +
"-tdep.dat").c_str());
double tfac = nanobarn*GeV2/(sqr(hbarc)*4.0*Constants::pi);
for (int i = 0; i < nstr; ++i ) {
string I(1, i < 10? '0' + i: 'a' + i - 10);
os << "# tdep" << I << endl;
for ( int it = 0; it < Nt; ++it )
os << (double(it) + 0.5)*(tMax/GeV2)/double(Nt) << '\t'
<< tel[it].average()[i]*tfac << '\t' << tel[it].err()[i]*tfac << endl;
os << endl;
os << "# t0dep" << I << endl;
for ( int it = 0; it < Nt; ++it ) {
Bin< valarray<double> > b = tel0[it].sqr();
os << (double(it) + 0.5)*(tMax/GeV2)/double(Nt) << '\t'
<< b.average()[i]*tfac << '\t' << b.err()[i]*tfac << endl;
}
os << endl;
os << "# tdepr" << I << endl;
for ( int it = 0; it < Nt; ++it )
os << (double(it) + 0.5)*(tMax/GeV2)/double(Nt) << '\t'
<< telqr[it].average()[i]*tfac << '\t' << telqr[it].err()[i]*tfac << endl;
os << endl;
os << "# tdepl" << I << endl;
for ( int it = 0; it < Nt; ++it )
os << (double(it) + 0.5)*(tMax/GeV2)/double(Nt) << '\t'
<< telql[it].average()[i]*tfac << '\t' << telql[it].err()[i]*tfac << endl;
os << endl;
hists[5*nstr + i]->scale(1.0/neve);
hists[6*nstr + i]->scale(1.0/neve);
hists[7*nstr + i]->scale(1.0/neve);
}
for (int i = 0; i < 20; ++i ) {
if ( hdyl[i]->sumAllBinHeights() > 0.0 )
hdyl[i]->scale(1.0/hdyl[i]->sumAllBinHeights());
if ( hdyr[i]->sumAllBinHeights() > 0.0 )
hdyr[i]->scale(1.0/hdyr[i]->sumAllBinHeights());
}
ofstream grid((generator()->filename() + ".grid").c_str());
for ( int ia = 0, Na = 50; ia < Na; ++ia ) {
for ( int id = 0, Nd = 50; id < Nd; ++id )
grid << ia << '\t' << id << '\t' << npartdiff[id][ia]/neve << endl;
grid << endl;
}
}
double GlauberAnalysis::avpart(FactoryBase::tH1DPtr h, double frac) {
double tot = 0.0;
for ( int n = 1; n < 200; ++n ) tot += h->binHeight(h->coordToIndex(n));
double sw = 0.0;
double snwf = 0.0;
double swf = 0.0;
for ( int n = 1; n < 200; ++n ) {
double w = h->binHeight(h->coordToIndex(n));
sw += w;
if ( sw >= frac*tot) {
snwf += n*w;
swf += w;
}
}
return swf > 0? snwf/swf: 0.0;
}
void GlauberAnalysis::ParticipantParton::decayAll() {
- map<double, Parton *> children;
+ multimap<double, Parton *> children;
for ( set<tPartonPtr>::const_iterator cit = p->children().begin();
cit != p->children().end(); ++cit )
- children[(**cit).oY()] = *cit;
+ children.insert(make_pair((**cit).oY(), *cit));
decay(children);
}
void GlauberAnalysis::ParticipantParton::
-decay(map<double, Parton *> & children) {
+decay(multimap<double, Parton *> & children) {
if ( children.empty() ) {
final[p] = this;
return;
}
- map<double, Parton *>::iterator it = children.begin();
+ multimap<double, Parton *>::iterator it = children.begin();
double cy = it->first;
Parton * c = it->second;
children.erase(it);
+ next = create(*p, cy);
+ child = createInitial(*c, cy);
+ next->mother = this;
+ next->index = index;
if ( c->mainParent() == p ) {
- next = create(*p, cy);
- child = create(*c, cy);
- next->mother = child->mother = this;
- next->index = child->index = index;
- next->decay(children);
- child->decayAll();
- }
- else
- decay(children);
+ child->mother = this;
+ child->index = index;
+ } else {
+ child->father = this;
+ }
+ next->decay(children);
+ child->decayAll();
}
map<const Parton *, GlauberAnalysis::ParticipantParton *>
GlauberAnalysis::ParticipantParton::final;
+map<const Parton *, GlauberAnalysis::ParticipantParton *>
+GlauberAnalysis::ParticipantParton::initial;
vector<GlauberAnalysis::ParticipantParton *>
GlauberAnalysis::ParticipantParton::all;
pair< vector<GlauberAnalysis::Nucleon>, vector<GlauberAnalysis::Nucleon> >
GlauberAnalysis::
dipsyParticipants(const DipoleState & dr, const DipoleState & dl,
const ImpactParameters & b, const DipoleXSec & xsec,
double fsum) const {
pair< vector<Nucleon>, vector<Nucleon> > ret;
int Ndl = dl.initialDipoles().size();
if ( Ndl%3 && Ndl != 1 ) return ret;
int Ndr = dr.initialDipoles().size();
if ( Ndr%3 && Ndr != 1 ) return ret;
int Nl = dl.initialDipoles().size()/3;
if ( Ndl == 1 ) Nl = 1;
int Nr = dr.initialDipoles().size()/3;
if ( Ndr == 1 ) Nr = 1;
ParticipantParton::clear();
ret = make_pair(vector<Nucleon>(Nl), vector<Nucleon>(Nr));
vector<Nucleon> & LN = ret.first;
for ( int il = 0; il < Nl; ++il ) {
if ( Ndl > 1 ) {
LN[il].add(dl.initialDipoles()[il*3], il);
LN[il].add(dl.initialDipoles()[il*3 + 1], il);
LN[il].add(dl.initialDipoles()[il*3 + 2], il);
} else {
LN[il].add(dl.initialDipoles()[0], il);
}
}
vector<Nucleon> & RN = ret.second;
for ( int ir = 0; ir < Nr; ++ir ) {
if ( Ndr > 1 ) {
RN[ir].add(dr.initialDipoles()[ir*3], ir);
RN[ir].add(dr.initialDipoles()[ir*3 + 1], ir);
RN[ir].add(dr.initialDipoles()[ir*3 + 2], ir);
} else {
RN[ir].add(dr.initialDipoles()[0], ir);
}
}
DipoleXSec::FList fl = xsec.flist(dl, dr, b);
for ( DipoleXSec::FList::const_iterator it = fl.begin();
it != fl.end(); ++it ) {
if ( it->first.second < UseRandom::rnd() ) continue;
pair<bool,bool> psel =
xsec.int0Partons(it->second.first->partons().first,
it->second.first->partons().second,
it->second.second->partons().first,
it->second.second->partons().second, b);
ParticipantParton * pl =
ParticipantParton::getFinal(psel.first?
it->second.first->partons().first:
it->second.first->partons().second);
ParticipantParton * pl2 =
ParticipantParton::getFinal(!psel.first?
it->second.first->partons().first:
it->second.first->partons().second);
Nucleon & lnuc = LN[pl->index];
// double yl = -pl->interact();
double yl = -min(pl->interact(), pl2->interact());
ParticipantParton * pr =
ParticipantParton::getFinal(psel.second?
it->second.second->partons().first:
it->second.second->partons().second);
ParticipantParton * pr2 =
ParticipantParton::getFinal(!psel.second?
it->second.second->partons().first:
it->second.second->partons().second);
Nucleon & rnuc = RN[pr->index];
// double yr = -pr->interact();
double yr = -min(pr->interact(), pr2->interact());
if ( !lnuc.participating && !rnuc.participating )
rnuc.maxdeltay = lnuc.maxdeltay = 10000.0;
else {
rnuc.maxdeltay = max(yl - pr->origin()->y, rnuc.maxdeltay);
lnuc.maxdeltay = max(yr - pl->origin()->y, lnuc.maxdeltay);
}
rnuc.participating = lnuc.participating = true;
}
return ret;
}
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)
<< ounit(nnSigSD, nb) << ntot << sigtot << signd << sigel
<< sigdt << sigdl << sigdr << sigdd << sigql << sigqr
<< sigtot2 << signd2 << sigel2
<< 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)
>> iunit(nnSigSD, nb) >> ntot >> sigtot >> signd >> sigel
>> sigdt >> sigdl >> sigdr >> sigdd >> sigql >> sigqr
>> sigtot2 >> signd2 >> sigel2
>> sigdt2 >> sigdl2 >> sigdr2 >> sigdd2 >> sigql2 >> sigqr2
>> Nt >> iunit(tMax, GeV2);
hists.resize(8*nstr);
hdyl.resize(20);
hdyr.resize(20);
// 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,CrossSection> interfaceElasticnnXSecSD
("SDnnXSec",
"The assumed single-diffractive nucleon-nucleon cross section "
"assumed in the calculation (in mb).",
&GlauberAnalysis::nnSigSD, 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,333 +1,351 @@
// -*- 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>
#include "Bin.h"
#include "GlauberModels.h"
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:
struct ParticipantParton {
ParticipantParton():
- p(0), mother(0), next(0), child(0), y(-10000), interacted(false),
- index(-1) {}
+ p(0), mother(0), father(0), next(0), child(0),
+ y(-10000), interacted(false), index(-1) {}
const Parton * p; // The original parton;
ParticipantParton * mother; // the main parent (if any. ie. not initial).
+ ParticipantParton * father; // the main parent (if any. ie. not initial).
ParticipantParton * next; // The same parton after an emission.
ParticipantParton * child; // The emitted parton (if any, ie. has emitted).
double y; // The rapidity when this parton instance was produced
bool interacted; // true if this parton has already interacted
int index; // The index of the nucleon to which this parton belongs.
const ParticipantParton * origin() const {
return mother? mother->origin(): this;
}
double interact() {
interacted = true;
if ( !mother || mother->interacted ) return y;
- return mother->interact();
+ if ( !father || father->interacted ) return mother->interact();
+ return min(mother->interact(), father->interact());
}
static map<const Parton *, ParticipantParton *> final;
+ static map<const Parton *, ParticipantParton *> initial;
static vector<ParticipantParton *> all;
static void clear() {
for ( int i = 0, N = all.size(); i < N; ++i ) delete all[i];
all.clear();
final.clear();
+ initial.clear();
}
static ParticipantParton * create(const Parton & pin, double yin) {
ParticipantParton * p = new ParticipantParton();
p->p = &pin;
p->y = yin;
all.push_back(p);
return p;
}
+ static ParticipantParton * createInitial(const Parton & pin, double yin) {
+ if ( ParticipantParton * pold = getInitial(&pin) ) return pold;
+ ParticipantParton * pnew = create(pin, yin);
+ initial[&pin] = pnew;
+ return pnew;
+ }
+
static ParticipantParton * getFinal(const Parton * p) {
map<const Parton *, ParticipantParton *>::iterator pit = final.find(p);
if ( pit == final.end() ) return 0;
return pit->second;
}
+ static ParticipantParton * getInitial(const Parton * p) {
+ map<const Parton *, ParticipantParton *>::iterator pit = initial.find(p);
+ if ( pit == initial.end() ) return 0;
+ return pit->second;
+ }
+
+
void decayAll();
- void decay(map<double, Parton *> & children);
+ void decay(multimap<double, Parton *> & children);
};
struct Nucleon {
Nucleon(): maxdeltay(-10000.0), participating(false) {}
void add(const Dipole * d, int index) {
add(d->partons().first, index);
add(d->partons().second, index);
}
void add(const Parton * p, int index) {
if ( valence.find(p) == valence.end() ) {
- ParticipantParton * pp = ParticipantParton::create(*p, p->oY());
+ ParticipantParton * pp = ParticipantParton::createInitial(*p, p->oY());
valence[p] = pp;
pp->index = index;
pp->decayAll();
}
}
map<const Parton *, ParticipantParton *> valence;
double maxdeltay;
bool participating;
};
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);
virtual void analyze2(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);
virtual void finalize2(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;
/**
* The assumed single diffractive nucleon - nucleon cross section.
*/
CrossSection nnSigSD;
/**
* Number of strategies.
*/
static const int nstr = 14;
/**
* The sums.
*/
valarray<double> sigtot, signd, sigel,
sigdt, sigdr, sigdl, sigdd, sigql, sigqr;
Bin< valarray<double> > sumlr, sumnd, sum2lr, sumlr2, suml2r, sumr2l;
/**
* The sum of squares.
*/
valarray<double> sigtot2, signd2, sigel2,
sigdt2, sigdr2, sigdl2, sigdd2, sigql2, sigqr2;
/**
* 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, hdyl, hdyr;
vector< Bin< valarray<double> > > tel, tel0, telql, telqr;
vector< vector<double> > npartdiff;
/**
* 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;
pair< vector<Nucleon>, vector<Nucleon> >
dipsyParticipants(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;
void print(valarray<double> sig, valarray<double> sig2,
string xstype) const;
void print(const Bin< valarray<double> > & sig, string xstype) const;
/**
* helper function for histogram booking;
*/
void bookHistos();
/**
* Calculate average number ofparticipants above a certain fraction given a histogram.
*/
static double avpart(FactoryBase::tH1DPtr h, double frac);
/**
* helper function for printing
*/
string getStrat(int i) const;
static const bool a2 = 1;
Glauber2R grey2R;
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:52 AM (1 d, 3 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983171
Default Alt Text
(45 KB)

Event Timeline