{

Double_t Pxdp1 = 0.99677;
Double_t Pmuxi = 0.9960;
Double_t mL = 80.425;  // Mass of the left-handed W
Double_t mRlow = 380;
Double_t nmR = 8000;
Double_t nzeta = 10000;
Double_t Maxzeta = 0.40;
Double_t ZeroDef = 0.0001;

// Coordinates of curve
Double_t mR[nmR];
Double_t zetaNeg[nmR];
Double_t zetaPos[nmR];

Double_t m1;
Double_t m2;
Double_t epsilon;
Double_t zeta;
Double_t xi;
Double_t Pmu;
Int_t imin;
imin = 0;

for(Int_t imR = 0; imR < nmR; imR++){
  mR[imR] = (imR)*(1200-mRlow)/(nmR-1)+mRlow;
  //  epsilon = pow(mL/mR[imR],2);
  for(Int_t izeta = -nzeta+1; izeta <= nzeta-1; izeta++){
    zeta = izeta*(Maxzeta/(nzeta-1));
    m1 = mL*pow(cos(zeta),2) + mR[imR]*pow(sin(zeta),2);
    m2 = mL*pow(sin(zeta),2) + mR[imR]*pow(cos(zeta),2);
    //    printf("mL = %f  m1 = %f  mR = %f  m2 = %f\n",mL,m1,mR[imR],m2);
    epsilon = pow(m1/m2,2);
    xi = 1 - 2*pow(epsilon,2) - 2*pow(zeta,2);
    Pmu = 1 - 2*pow((epsilon + zeta),2);
    if(sqrt(pow(Pmu*xi - Pmuxi,2)) < ZeroDef){
      // on 90% confidence limit curve
      zetaNeg[imR] = zeta;
      break;
    }
    else if(izeta == 0){
      // mR too low
      imin = imR;
    }
  }
  for(Int_t jzeta = nzeta-1; jzeta >= -nzeta+1; jzeta--){
    zeta = jzeta*(Maxzeta/(nzeta-1));
    m1 = mL*pow(cos(zeta),2) + mR[imR]*pow(sin(zeta),2);
    m2 = mL*pow(sin(zeta),2) + mR[imR]*pow(cos(zeta),2);
    epsilon = pow(m1/m2,2);
    xi = 1 - 2*pow(epsilon,2) - 2*pow(zeta,2);
    Pmu = 1 - 2*pow((epsilon + zeta),2);
    if(sqrt(pow(Pmu*xi - Pmuxi,2)) < ZeroDef){
      // on 90% confidence limit curve
      zetaPos[imR] = zeta;
      break;
    }
  }
}

printf("Lower limit on mass of WR is %f\n",mR[imin+1]);

Int_t nmR2 = nmR - imin - 1;
Double_t mR2[nmR2];
Double_t zN[nmR2];
Double_t zP[nmR2];

{
  Int_t kmR;
  for(Int_t jmR = 0; jmR < nmR2; jmR++){
    kmR = jmR + imin;
    mR2[jmR] = mR[kmR];
    zN[jmR] = zetaNeg[kmR];
    zP[jmR] = zetaPos[kmR];
  }
}

TGraph *gTWNeg;
gTWNeg = new TGraph(nmR2, mR2, zN);

TGraph *gTWPos;
gTWPos = new TGraph(nmR2, mR2, zP);

TH2D *hex1;
hex1 = new TH2D("hex1","",6,0.0,1200.,10,-0.1,0.1);
hex1->SetStats(kFALSE);
hex1->GetXaxis()->SetTitle("Mass of right-handed W (GeV/c)");
hex1->GetXaxis()->CenterTitle();
hex1->GetXaxis()->SetTitleSize(0.05);
hex1->GetXaxis()->SetLabelSize(0.05);
hex1->GetYaxis()->SetTitle("Left-right mixing angle, #zeta");
hex1->GetYaxis()->CenterTitle();
hex1->GetYaxis()->SetTitleSize(0.05);
hex1->GetYaxis()->SetLabelSize(0.05);
hex1->GetYaxis()->SetTitleOffset(1.2);

TCanvas *cTWIST;
cTWIST = new TCanvas("cTWIST","Left-right exclusion plot,  TWIST limits",600,400);
cTWIST->SetFillColor(0);
cTWIST->SetFrameBorderMode(0);
cTWIST->SetBorderMode(0);
cTWIST->SetBorderSize(0);
cTWIST->SetRightMargin(0.04);
cTWIST->SetLeftMargin(0.12);
cTWIST->SetTopMargin(0.03);
cTWIST->SetBottomMargin(0.12);
hex1->Draw();
gTWNeg->Draw("P");
gTWPos->Draw("P");

}