Ticket #7: CalcMultiPhaseSpace.C

// Calculate the multibody phase space integral (4pi or "rho,rho") 
// for the K-matrix

double mPi = 0.13957018;
double mPiSq = mPi*mPi;
double mPiSq4 = 4.0*mPiSq;
double mPiSq16 = 16.0*mPiSq;
double mPi2 = 2.0*mPi;
double mPi4 = 4.0*mPi;
double pi = acos(-1.0);
double piSq = pi*pi;

TCanvas* theCanvas = new TCanvas("theCanvas", "", 800, 600);

void CalcMultiPhaseSpace() {

    // Here, s is the invariant mass-squared of the whole system
    // while s1 and s2 are the invariant mass-squared of the two di-pion
    // parts
    // Width Gamma(s) = rho_1^3, where rho_1  = sqrt(1.0 - 4*m_pi^2/s)
    // For a given s, s1 can vary between (2m_pi)^2 and s, while s2 varies
    // between (2m_pi)^2 and sqrt(s1) - sqrt(s), since m_0 >= m_1 + m_2
    // and s = m_0^2, s1 = m_1^2, s2 = m_2^2

    // Find the phase space multi-integral for s between ~0 and 1
    double ds = 1e-3;
    double sMin = 0.0; //mPiSq16;
    double sMax = 1.0 + 2.0*ds;
    int N0 = int((sMax - sMin)/ds);
    cout<<"N0 = "<<N0<<endl;

    // The nominal pole mass of the rho meson
    double M = 0.775;
    double MSq = M*M;

    // Do the double integral
    int i, j, k;

    double ds1 = 1e-3;
    double ds2 = 1e-3;

    double s1Min = mPiSq4; // (2*m_pi)^2
    double s2Min = mPiSq4;

    // Create a graph for the values
    theCanvas->Clear();
    gROOT->SetStyle("Plain");
    gStyle->SetOptStat(0);
    theCanvas->UseCurrentStyle();

    TGraph* graph = new TGraph(N0);

    for (i = 0; i < N0; i++) {

        if (i%100 == 0) {cout<<"i = "<<N0 - i <<endl;}

        double s = sMin + ds*i;
        double integral(0.0);

        if (s < 0.0) {continue;}
 
        // sqrt(s1) will vary between 2*mpi and sqrt(s) - 2*mpi.
        // The 2*mpi factor is due to the minimum mass needed for s2.
        double roots_2mpi = sqrt(s) - mPi2;
        double s1Max = roots_2mpi*roots_2mpi;
        //cout<<"s1Min = "<<s1Min<<", s1Max = "<<s1Max<<endl;

        int N1 = int((s1Max - s1Min)/ds1);
        //cout<<"N1 = "<<N1<<endl;

        for (j = 0; j < N1; j++) {

            double s1 = s1Min + ds1*j;
            double gamma1 = getGamma(s1);

            double Ms1 = MSq - s1;
            double denom1 = Ms1*Ms1 + MSq*gamma1*gamma1;

            // Find the maximum s2 value given s and s1
            double s2Max(s2Min);

            if (fabs(s) > 1e-10 && fabs(s1) > 1e-10) {
                s2Max = sqrt(s) - sqrt(s1);
                s2Max *= s2Max;
            }

            int N2 = int((s2Max - s2Min)/ds2);

            for (k = 0; k < N2; k++) {
                
                double s2 = s2Min + ds2*k;
                double gamma2 = getGamma(s2);

                double Ms2 = MSq - s2;
                double denom2 = Ms2*Ms2 + MSq*gamma2*gamma2;

                double denom = s*denom1*denom2;
                
                if (fabs(denom) > 1e-10) {

                    double sTerm = s + s1 - s2;
                    double sqrtTerm = sTerm*sTerm - 4.0*s*s1;
                    
                    double num(0.0);
                    
                    if (sqrtTerm > 0.0) {
                        num = gamma1*gamma2*sqrt(sqrtTerm);
                    }

                    integral += num*ds1*ds2*MSq/(denom*piSq);


                } // denom is not zero

            } // k

        } // j
        
        cout<<"For i = "<<i<<", s = "<<s<<", integral = "<<integral<<endl;
        graph->SetPoint(i, s, integral);

    } // i

    // Calculate the phase space factor at s = 1
    double rho1 = sqrt(1.0 - mPiSq16);
    // Find the same factor from the graph
    int pointNo = int((1.0 - sMin)/ds);
    cout<<"pointNo = "<<pointNo<<endl;

    double sGraph(0.0), rhoGraph(0.0);
    graph->GetPoint(pointNo, sGraph, rhoGraph);

    // Find the scaling factor rho0
    double rho0 = rho1/rhoGraph;

    cout<<"For sGraph = "<<sGraph<<", rhoGraph = "
        <<rhoGraph<<"; rho0 = "<<rho0<<endl;

    // Rescale the graph points
    for (i = 0; i < N0; i++) {

        double x(0.0);
        double y(0.0);
        graph->GetPoint(i, x, y);
        y *= rho0;
        graph->SetPoint(i, x, y);

    }

    // Plot the curves

    theCanvas->cd();

    TH2F* nullH = new TH2F("nullH", "", 2, 0.0, 2.0, 2, 0.0, 1.1);
    nullH->SetXTitle("s (GeV^{2}/c^{4})");
    nullH->Draw();

    graph->SetMarkerStyle(24);
    graph->Draw("p");

    // Fit the function to this distribution
    TF1* fun = new TF1("fun", polFunction, 0.0, 2.0, 7);
    fun->SetLineColor(kRed);
    fun->SetParameters(1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0);

    graph->Fit("fun");

    // Overlay "old" function
    TF1* oldFun = new TF1("oldFun", rhoFun, 0.0, 1.0, 1);
    oldFun->SetParameter(0, 1.0);
    oldFun->SetLineColor(kBlue);
    oldFun->Draw("same");

    // Add the continuation from s = 1
    TF1* endFun = new TF1("endFun", endTerm, 1.0, 2.0, 1);
    endFun->SetParameter(0, 1.0);
    endFun->SetLineColor(kGreen);
    endFun->Draw("same");
 
    // Add the continuation from s = 1
    TF1* endFun2 = new TF1("endFun2", endTerm2, 1.0, 2.0, 1);
    endFun2->SetParameter(0, 1.0);
    endFun2->SetLineColor(kOrange+7);
    endFun2->Draw("same");
    
    // Text labels
    TLatex* latex = new TLatex();
    latex->SetTextSize(0.03);

    latex->SetTextColor(kBlue);
    latex->DrawText(1.2, 0.6, "Old Laura++ implementation");
    latex->SetTextColor(kBlack);
    latex->DrawText(1.2, 0.5, "New integral calculation");
    latex->SetTextColor(kRed);
    latex->DrawText(1.2, 0.4, "New parameterisation");

    latex->SetTextColor(kGreen);
    latex->DrawText(1.2, 0.3, "Sqrt continuation for s > 1");

    latex->SetTextColor(kOrange+7);
    latex->DrawText(1.2, 0.2, "Non-sqrt continuation for s > 1");

    theCanvas->Print("Laura_rho4pi.png");

}

double expFunction(double* x, double* par) {

    double s = x[0];
    double c0 = par[0];
    double c1 = par[1];
    double c2 = par[2];
    double c3 = par[3];
    double c4 = par[4];
    
    double expTerm = (c0*s + c1)*s + c2;

    double rho(0.0);

    if (fabs(expTerm) < 30.0) {

        rho = c3*exp(expTerm) + c4;

    }

    return rho;

}

double polFunction(double* x, double* par) {

    double s = x[0];
    double c0 = par[0];
    double c1 = par[1];
    double c2 = par[2];
    double c3 = par[3];
    double c4 = par[4];
    double c5 = par[5];
    double c6 = par[6];

    double rho(0.0);

    if (s < 1.0) {

        //rho = (((((c0*s + c1)*s + c2)*s + c3)*s + c4)*s + c5)*s + c6;
        rho = ((c0*s + c1)*s + c2)*s + c3;
        rho = ((rho*s + c4)*s + c5)*s + c6;

    } else {

        double ratio(0.0);
        if (fabs(s) > 1e-10) {
            ratio = mPiSq16/s;
        }
        
        if (ratio > 0.0 && ratio < 1.0) {
            rho = sqrt(1.0 - ratio);
        }

    }

    return rho;

}

double endTerm(double* x, double* par) {

    double s = x[0];
    double c0 = par[0];

    double ratio(0.0);
    if (fabs(s) > 1e-10) {
        ratio = mPiSq16/s;
    }

    double result(0.0);
    if (ratio > 0.0 && ratio < 1.0) {
        result = sqrt(1.0 - ratio);
    }

    return result;


}

double endTerm2(double* x, double* par) {

    double s = x[0];
    double c0 = par[0];

    double ratio(0.0);
    if (fabs(s) > 1e-10) {
        ratio = mPiSq16/s;
    }

    double result(0.0);
    if (ratio > 0.0 && ratio < 1.0) {
        result = 1.0 - ratio;
    }

    return result;


}

double getGamma(double s) {

    // The gamma0 factor is the "width" of the 4pi state.
    // The PDG tables give 300 MeV as a typical 4pi width, which 
    // is ~75% of the total width, with large uncertainties
    double gamma0 = 0.3; // 0.3 GeV
    double rho1 = getRho1(s);
    double gamma = gamma0*rho1*rho1*rho1;
    return gamma;

}

double getRho1(double s) {

    double rho(0.0);
    
    double term(0.0);
    if (fabs(s) > 1e-10) {
        term = 1.0 - (mPiSq4/s);
    }

    if (fabs(term) > 1e-10) {
        rho = sqrt(term);
    }

    return rho;

}


double old_rhoFun(double* x, double* par) {

  double rho(0.0);
  double s = x[0];
  double A = par[0];
  
  if (s <= 1.0) {
    double term1 = (((11.3153*s - 21.8845)*s + 16.8358)*s - 6.39017)*s + 1.2274;
    term1 += 0.00370909/(s*s);
    term1 -= 0.111203/s;
    double term2 = 0.82965;
    rho = term1*term2;
  } else {
    rho = sqrt(1.0 - (0.311677/s));
  }

  rho *= A;
  return rho;

}

Double_t rhoFun(Double_t*x, Double_t* par) {

    double mSqAB = x[0];
    double A = par[1];
    double mPi = 0.139570;
 
    if (mSqAB > 1.0) {
        return rho( 4.0 * mPi, mSqAB );
    }
 
    double m4AB = std::pow( mSqAB, 2 );
    double m6AB = std::pow( mSqAB, 3 );
    double m8AB = std::pow( mSqAB, 4 );
    
    double term = 0.0;
    term += 0.00370909 / m4AB;
    term -= 0.111203   / mSqAB;
    term += 1.2274;
    term -= 6.39017    * mSqAB;
    term += 16.8358    * m4AB;
    term -= 21.8845    * m6AB;
    term += 11.3153    * m8AB;
    
    return rho( 4.0 * mPi, 1.0 ) * term;

}

double rho(double mCh, double mSqAB) {

    double rhoSq = 1.0 - std::pow(mCh, 2)/mSqAB;
 
    return std::sqrt(rhoSq);

}