{

Double_t lmar = 0.25;
Double_t rmar = 0.01;
Double_t bmar = 0.25;
Double_t tmar = 0.05;
Double_t lbl = 0.1;
Double_t hpanel = 1.0;
Double_t wpanel = 1.0;
Double_t height;
height = 2*hpanel+bmar+tmar;
Double_t width;
width = 2*wpanel+lmar+rmar;
Double_t hxScale = hpanel/(hpanel+tmar);
Double_t hdxScale = hpanel/(hpanel+bmar);
Double_t Tlbl;
Tlbl = 1.5*lbl;

// This root macro will produce beam characterization plots from geant beam tune
// files found in CAL_DB

Int_t istat;
char CommentLine[100];

// // To plot pre Oct 23 surface beam
// FILE *fBeam;
// fBeam = fopen("preOct23_950_mu_trimmed.dat","r");
// Int_t CountMax = 110;

// // To plot post Oct 23 surface beam
// FILE *fBeam;
// fBeam = fopen("postOct23_950_mu_trimmed.dat","r");
// Int_t CountMax = 20;

// To plot cloud beam
FILE *fBeam;
fBeam = fopen("postOct23_cld_mu_trimmed.dat","r");
Int_t CountMax = 10;

istat = fscanf(fBeam,"%s",&CommentLine);

Int_t nX;
istat = fscanf(fBeam,"%i",&nX);

printf("nX = %i \n",nX);

Float_t x[nX];
Float_t xCount[nX];

for(Int_t iX = 0;iX<nX;iX++){
  istat = fscanf(fBeam,"%f %f",&x[iX],&xCount[iX]);
  //    printf("x = %f and xCount = %f \n", x[iX], xCount[iX]);
}

istat = fscanf(fBeam,"%s",&CommentLine);

Int_t ndX;
istat = fscanf(fBeam,"%i",&ndX);

printf("ndX = %i \n",ndX);

Float_t x2[ndX];
Float_t dx[ndX];
Float_t x2Err[ndX];
Float_t dxErr[ndX];

for(Int_t iX = 0;iX<ndX;iX++){
  istat = fscanf(fBeam,"%f %f %f",&x2[iX],&dx[iX],&dxErr[iX]);
  x2Err[iX] = 2.81;
}

istat = fscanf(fBeam,"%s",&CommentLine);

Int_t nY;
istat = fscanf(fBeam,"%i",&nY);

printf("nY = %i \n",nY);

Float_t y[nY];
Float_t yCount[nY];

for(Int_t iY = 0;iY<nY;iY++){
  istat = fscanf(fBeam,"%f %f",&y[iY],&yCount[iY]);
}

istat = fscanf(fBeam,"%s",&CommentLine);

Int_t ndY;
istat = fscanf(fBeam,"%i",&ndY);

printf("ndY = %i \n",ndY);

Float_t y2[ndY];
Float_t dy[ndY];
Float_t y2Err[ndY];
Float_t dyErr[ndY];

for(Int_t iY = 0;iY<ndY;iY++){
  istat = fscanf(fBeam,"%f %f %f",&y2[iY],&dy[iY],&dyErr[iY]);
  y2Err[iY] = 2.81;
}

gStyle->SetEndErrorSize(0.0);

TGraph *gx;
gx = new TGraph(nX,x,xCount);
gx->SetMarkerStyle(22);
TGraphErrors *gdx;
gdx = new TGraphErrors(ndX,x2,dx,x2Err,dxErr);
TGraph *gy;
gy = new TGraph(nY,y,yCount);
gy->SetMarkerStyle(22);
TGraphErrors *gdy;
gdy = new TGraphErrors(ndY,y2,dy,y2Err,dyErr);

TH2F *hx;
hx = new TH2F("hx","",9,-45.,45.,CountMax,0.001,100.*CountMax);
hx->SetStats(kFALSE);
hx->GetYaxis()->SetTitle("Count");
hx->GetYaxis()->CenterTitle();
hx->GetYaxis()->SetTitleSize(Tlbl*hxScale);
hx->GetYaxis()->SetLabelSize(lbl*hxScale);
hx->GetXaxis()->SetLabelSize(0.0);
hx->GetYaxis()->SetNdivisions(504,kTRUE);
hx->GetXaxis()->SetNdivisions(208,kTRUE);
hx->GetYaxis()->SetTitleOffset(0.6/hxScale);

TH2F *hy;
hy = new TH2F("hy","",9,-45.,45.,CountMax,0.001,100.*CountMax);
hy->SetStats(kFALSE);
hy->GetYaxis()->SetLabelSize(0.0);
hy->GetXaxis()->SetLabelSize(0.0);
hy->GetYaxis()->SetNdivisions(504,kTRUE);
hy->GetXaxis()->SetNdivisions(208,kTRUE);
TH2F *hdx;
hdx = new TH2F("hdx","",9,-45.,45.,9,-45.,45.);
hdx->SetStats(kFALSE);
hdx->GetYaxis()->SetTitle("Angle (mrad)");
hdx->GetYaxis()->CenterTitle();
hdx->GetYaxis()->SetTitleSize(Tlbl*hdxScale);
hdx->GetYaxis()->SetLabelSize(lbl*hdxScale);
hdx->GetXaxis()->SetTitle("x position (mm)");
hdx->GetXaxis()->CenterTitle();
hdx->GetXaxis()->SetTitleSize(Tlbl*hdxScale);
hdx->GetXaxis()->SetLabelSize(lbl*hdxScale);
hdx->GetYaxis()->SetNdivisions(208,kTRUE);
hdx->GetXaxis()->SetNdivisions(208,kTRUE);
hdx->GetXaxis()->SetTitleOffset(0.55/hdxScale);
hdx->GetYaxis()->SetTitleOffset(0.6/hdxScale);

TH2F *hdy;
hdy = new TH2F("hdy","",9,-45.,45.,9,-45.,45.);
hdy->SetStats(kFALSE);
hdy->GetYaxis()->SetLabelSize(0.0);
hdy->GetXaxis()->SetTitle("y position (mm)");
hdy->GetXaxis()->CenterTitle();
hdy->GetXaxis()->SetTitleSize(Tlbl*hdxScale);
hdy->GetXaxis()->SetLabelSize(lbl*hdxScale);
hdy->GetYaxis()->SetNdivisions(208,kTRUE);
hdy->GetXaxis()->SetNdivisions(208,kTRUE);
hdy->GetXaxis()->SetTitleOffset(0.55/hdxScale);

TCanvas *cbeam;
cbeam = new TCanvas("cbeam","",600,400);
cbeam->SetFillColor(0);
cbeam->SetBorderSize(0.0);
cbeam->SetBorderMode(0);
cbeam->Divide(2,2);

cbeam->cd(1);
TPad *p1;
p1 = gPad;
p1->SetPad("p1","",0.0,(hpanel+bmar)/height,(wpanel+lmar)/width,1.0);
p1->SetFillColor(0);
p1->SetBorderMode(0);
p1->SetBorderSize(0.0);
p1->SetTopMargin(tmar*hpanel/(hpanel+tmar));
p1->SetRightMargin(0.0);
p1->SetLeftMargin(lmar*wpanel/(wpanel+lmar));
p1->SetBottomMargin(0.0);
p1->SetFrameBorderMode(0);
hx->Draw();
gx->Draw("P");

cbeam->cd(2);
TPad *p2;
p2 = gPad;
p2->SetPad("p2","",(wpanel+lmar)/width,(hpanel+bmar)/height,1.0,1.0);
p2->SetFillColor(0);
p2->SetBorderMode(0);
p2->SetBorderSize(0.0);
p2->SetTopMargin(tmar*hpanel/(hpanel+tmar));
p2->SetRightMargin(rmar*wpanel/(wpanel+rmar));
p2->SetLeftMargin(0.0);
p2->SetBottomMargin(0.0);
p2->SetFrameBorderMode(0);
hy->Draw();
gy->Draw("P");

cbeam->cd(3);
TPad *p3;
p3 = gPad;
p3->SetPad("p3","",0.0,0.0,(wpanel+lmar)/width,(hpanel+bmar)/height);
p3->SetFillColor(0);
p3->SetBorderMode(0);
p3->SetBorderSize(0.0);
p3->SetTopMargin(0.0);
p3->SetRightMargin(0.0);
p3->SetLeftMargin(lmar*wpanel/(wpanel+lmar));
p3->SetBottomMargin(bmar*hpanel/(hpanel+bmar));
p3->SetFrameBorderMode(0);
hdx->Draw();
gdx->Draw("P");

cbeam->cd(4);
TPad *p4;
p4 = gPad;
p4->SetPad("p4","",(wpanel+lmar)/width,0.0,1.0,(hpanel+bmar)/height);
p4->SetFillColor(0);
p4->SetBorderMode(0);
p4->SetBorderSize(0.0);
p4->SetTopMargin(0.0);
p4->SetRightMargin(rmar*wpanel/(wpanel+rmar));
p4->SetLeftMargin(0.0);
p4->SetBottomMargin(bmar*hpanel/(hpanel+bmar));
p4->SetFrameBorderMode(0);
hdy->Draw();
gdy->Draw("P");
}