//
void magbeam_FlatTop_xy(){
  
  gROOT->Reset();
  
#define MAX_nentries 100000;
  
  TString pstit;
  
  // define style here 
  // general parameters
  gStyle->SetOptDate(0);     gStyle->SetOptTitle(0);
  gStyle->SetStatColor(10);  gStyle->SetStatH(0.2);
  gStyle->SetStatW(0.3);     
  gStyle->SetOptStat(111); 
  
  // canvas parameters
  gStyle->SetFrameBorderMode(0);
  gStyle->SetFrameBorderSize(0);
  gStyle->SetFrameFillColor(0);
  gStyle->SetCanvasColor(0);
  
  //
  
  // pads parameters
  //  gStyle->SetPadColor(39); 
  gStyle->SetPadColor(0); 
  gStyle->SetPadBorderMode(0);
  gStyle->SetPadBorderSize(0);
  gStyle->SetPadBottomMargin(0.18);
  gStyle->SetPadRightMargin(0.10);
  gStyle->SetPadLeftMargin(0.14);
  gStyle->SetLabelSize(0.05,"x");
  gStyle->SetLabelSize(0.05,"y");
  gStyle->SetTitleXSize(0.08);
  gStyle->SetPaperSize(10,12);
  
  gStyle->SetOptLogx(0);
  gStyle->SetOptLogy(0);
  
  gStyle->SetTitleYOffset(0.8);
  gStyle->SetTitleYSize(0.10);
  gROOT->ForceStyle();  
  
  ////  
  
  Double_t c = 299792458.0; // Speed of Light (m/s)
  Double_t e = 1.602176e-19; // Electron Charge (C)
  Double_t mec2 = 0.510998e+6; // Electron Mass (eV)
  Double_t B_cath = 0.2; // Cathode Solenoid Field (T)
  
  Double_t eOverme = 1.758820088e+11; // Electron Cyclotron Frequency/Field ( rad/(s T) )
  
  //// 
  
  TH1F* h1 = new TH1F("h1","h1", 100,0.0,2.0);
  TH1F* h2 = new TH1F("h2","h2", 100,0.0,5.0);
  TH2F* h3 = new TH2F("h3","h3", 100,-2.0,2.0, 100,-2.0,2.0);
  TH1F* h4 = new TH1F("h4","h4", 100,0.0,100.0);
  TH1F* h5 = new TH1F("h5","h5", 100,-80.0,80.0);
  TH1F* h6 = new TH1F("h6","h6", 100,-80.0,80.0);
  TH1F* h7 = new TH1F("h7","h7", 100,0.0,300.0);
  TH1F* h8 = new TH1F("h8","h8", 100,0.0,100.0);
  TH1F* h9 = new TH1F("h9","h9", 100,0.0,100.0);

  TH1F* eh1 = new TH1F("eh1","eh1", 100000,0.0,1.0e-3);
  TH1F* eh2 = new TH1F("eh2","eh2", 100000,0.0,1.0e+6);
  TH1F* eh3 = new TH1F("eh3","eh3", 100000,-10.0,10.0);

  ////
  
  ofstream outfile;    
  char ofile[200];
  sprintf(ofile,"MagBeam_FlatTop_xy.txt");
  cout << "Data file is " << ofile << endl;
  
  outfile.open(ofile);
  outfile << " No. \t" << " x (mm) \t" << " y (mm) \t" << " px (keV/c) \t" << " py (keV/c) \t" << endl;
 
  ////

  Double_t a0 = 0.00156; // Flat-top beam radius (m)
  Double_t r; // Electron Radius
  
  for (Int_t i=0; i<MAX_nentries; i++){
    Double_t u = gRandom->Uniform(0.0,a0) + gRandom->Uniform(0.0,a0); 
    if (u>a0){r=2.0*a0-u;} else {r=u;} // Generate a random point within a circle (uniformly)
    //    http://stackoverflow.com/questions/5837572/generate-a-random-point-within-a-circle-uniformly

    Double_t phi = gRandom->Uniform( 0.0, 2.0*TMath::Pi() ); // Generate electron angle
    
    Double_t x = r*TMath::Cos(phi);
    Double_t y = r*TMath::Sin(phi);
    
    h1->Fill(r*1.0e+3); // mm, <r> = 2/3 a0
    h2->Fill(r*1.0e+3*r*1.0e+3); // mm, <r*r> = a0*a0/2
    h3->Fill(x*1.0e+3,y*1.0e+3); // mm
    
    //
    
    // Phi Momentum (keV/c)
    Double_t phi_mom = 1.0/2.0 * eOverme * mec2/c**2 * B_cath * r * c * 1.0e-3;
    //    cout << " Phi Momentum: phi_mom = " << phi_mom << " (eV/c)" << endl;

    h4->Fill(phi_mom);
    
    ////
    
    Double_t px = -1.0 * phi_mom*TMath::Sin(phi);
    Double_t py = phi_mom*TMath::Cos(phi);
      
    h5->Fill(px);
    h6->Fill(py);

    ////

    outfile << i << "\t"  << x*1.0e+3 << "   \t" << y*1.0e+3 << "   \t" << px << "   \t" << py << endl;

    // Canonical Angular Momentum at Cathode (neV s)
    Double_t L_mom = 1.0e+9 * 1.0/2.0 * eOverme * mec2/c**2 * B_cath * r**2;
    //   cout << " Canonical Angular Momentum at Cathode: L_mom = " << L_mom << " (neV s)" << endl;
    
    h7->Fill(L_mom);
    
    // Cathode Magnetic Emittance (um)
    Double_t mag_emit_cathode = 1.0e+6 * (eOverme * B_cath)/(4.0*c) * r**2;
    //    cout << " Cathode Magnetic Emittance: mag_emit_cathode = " << mag_emit_cathode << " um" << endl;
    
    h8->Fill(mag_emit_cathode);
 
    // Beam Emittance x (um)
   
    eh1->Fill(x*x);
    eh2->Fill(px*px);
    eh3->Fill(x*px);
    
  }

  ///

  Double_t x2_ave  = eh1->GetMean();
  Double_t px2_ave = eh2->GetMean() * 1.0e+6; // keV/c to eV/c
  Double_t xpx_ave = eh3->GetMean() * 1.0e+3; // keV/c to eV/c
 
  cout << " <x2> : x2_ave = " << x2_ave << " m2" << endl;
  cout << " <px2> : px2_ave = " << px2_ave << " eV/c2" << endl;
  cout << " <xpx> : xpx_ave = " << xpx_ave << " m.eV/c" << endl;

  Double_t emit_x = 1.0e+6 * 1.0/mec2 * TMath::Sqrt(x2_ave * px2_ave - xpx_ave * xpx_ave);
  cout << " Emittance x : emit_x = " << emit_x << " um" << endl;
  h9->Fill(emit_x);

  cout << " Cathode Magnetic Emittance: mag_emit_cathode = " <<  h8->GetMean() << " um" << endl;
  
  ///////
  
  a1 = new TCanvas("a1","Magnetized Flat-top Beam",40,40,1000,1000);
  a1->Divide(3,3);
  
  a1->cd(1); 
  
  gPad->SetLogy(0);
  
  h1->SetMarkerColor(kRed);
  h1->SetLineColor(kRed);
  h1->SetLineWidth(1);
  h1->SetMarkerSize(1.0);
  h1->SetMarkerStyle(20);
  
  h1->GetXaxis()->CenterTitle();
  h1->GetXaxis()->SetTitleOffset(1.2);
  h1->GetXaxis()->SetTitleSize(0.06);
  TString tit = "r (mm)";
  h1->GetXaxis()->SetTitle(tit.Data());
  
  h1->GetYaxis()->CenterTitle();
  h1->GetYaxis()->SetTitleOffset(1.0);
  h1->GetYaxis()->SetTitleSize(0.06);
  TString tit = "";
  h1->GetYaxis()->SetTitle(tit.Data());
  
  h1->Draw("histo");
  
  // //  
  
  a1->cd(2); 
  
  gPad->SetLogy(0);
  
  h2->SetMarkerColor(kGreen);
  h2->SetLineColor(kGreen);
  h2->SetLineWidth(1);
  h2->SetMarkerSize(1.0);
  h2->SetMarkerStyle(20);
  
  h2->GetXaxis()->CenterTitle();
  h2->GetXaxis()->SetTitleOffset(1.2);
  h2->GetXaxis()->SetTitleSize(0.06);
  TString tit = "r^{2} (mm^{2})";
  h2->GetXaxis()->SetTitle(tit.Data());
  
  h2->GetYaxis()->CenterTitle();
  h2->GetYaxis()->SetTitleOffset(1.0);
  h2->GetYaxis()->SetTitleSize(0.06);
  TString tit = "";
  h2->GetYaxis()->SetTitle(tit.Data());
  
  h2->Draw("histo");
  
  // //  
  
  a1->cd(3); 
  gPad->SetGrid();
  
  h3->SetLineWidth(1);
  h3->SetMarkerSize(1.0);
  h3->SetMarkerStyle(20);
  
  h3->GetXaxis()->CenterTitle();
  h3->GetXaxis()->SetTitleOffset(1.2);
  h3->GetXaxis()->SetTitleSize(0.06);
  TString tit = "x (mm)";
  h3->GetXaxis()->SetTitle(tit.Data());
  
  h3->GetYaxis()->CenterTitle();
  h3->GetYaxis()->SetTitleOffset(1.0);
  h3->GetYaxis()->SetTitleSize(0.06);
  TString tit = "y (mm)";
  h3->GetYaxis()->SetTitle(tit.Data());
  
  h3->Draw("COL");
  
  // //  
  
  a1->cd(4); 
  
  gPad->SetLogy(0);
  
  h4->SetMarkerColor(kCyan);
  h4->SetLineColor(kCyan);
  h4->SetLineWidth(1);
  h4->SetMarkerSize(1.0);
  h4->SetMarkerStyle(20);
  
  h4->GetXaxis()->CenterTitle();
  h4->GetXaxis()->SetTitleOffset(1.2);
  h4->GetXaxis()->SetTitleSize(0.06);
  TString tit = "P_{#phi} (keV/c)";
  h4->GetXaxis()->SetTitle(tit.Data());
  
  h4->GetYaxis()->CenterTitle();
  h4->GetYaxis()->SetTitleOffset(1.0);
  h4->GetYaxis()->SetTitleSize(0.06);
  TString tit = "";
  h4->GetYaxis()->SetTitle(tit.Data());
  
  h4->Draw("histo");

  // //  
  
  a1->cd(5); 
  
  gPad->SetLogy(0);
  
  h5->SetMarkerColor(kCyan);
  h5->SetLineColor(kCyan);
  h5->SetLineWidth(1);
  h5->SetMarkerSize(1.0);
  h5->SetMarkerStyle(20);
  
  h5->GetXaxis()->CenterTitle();
  h5->GetXaxis()->SetTitleOffset(1.2);
  h5->GetXaxis()->SetTitleSize(0.06);
  TString tit = "p_{x} (keV/c)";
  h5->GetXaxis()->SetTitle(tit.Data());
  
  h5->GetYaxis()->CenterTitle();
  h5->GetYaxis()->SetTitleOffset(1.0);
  h5->GetYaxis()->SetTitleSize(0.06);
  TString tit = "";
  h5->GetYaxis()->SetTitle(tit.Data());
  
  h5->Draw("histo");

  // //  
  
  a1->cd(6); 
  
  gPad->SetLogy(0);
  
  h6->SetMarkerColor(kCyan);
  h6->SetLineColor(kCyan);
  h6->SetLineWidth(1);
  h6->SetMarkerSize(1.0);
  h6->SetMarkerStyle(20);
  
  h6->GetXaxis()->CenterTitle();
  h6->GetXaxis()->SetTitleOffset(1.2);
  h6->GetXaxis()->SetTitleSize(0.06);
  TString tit = "p_{y} (keV/c)";
  h6->GetXaxis()->SetTitle(tit.Data());
  
  h6->GetYaxis()->CenterTitle();
  h6->GetYaxis()->SetTitleOffset(1.0);
  h6->GetYaxis()->SetTitleSize(0.06);
  TString tit = "";
  h6->GetYaxis()->SetTitle(tit.Data());
  
  h6->Draw("histo");
  
  // //  
  
  a1->cd(7); 
  
  gPad->SetLogy(0);

  h7->SetMarkerColor(kMagenta);
  h7->SetLineColor(kMagenta);
  h7->SetLineWidth(1);
  h7->SetMarkerSize(1.0);
  h7->SetMarkerStyle(20);
  
  h7->GetXaxis()->CenterTitle();
  h7->GetXaxis()->SetTitleOffset(1.2);
  h7->GetXaxis()->SetTitleSize(0.06);
  TString tit = "Canonical Angular Momentum (neV s)";
  h7->GetXaxis()->SetTitle(tit.Data());
  
  h7->GetYaxis()->CenterTitle();
  h7->GetYaxis()->SetTitleOffset(1.0);
  h7->GetYaxis()->SetTitleSize(0.06);
  TString tit = "";
  h7->GetYaxis()->SetTitle(tit.Data());
  
  h7->Draw("histo");
  
  // //  
  
  a1->cd(8); 
  
  gPad->SetLogy(0);
  
  h8->SetMarkerColor(kBlue);
  h8->SetLineColor(kBlue);
  h8->SetLineWidth(1);
  h8->SetMarkerSize(1.0);
  h8->SetMarkerStyle(20);
  
  h8->GetXaxis()->CenterTitle();
  h8->GetXaxis()->SetTitleOffset(1.2);
  h8->GetXaxis()->SetTitleSize(0.06);
  TString tit = "Cathode Magnetic Emittance (#mum)";
  h8->GetXaxis()->SetTitle(tit.Data());
  
  h8->GetYaxis()->CenterTitle();
  h8->GetYaxis()->SetTitleOffset(1.0);
  h8->GetYaxis()->SetTitleSize(0.06);
  TString tit = "";
  h8->GetYaxis()->SetTitle(tit.Data());
  
  h8->Draw("histo");

  // //  
  
  a1->cd(9); 
  
  gPad->SetLogy(0);
  
  h9->SetMarkerColor(kBlue);
  h9->SetLineColor(kBlue);
  h9->SetLineWidth(1);
  h9->SetMarkerSize(1.0);
  h9->SetMarkerStyle(20);
  
  h9->GetXaxis()->CenterTitle();
  h9->GetXaxis()->SetTitleOffset(1.2);
  h9->GetXaxis()->SetTitleSize(0.06);
  TString tit = "Beam Emittance #epsilon_{x} (#mum)";
  h9->GetXaxis()->SetTitle(tit.Data());
  
  h9->GetYaxis()->CenterTitle();
  h9->GetYaxis()->SetTitleOffset(1.0);
  h9->GetYaxis()->SetTitleSize(0.06);
  TString tit = "";
  h9->GetYaxis()->SetTitle(tit.Data());
  
  h9->Draw("histo");

  // //   
  // //
  
  a1->cd();
  pstit = "magbeam_FlatTop_xy.gif";
  a1->Print(pstit.Data());
  ///////
}