[Halld-npp] Primakoff cross section
Elton Smith
elton at jlab.org
Mon Mar 2 13:27:44 EST 2020
Hi Ilya,
Find attached the ROOT script that I used to produce all the plots. The
first plot shows the Crystal Ball data with my parameterization. You can
use the parameterization directly or take values of the points from the
plot if you run the script interactively. I will come by your office and
see if you need something more specific.
Cheers, Eton.
Elton Smith
Jefferson Lab MS 12H3
12000 Jefferson Ave STE 4
Newport News, VA 23606
(757)269-7625
(757)269-6331 fax
On 3/2/20 10:48 AM, Ilya Larin wrote:
> Thank you, Elton and Jose.
>
> The point I missed is that Mpipi doesn't have a spread in that case,
> unlike photoproduction
> on fixed target and defined by cm energy of colliding gammas.
> If Elton can send his cross-section points, I can redo my calculations
> for them or
> I can do it with the experimental points from that article. I will also
> double-check everything in my calculations today. It would be also
> interesting to compare
> with Elton's value obtained with form factors that I calculated last week.
>
> All the best,
> Ilya
>
> ------------------------------------------------------------------------
> *От:* Elton Smith <elton at jlab.org>
> *Отправлено:* 2 марта 2020 г. 7:48
> *Кому:* Ilya Larin <ilarin at jlab.org>
> *Копия:* Elton Smith <elton at jlab.org>; Halld-npp at jlab.org
> <Halld-npp at jlab.org>
> *Тема:* Re: Primakoff cross section
> Hi Ilya,
>
> The cross section sigma(gg -> pipi) is a two body reaction, so the
> energy dependence can be expressed as a function of the cm energy
> (Wpipi). It is the cross section as a function of energy but not
> differential in that variable. It reproduces the Crystal Ball Fig. 8
> data (see attachment). Note the dimensions of nb. From Eq. 4 of the
> proposal (see attached) you can check the dimensions, which is nb.
>
> Hope this helps. We need to agree on this cross section. Thanks, Elton.
>
> Elton Smith
> Jefferson Lab MS 12H3
> 12000 Jefferson Ave STE 4
> Newport News, VA 23606
> (757)269-7625
> (757)269-6331 fax
>
> On 3/1/20 6:06 AM, Ilya Larin wrote:
> >
> > Hi Elton,
> >
> > I can just comment that the 200 nb value I got has already a factor of
> > 1.25 taken into account.
> > I also do not quite understand (maybe because it's 6 am already) why
> > in Fig. 3 the units
> > are nb, not nb/GeV, maybe it is nb per histogram bin size, which is
> > 50MeV I assume, or something else. Having Fig. 3 understood is
> > essential for the cross-section results comparison,
> > because I used values from Jose and can/should put Fig. 3 values
> > instead to compare the results directly.
> >
> > Thank you,
> > Ilya
> >
> >
> > ------------------------------------------------------------------------
> > *От:* Elton Smith <elton at jlab.org>
> > *Отправлено:* 28 февраля 2020 г. 15:15
> > *Кому:* Elton Smith <elton at jlab.org>
> > *Копия:* Halld-npp at jlab.org <Halld-npp at jlab.org>; Ilya Larin
> > <ilarin at jlab.org>
> > *Тема:* Primakoff cross section
> > Hi Ilya,
> >
> > I compared your calculation of the cross section in the LOI. See figs.
> > 3-5. I get a cross section/nucleus of 345 nb. If you took Jose's cross
> > section directly, you need to divide your value 200 nb by 0.8 since it
> > is for |costhe|<0.8, which gives 250 nb, or 70% of my value.
> >
> > You can check my parameterization of sigma(gg->pi0pi0) (Fig. 3) with
> > Jose's. I assumed it is essentially flat above the threshold instead of
> > having a shape. This may explain most of the remaining discrepancy. We
> > may wish to update the document with your value since Jose's
> > parameterization is more constrained, but we can discuss this.
> >
> > Cheers, Elton.
> >
> > Elton Smith
> > Jefferson Lab MS 12H3
> > 12000 Jefferson Ave STE 4
> > Newport News, VA 23606
> > (757)269-7625
> > (757)269-6331 fax
> >
> > On 2/28/20 2:36 PM, Elton Smith wrote:
> > > Dear collaborators,
> > >
> > > I have posted the minutes from the NPP meeting this morning on the
> > > wiki at
> > >
> >
> https://halldweb.jlab.org/wiki/index.php/Pi0_Polarizability_Meeting_Feb_28,_2020#Minutes.
>
> >
> > >
> > >
> > > Ilya: You can go ahead an link your presentation directly to the
> > > meeting page.
> > >
> > > Thanks, Elton.
> > >
> >
>
-------------- next part --------------
Double_t sigma_ggpi0pi0_func (Double_t *x, Double_t *par){
// Parameterization for data from Masiske Phys Rev D 41 (1990) 3324.
// Returns cross section in units of nb
// constants
// Double_t const PI = 3.14159;
// parin[0] = A; // parameter 1: normalization
// parin[1] = mu; // parameter 2: mean of Fermi Function
// parin[2] = kT // parameter 3: transition parameter
// Double_t MPI =0.139570;
Double_t MPI =0.1349766;
Double_t PI=3.14159;
Double_t A = par[0];
Double_t mu = par[1];
Double_t kT = par[2];
Double_t Wpipi = x[0] ;
Double_t f;
if (Wpipi < 2*MPI) {
f = 0.5;
}
else {
f = 1 / (exp((Wpipi-mu)/kT) + 1); // Use Fermi function. Data approx flat out to about 0.8 GeV, ignore beyond that. Data for costhe<0.8
}
// cout << " Wpipi=" << Wpipi << " A=" << A << " mu=" << mu << " kT=" << kT << " f=" << f << endl;
return 2*A*(0.5-f)/0.8; // Convert to nb, assuming isotropic angular distribution of pions. Correct for 80% coverage
}
Double_t ff_func (Double_t *x, Double_t *par){
// return the nuclear form factor accourding to 3 parameter Fermi distribution (3pF). For the 2 parameter (2pF), w = 0.
// Parameters taken from H. De Vries, C.W. De Jager, and C. De Vries.
// Nuclear charge-density-distribution parameters from elastic electron scattering. Atomic Data and Nuclear Data Tables, 36(3):495 – 536, 1987
// For the Fermi Model
// See Journal of Research of the National Bureau of Standards - B Mathenatics and Mathematical Physics
// Vol. 70B, No. 1, Jan-Mar 1966. "The Form Factor of the Fermi Model Spatial Distribution," by Maximon and Schrack
//
// Function is a function of q, the three-momentum transfer to the nucleus.
// Input argument is -t
// Note that q is the 3-vector momentum, but for low -t, q ~ sqrt(-t).
// constants
Double_t alpha = 1/137.;
Double_t pi = 3.14159;
Double_t hbarc = 0.19733; // GeV*fm
Double_t q = sqrt(x[0])/hbarc; // take q to be in fm^-1. Input variable is positive (-t)
Double_t c0 = par[0]; // half-density radius
Double_t z0 = par[1]; // skin or diffuseness parameter
Double_t w = par[2]; // additional parameter for the 3pF parameterization.
Double_t rho0;
Double_t sum=0;
Int_t jmax=4;
for (Int_t j=1;j<jmax;j++) { // truncate after 3 terms, first term dominates.
Double_t sum1 =pow(-1.,j-1)*exp(-j*c0/z0)/(j*j*j);
sum += sum1;
// cout << "jmax=" << jmax << " j=" << j << " c0=" << c0 << " z0=" << z0 << " sum1=" << sum1 << " sum=" << sum << endl;
}
rho0 = (4*pi*c0/3)*( pi*z0*pi*z0 + c0*c0 + 8*pi*z0*z0*z0*sum);
rho0 = 1/rho0;
Double_t f = 0;
f = (4*pi*pi*rho0*z0*z0*z0)/(q*q*z0*z0*sinh(pi*q*z0)*sinh(pi*q*z0))
* (pi*q*z0 * cosh(pi*q*z0) * sin(q*c0) - q*c0 *cos(q*c0) * sinh(pi*q*z0) )
+ 8*pi*rho0*z0*z0*z0*sum;
// cout << " q=" << q << " f=" << f << endl;
if (q > 0) {
return f;
}
else {
cout << "*** ff_func, q=" << q << " return 1" << endl;
return 1; // FF is 1 at zero
}
}
Double_t ff_HO_func (Double_t *x, Double_t *par){
// return the nuclear form factor accourding to Harmonic Oscillator parameterization.
// Parameters taken from H. De Vries, C.W. De Jager, and C. De Vries.
// Nuclear charge-density-distribution parameters from elastic electron scattering. Atomic Data and Nuclear Data Tables, 36(3):495 – 536, 1987
//
// The Fourier transfrom is obtained as Transform[exp(-(x/a)^2) + alpha(x/a)^2 exp(-(x/a)^2)] = (a/2sqrt(2) exp(-a^2q^2/4) + alpha*a/2sqrt(2) exp(-a^2q^2/4) (1-aq)
//
// Function is a function of q, the three-momentum transfer to the nucleus.
// Input argument is -t
// Note that q is the 3-vector momentum, but for low -t, q ~ sqrt(-t).
// constants
Double_t pi = 3.14159;
Double_t hbarc = 0.19733; // GeV*fm
Double_t q = sqrt(x[0])/hbarc; // take q to be in fm^-1. Input variable is positive (-t)
Double_t rHO = par[0]; // proton: mean-square radius
Double_t aHO = par[1]; // HO scaling factor a
Double_t alphaHO = par[2]; // parameter unconstrained for MHO
Double_t ff0 = 1/(1+alphaHO); // form factor nomalization - ignore factor of a/2sqrt(2)
Double_t f = ff0*exp(-aHO*aHO*q*q/4.) * (1 + alphaHO*(1-aHO*aHO*q*q/2.)); // Fourier transform of Harmonic Oscillator charge density
cout << "-t=" << x[0] << " q=" << q << " f=" << f << " rHO=" << rHO << " aHO=" << aHO << " alphaHO=" << alphaHO << " ff0=" << ff0 << endl;
if (q > 0) {
return f;
}
else {
cout << "*** ff_HO_func, q=" << q << " return 1" << endl;
return 1;
}
}
Double_t sigmat_func (Double_t *x, Double_t *par){
// return the cross section for Primakoff production of pi+pi-. CPP proposal PR12-13-008 Eq 8.
// independent variable is momentum transfer -t, in GeV2;
// constants
Double_t alpha = 1/137.;
Double_t pi = 3.14159;
Double_t betapipi = 0.999; // beta=0.999 (W=0.3), beta=0.997 (W=0.4), beta= 0.986 (W=1.0 GeV)
// Double_t sigmagg = 1; // take sigma (gg -> pi+ pi-) = 1 for now
// Double_t Z = 50; // Z of Sn, target
Double_t Z = 82; // Z of Pb, target
Double_t coef = 4*alpha*Z*Z/(pi);
Double_t Wpipi = par[0];
Double_t Eg = par[1];
Double_t t = -x[0] ; // theta of pipi system in the lab. Input Variable is positive (i.e. -t)
Double_t xin[2];
Double_t parin[3];
xin[0] = -t; // input variable to ff is positive (-t)
parin[0] = par[2]; // half-density radius, fm
parin[1] = par[3]; // diffuseness paramter, fm
parin[2] = par[4]; // exptra parameter for 3pF
Double_t FF = ff_func(xin,parin);
// include other masses here
Double_t m1 = 0; // mass of photon, incident beam
// Double_t m2 = 108.*0.931494; // mass of 116Sn, target. Conversion to MeV
Double_t m2 = 208.*0.931494; // mass of 208Pb, target, Conversion to MeV
Double_t m3 = Wpipi; // mass of 2pi system, scattered system
Double_t m4 = m2; // recoil target
xin[0] = Wpipi; // W, 2pi mass
parin[0] = par[5]; // parameter 1: A normalization
parin[1] = par[6]; // parameter 1: mu threshold
parin[2] = par[7]; // parameter 1: kT transition
Double_t sigmagg = sigma_ggpi0pi0_func(xin,parin);
// cout << " Wpipi=" << Wpipi << " par0=" << parin[0] << " par1=" << parin[1] << " par2=" << parin[2] << " sig=" << sigmagg << endl;
Double_t f = 0;
Double_t s = m2*m2 + 2*Eg*m2;
Double_t sqrts = s > 0? sqrt(s) : 0;
if (s < 0) {
cout << "*** sigma_func: s =" << s << " < 0!" << endl;
return f;
}
Double_t E1cm = (s+m1*m1-m2*m2)/(2*sqrts);
Double_t E2cm = (s+m2*m2-m1*m1)/(2*sqrts);
Double_t E3cm = (s+m3*m3-m4*m4)/(2*sqrts);
Double_t E4cm = (s+m4*m4-m3*m3)/(2*sqrts);
Double_t p1cm = E1cm*E1cm - m1*m1? sqrt(E1cm*E1cm - m1*m1) : 0;
Double_t p2cm = E2cm*E2cm - m2*m2? sqrt(E2cm*E2cm - m2*m2) : 0;
Double_t p3cm = E3cm*E3cm - m3*m3? sqrt(E3cm*E3cm - m3*m3) : 0;
Double_t p4cm = E4cm*E4cm - m4*m4? sqrt(E4cm*E4cm - m4*m4) : 0;
Double_t arg = (m1*m1-m3*m3-m2*m2+m4*m4)/(2*sqrts);
Double_t t0 = arg*arg - (p1cm - p3cm)*(p1cm - p3cm);
Double_t t1 = arg*arg - (p1cm + p3cm)*(p1cm + p3cm);
Double_t betastar = Eg/(Eg + m2);
Double_t gammastar = (Eg + m2)/sqrts;
Double_t betapipicm = p3cm/E3cm;
Double_t thepipicm = t0 -t > 0? sqrt((t0 -t)/(p1cm*p3cm)) : 0;
Double_t conv = 1./(gammastar*(1 + betastar/betapipicm));
if (-t > -t0) {
f = (coef/2)* (sigmagg/Wpipi)*Eg*Eg*Eg*Eg * (t0-t)* betapipi*betapipi * (FF*FF/(t*t))*conv*conv*conv*conv/(p1cm*p3cm*p1cm*p3cm);
}
else {
f = 0;
}
// cout << " t=" << t << " betastar=" << betastar << " gammastar=" << gammastar << " betapipicm=" << betapipicm << " t0=" << t0 << " thepipicm=" << thepipicm << " f=" << f << endl;
// cout << " t=" << t << " FF=" << FF << " f=" << f << endl;
return f;
}
void sigma_2pi0(void)
{
// File: sigma_2pi0.C
// Compute the cross section for Primakoff production of pi+pi- in gamma A -> A pi0 pi0 reaction. Copied from sigma_2pi and changed sigma_gg
// See Eq. 8 from CPP proposal PR12-13-008
//
gStyle->SetPalette(1,0);
gStyle->SetOptStat(111111);
gStyle->SetOptFit(111111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
char string[256];
TString filename = "sigma_2pi0_figs";
// define function
// Define historgrams
Int_t nbins=20;
// define sigma_2pi distribution curve to plot
Int_t npar = 8;
Double_t Wpipi = 0.32;
// Double_t MPI =0.139570;
Double_t MPI =0.1349766;
Double_t Eg = 6;
Double_t c0 = 6.62; // Pb half-density radius, fm
Double_t z0 = 0.546; // Pb difuseness parameter, fm
Double_t w0 = 0; // extra parameter for 3pF
// Double_t R0 = 5.358; // Sn half-density radius, fm
// Double_t a0 = 0.550; // Sn difuseness parameter, fm
Double_t A= 9.6; // fit to data for costhe<0.8 Marsiske Phys Rev D 41 (1990) 3324
Double_t mu= 2*MPI; // should go to zero at threshold
Double_t kT= 0.028; // transition over about 100 MeV.
Double_t xmin=0;
Double_t xmax=1;
Double_t tmin=0; // use -t (positive) as input
Double_t tmax=0.005;
TF1 *sigmat = new TF1 ("sigma",sigmat_func,tmin,tmax,npar);
sigmat->SetParameters(Wpipi,Eg,c0,z0,w0,A,mu,kT);
sigmat->SetParNames("Wpipi","Eg","c0","z0","w0","A","mu","kT");
tmin=0; // use -t (positive) as input
tmax=0.005;
npar = 3;
// tmax=0.2;
TF1 *ff = new TF1("ff",ff_func,tmin,tmax,npar);
ff->SetParameters(c0,z0,w0);
ff->SetParNames("c0","z0","w0");
Double_t aHO=2.50; // fm. 9Be
Double_t alphaHO =1.77; // fm. 9Be
Double_t rHO = 0.631; // fm. 9Be
npar = 3;
TF1 *ff_HO = new TF1("ff_HO",ff_HO_func,tmin,tmax,npar); // this is the form factor, but same inputs as charge distribution.
ff_HO->SetParameters(rHO,aHO,alphaHO);
ff_HO->SetParNames("rHO","aHO","alphaHO");
npar = 3;
Double_t Wmin=0.;
Double_t Wmax=0.8;
TF1 *sigma_ggpi0pi0 = new TF1("sigma_ggpi0pi0",sigma_ggpi0pi0_func,Wmin,Wmax,npar);
sigma_ggpi0pi0->SetParameters(A,mu,kT);
sigma_ggpi0pi0->SetParNames("A","mu","kT");
const Int_t npts=9;
Double_t yMarsiske[npts]={60,60,130,150,140,155,135,190,150};
Double_t yMarsiske_error[npts]={30,25,50,50,50,50,50,90,75}; // full error
Double_t xMarsiske[npts]={0.275,0.325,0.375,0.425,0.475,0.525,0.575,0.625,0.675};
Double_t xMarsiske_error[npts]={0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05}; // bin width
Double_t scale = 15./237./0.8; // 15 nb is 237 on ruler scale.0.8 factor to scale to full coverage
for (Int_t jj=0; jj<npts; jj++) {
yMarsiske[jj] *= scale; //
yMarsiske_error[jj] *= scale/2; // divide full error by 2
xMarsiske_error[jj] /= 2; // divide bin width by 2
}
TGraphErrors *gr_Marsiske = new TGraphErrors(npts,xMarsiske,yMarsiske,xMarsiske_error,yMarsiske_error);
TCanvas *c1 = new TCanvas("c1", "c1",200,10,1000,700);
gPad->SetGridx();
gPad->SetGridy();
Double_t ymin=0.;
Double_t ymax=20;
sigma_ggpi0pi0->SetTitle("#gamma #gamma #rightarrow #pi^{0} #pi^{0}");
// sigma_ggpi0pi0->GetXaxis()->SetRangeUser(xmin,xmax);
sigma_ggpi0pi0->GetYaxis()->SetRangeUser(ymin,ymax);
sigma_ggpi0pi0->GetXaxis()->SetTitleSize(0.05);
sigma_ggpi0pi0->GetYaxis()->SetTitleSize(0.05);
sigma_ggpi0pi0->GetXaxis()->SetTitle("M_{#pi#pi} (GeV)");
sigma_ggpi0pi0->GetYaxis()->SetTitle("#sigma (nb)");
sigma_ggpi0pi0->SetLineColor(4);
// sigma_ggpi0pi0->FixParameter(1,mu);
sigma_ggpi0pi0->Draw();
// gr_Marsiske->Fit("sigma_ggpi0pi0");
gr_Marsiske->SetMarkerStyle(20);
gr_Marsiske->SetMarkerSize(1.0);
gr_Marsiske->SetMarkerColor(2);
gr_Marsiske->Draw("psame");
TLegend *leg0 = new TLegend(0.2,0.7,0.45,0.85);
leg0->AddEntry(gr_Marsiske,"Crystal Ball 1990 / 0.8","p");
leg0->AddEntry(sigma_ggpi0pi0,"Parameterization","l");
leg0->Draw();
TCanvas *c2 = new TCanvas("c2", "c2",200,10,1000,700);
// c1->cd(2);
gPad->SetLogy();
sigmat->SetTitle("");
// sigmat->GetXaxis()->SetRangeUser(xmin,xmax);
// sigmat->GetYaxis()->SetRangeUser(ymin,ymax);
sigmat->GetXaxis()->SetTitleSize(0.05);
sigmat->GetYaxis()->SetTitleSize(0.05);
sigmat->GetXaxis()->SetTitle("-t (GeV^{2})");
sigmat->GetYaxis()->SetTitle("#sigma");
sigmat->SetLineColor(4);
sigmat->Draw();
TF1 *sigmat1 = sigmat->DrawCopy("same");
sigmat1->SetLineColor(4);
Wpipi=0.32;
sigmat1->SetParameters(Wpipi,Eg,c0,z0,w0,A,mu,kT);
TF1 *sigmat2 = sigmat->DrawCopy("same");
sigmat2->SetLineColor(2);
Wpipi=0.39;
sigmat2->SetParameters(Wpipi,Eg,c0,z0,w0,A,mu,kT);
TF1 *sigmat3 = sigmat->DrawCopy("same");
sigmat3->SetLineColor(1);
Wpipi=0.47;
sigmat3->SetParameters(Wpipi,Eg,c0,z0,w0,A,mu,kT);
Double_t int1 = sigmat1->Integral(tmin,tmax);
Double_t int2 = sigmat2->Integral(tmin,tmax);
Double_t int3 = sigmat3->Integral(tmin,tmax);
cout << "int1(W=0.32)=" << int1 << " int2(W=0.39)=" << int2 << " int3(W=0.47)=" << int3 << endl;
TLegend *leg = new TLegend(0.4,0.7,0.6,0.9);
leg->AddEntry(sigmat,"W=0.32 GeV","l");
leg->AddEntry(sigmat2,"W=0.39 GeV","l");
leg->AddEntry(sigmat3,"W=0.47 GeV","l");
leg->Draw();
sprintf (string,"E_{#gamma}=%.2f GeV\n",Eg);
printf("string=%s",string);
TLatex *t1 = new TLatex(0.7,0.84,string);
t1->SetNDC();
t1->SetTextSize(0.04);
t1->Draw();
sprintf (string,"c=%.3f fm\n",c0);
printf("string=%s",string);
t1->DrawLatex(0.7,0.8,string);
sprintf (string,"a=%.3f fm\n",z0);
printf("string=%s",string);
t1->DrawLatex(0.7,0.76,string);
const Int_t nsigma=700;
Double_t sigmaP[nsigma];
Double_t sigmaW[nsigma];
Double_t Wpipi_min=0.2;
Double_t Wpipi_max=0.9;
Eg = 6;
c0 = 6.62; // Pb half-density radius, fm
z0 = 0.546; // Pb difuseness parameter, fm
w0 = 0;
for (Int_t jj=0; jj<nsigma; jj++) {
Wpipi = Wpipi_min + jj*0.001;
sigmat->SetParameters(Wpipi,Eg,c0,z0,w0,A,mu,kT);
Double_t sigma = sigmat->Integral(tmin,tmax);
// cout << " jj=" << jj << " Wpipi=" << Wpipi << " sigma=" << sigma << endl;
// sigmaP[jj]=sigma/1.e6; // convert nb to mb
sigmaP[jj]=sigma; // units of nb
sigmaW[jj]=Wpipi;
}
TGraph *gr_sigmaP = new TGraph(nsigma,sigmaW,sigmaP);
TCanvas *c4 = new TCanvas("c4", "c4",200,10,1000,700);
// c1->cd(2);
gPad->SetLogy();
gPad->SetGridx();
gPad->SetGridy();
ymin = 1;
ymax = 1e6;
gr_sigmaP->SetTitle("#gamma Pb #rightarrow Pb #pi^{0} #pi^{0}");
// gr_sigmaP->GetXaxis()->SetRangeUser(xmin,xmax);
gr_sigmaP->GetYaxis()->SetRangeUser(ymin,ymax);
gr_sigmaP->GetXaxis()->SetTitleSize(0.05);
gr_sigmaP->GetYaxis()->SetTitleSize(0.05);
gr_sigmaP->GetXaxis()->SetTitle("M_{#pi#pi} (GeV)");
gr_sigmaP->GetYaxis()->SetTitle("#sigma (nb/GeV)");
gr_sigmaP->SetLineColor(4);
gr_sigmaP->Draw();
const Int_t nCounts=14;
Double_t xCounts[nCounts]={0.225,0.275,0.325,0.375,0.425,0.475,0.525,0.575,0.625,0.675,0.725,0.775,0.825,0.875};
Double_t yCounts[nCounts];
Double_t xCounts_err[nCounts];
Double_t yCounts_err[nCounts];
Double_t flux= 1e7; // 10^7 Hz
// Double_t ttarg=7.3 * 0.06; // target thickness in g/cm2 5% Sn target
Double_t ttarg=11.35 * 0.028; // target thickness in g/cm2 5% Pb target
Double_t N0=6.02e23; //
Double_t A_Pb=207.2;
Double_t Rconst = flux*1e-33*ttarg*N0/A_Pb; // use A of lead for consistency with cross section. sigma is in nb
Double_t time=20*3600*24; // Use 20 PAC days
Double_t Int_sigmaP=0;
Double_t Delta = sigmaW[1]-sigmaW[0];
Double_t acceptance=0;
cout << " Delta=" << Delta << endl;
for (Int_t kk=0; kk<nCounts; kk++) {
Double_t Delta_sigma = 0;
for (Int_t jj=0; jj<50; jj++) {
Int_sigmaP += sigmaP[kk*50+jj]*Delta;
Delta_sigma += sigmaP[kk*50+jj]*Delta*Rconst*time;
// cout << " kk=" << kk << " jj=" << jj << " ndx=" << kk*50+jj << " Int_sigmaP=" << Int_sigmaP << " Delta_sigma=" << Delta_sigma << endl;
}
acceptance = -2.8 + 10*xCounts[kk]; // approximation. Linear increase from 0.28 to 0.34, then plateau at 0.6.
if (acceptance < 0) acceptance = 0;
if (acceptance > 0.6) acceptance = 0.6;
yCounts[kk] = Delta_sigma * acceptance;
}
Double_t Rate_sigmaP= Int_sigmaP*Rconst*acceptance;
Double_t Counts_sigmaP=0;
for (Int_t kk=0; kk<nCounts; kk++) {
Counts_sigmaP += yCounts[kk];
}
cout << " Int_sigmaP=" << Int_sigmaP << " Rate_sigmaP=" << Rate_sigmaP << " Rate_sigmaP_count=" << Rate_sigmaP*time << " Counts_sigmaP=" << Counts_sigmaP << endl;
for (Int_t jj=0; jj<nCounts; jj++) {
yCounts_err[jj] = sqrt(yCounts[jj]);
xCounts_err[jj] = 0.001;
}
TGraphErrors *gr_Counts = new TGraphErrors(nCounts,xCounts,yCounts,xCounts_err,yCounts_err);
TCanvas *c5 = new TCanvas("c5", "c5",200,10,1000,700);
ymin=1;
ymax=1e4;
// c1->cd(2);
gPad->SetLogy();
gPad->SetGridx();
gPad->SetGridy();
gr_Counts->SetTitle("#gamma Pb #rightarrow Pb #pi^{0} #pi^{0}");
// gr_Counts->GetXaxis()->SetRangeUser(xmin,xmax);
gr_Counts->GetYaxis()->SetRangeUser(ymin,ymax);
gr_Counts->GetXaxis()->SetTitleSize(0.05);
gr_Counts->GetYaxis()->SetTitleSize(0.05);
gr_Counts->GetXaxis()->SetTitle("M_{#pi#pi} (GeV)");
gr_Counts->GetYaxis()->SetTitle("Counts/50 MeV");
gr_Counts->SetMarkerStyle(20);
gr_Counts->SetMarkerColor(2);
gr_Counts->SetMarkerSize(1.2);
gr_Counts->Draw();
TString text;
text.Form("Flux=%.1g/s",flux);
TLatex *t2 = new TLatex(0.6,3000,text);
t2->SetTextSize(0.05);
t2->Draw();
text.Form("t= %.3f g/cm^{2}",ttarg);
t2->DrawLatex(0.6,1200,text);
text.Form("#sigma_{tot}= %.0f nb/nucleus",Int_sigmaP);
t2->DrawLatex(0.6,600,text);
text.Form("Time=20 PAC days");
t2->DrawLatex(0.6,200,text);
text.Form("Accepted Events=%.0f",Counts_sigmaP);
t2->DrawLatex(0.2,3000,text);
// estimate uncertainties on sigma_gg by scaling fractional uncertainties in sigma_gPb
const Int_t nExpect=12;
Double_t xExpect[nExpect]={0.225,0.275,0.325,0.375,0.425,0.475,0.525,0.575,0.625,0.675,0.725,0.775};
Double_t yExpect[nExpect];
Double_t xExpect_err[nExpect];
Double_t yExpect_err[nExpect];
Double_t offset = 0.1*(xExpect[1]-xExpect[0]);
for (Int_t kk=0; kk<nExpect; kk++) {
Double_t x = xExpect[kk] + offset;
xExpect[kk] = x;
xExpect_err[kk] = xCounts_err[kk];
yExpect[kk] = sigma_ggpi0pi0->Eval(x);
yExpect_err[kk] = yCounts[kk] > 0? yExpect[kk] * yCounts_err[kk]/yCounts[kk] : 0; // scale to uncertainty of measurements
cout << " x=" << x << " Expect=" << yExpect[kk] << " yExpect_err=" << yExpect_err[kk] << " offset=" << offset << endl; }
TGraphErrors *gr_Expect = new TGraphErrors(nExpect,xExpect,yExpect,xExpect_err,yExpect_err);
TCanvas *c6 = new TCanvas("c6", "c6",200,10,1000,700);
ymin=0.;
ymax=20;
gr_Expect->SetTitle("#gamma #gamma #rightarrow #pi^{0} #pi^{0}");
// gr_Expect->GetXaxis()->SetRangeUser(xmin,xmax);
gr_Expect->GetYaxis()->SetRangeUser(ymin,ymax);
gr_Expect->GetXaxis()->SetTitleSize(0.05);
gr_Expect->GetYaxis()->SetTitleSize(0.05);
gr_Expect->GetXaxis()->SetTitle("M_{#pi#pi} (GeV)");
gr_Expect->GetYaxis()->SetTitle("#sigma (nb)");
gr_Expect->SetLineColor(2);
gr_Expect->SetMarkerStyle(20);
gr_Expect->SetMarkerColor(2);
// gr_Expect->FixParameter(1,mu);
gr_Expect->Draw("Ap");
// gr_Marsiske->Fit("gr_Expect");
TGraphErrors *gr_Marsiske2 = (TGraphErrors *)gr_Marsiske->Clone("gr_Marsiske2");
gr_Marsiske2->SetMarkerColor(1);
gr_Marsiske2->SetMarkerStyle(24);
gr_Marsiske2->Draw("psame");
TLegend *leg1 = new TLegend(0.15,0.75,0.55,0.9);
leg1->AddEntry(gr_Marsiske2,"Crystal Ball 1990 / 0.8","p");
leg1->AddEntry(gr_Expect,"Expected Statistical Uncertainties","p");
leg1->Draw();
c1->SaveAs(filename+".pdf(");
// c3->SaveAs(filename+".pdf");
c4->SaveAs(filename+".pdf");
c5->SaveAs(filename+".pdf");
c6->SaveAs(filename+".pdf)");
}
More information about the Halld-npp
mailing list