Double_t ff_func (Double_t *x, Double_t *par){
    
    // return the nuclear form factor accourding to 2 parameter Fermi distribution
    // See Journall 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 R0  = par[0];  				 // half-density radius
    Double_t a0 = par[1];                    // skin or diffuseness parameter
    
    Double_t rho0;
    Double_t sum=0;
    Int_t jmax=4;
    for (j=1;j<jmax;j++) {                    // truncate after 3 terms, first term dominates.
        Double_t sum1 =pow(-1.,j-1)*exp(-j*R0/a0)/(j*j*j);
        sum += sum1;
        // cout << "jmax=" << jmax << " j=" << j << " R0=" << R0 << " a0=" << a0 << " sum1=" << sum1 << " sum=" << sum << endl;
    }
    
    rho0 = (4*pi*R0/3)*( pi*a0*pi*a0  + R0*R0 + 8*pi*a0*a0*a0*sum);
    rho0 = 1/rho0;
    
    Double_t f = 0;
    
    f = (4*pi*pi*rho0*a0*a0*a0)/(q*q*a0*a0*sinh(pi*q*a0)*sinh(pi*q*a0))
    	* (pi*q*a0 * cosh(pi*q*a0) * sin(q*R0) - q*R0 *cos(q*R0) * sinh(pi*q*a0) )
    	+ 8*pi*rho0*a0*a0*a0*sum;
    
    // cout << " q=" << q << " f=" << f << endl;
    return f;
    
}
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 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[2];
    
    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
    
    Double_t FF = ff_func(xin,parin);
    
    // include other masses here
    
    Double_t m1 = 0;      // mass of photon, incident beam
    Double_t m2 = 108.;    // mass of 116Sn, target
    Double_t m3 = Wpipi;   // mass of 2pi system, scattered system
    Double_t m4 = m2;     // recoil target
    
    
    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 << " thepipi=" << thepipi << " f=" << f << endl;
    // cout << " t=" << t << " FF=" << FF << " f=" << f << endl;
    return	f;
}
void sigma_2pi(void)
{
// File: sigma_2pi.C
    // Compute the cross section for Primakoff production of pi+pi- in gamma A -> A pi+ pi- reaction
    // 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_2pi_figs";
    // define function
    
    // Define historgrams
    Int_t nbins=20;
    
    TH1F *h1_perp_chi2 = new TH1F("h1_perp_chi2","PERP Chi2",nbins,0,100);
    
    // define sigma_2pi distribution curve to plot
    
    Int_t npar = 4;
    Double_t Wpipi = 0.32;
    Double_t Eg = 6;
    // Double_t R0 = 6.62;   // Pb half-density radius, fm
    // Double_t a0 = 0.546;   // Pb difuseness parameter, fm
    Double_t R0 = 5.358;   // Sn half-density radius, fm
    Double_t a0 = 0.550;   // Sn difuseness parameter, fm
    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,R0,a0);
    sigmat->SetParNames("Wpipi","Eg","R0","a0");
    
    tmin=0;               // use -t (positive) as input
    tmax=0.005;
    npar = 2;
    // tmax=0.2;
    TF1 *ff = new TF1("ff",ff_func,tmin,tmax,npar);
    ff->SetParameters(R0,a0);
    ff->SetParNames("R0","a0");
    
    Double_t hist_mean;
    Double_t chi2_mean;
    
    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 *sigmat2 = sigmat->DrawCopy("same");
    sigmat2->SetLineColor(2);
    Wpipi=0.39;
    sigmat2->SetParameters(Wpipi,Eg,R0,a0);
    TF1 *sigmat3 = sigmat->DrawCopy("same");
    sigmat3->SetLineColor(1);
    Wpipi=0.47;
    sigmat3->SetParameters(Wpipi,Eg,R0,a0);
    
    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,"c=%.3f fm\n",R0);
    printf("string=%s",string);
    TLatex *t1 = new TLatex(0.7,0.8,string);
    t1->SetNDC();
    t1->SetTextSize(0.04);
    t1->Draw();
    sprintf (string,"a=%.3f fm\n",a0);
    printf("string=%s",string);
    t1->DrawLatex(0.7,0.76,string);
    
    
    
    TCanvas *c3 = new TCanvas("c3", "c3",200,10,1000,700);
    
    // c1->cd(2);
    gPad->SetLogy();
    
    ff->SetParameters(R0,a0);
    ff->SetTitle("");
    // ff->GetXaxis()->SetRangeUser(xmin,xmax);
    // ff->GetYaxis()->SetRangeUser(ymin,ymax);
    ff->GetXaxis()->SetTitleSize(0.05);
    ff->GetYaxis()->SetTitleSize(0.05);
    ff->GetXaxis()->SetTitle("-t (GeV^{2})");
    ff->GetYaxis()->SetTitle("F(q)");
    ff->SetMarkerColor(4);
    ff->Draw();
    
    sprintf (string,"c=%.3f fm\n",R0);
    printf("string=%s",string);
    t1->DrawLatex(0.7,0.8,string);
    sprintf (string,"a=%.3f fm\n",a0);
    printf("string=%s",string);
    t1->DrawLatex(0.7,0.76,string);
    
    c2->SaveAs(filename+".pdf(");
    c3->SaveAs(filename+".pdf)");

}