{

// Root macro to produce difference plots between scaled OPERA field and measured field.

Double_t lmar = 0.3;
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 TScale = hpanel/(hpanel+tmar);
Double_t BScale = hpanel/(hpanel+bmar);
Double_t LScale = wpanel/(wpanel+lmar);
Double_t RScale = wpanel/(wpanel+rmar);
Double_t Tlbl;
Tlbl = 1.2*lbl;

FILE *fHigh;
fHigh = fopen("2002-04-26-15-52-59.dat","r");

FILE *fLow;
fLow = fopen("2002-04-24-17-39-16.dat","r");

FILE *fOPERA;
fOPERA = fopen("2002-04-19-16-56-15.dat","r");

Int_t istat;
char CommentLine[100];
Int_t nz = 261;

// Read the comment lines
istat = fscanf(fHigh,"%s",&CommentLine);
istat = fscanf(fLow,"%s",&CommentLine);
istat = fscanf(fOPERA,"%s",&CommentLine);

Float_t phi[24];

Float_t BHAxis[nz];
Float_t BHR16[nz];
Float_t BLAxis[nz];
Float_t BLR16[nz];
Float_t BOAxis[nz];
Float_t BOR16[nz];
Float_t zHB[nz];
Float_t zLB[nz];
Float_t zOB[nz];

Int_t iEntry;
Float_t z;
Float_t b1;
Float_t b2;
Float_t b3;
Float_t b4;
Float_t b5;
Float_t b6;
Float_t b7;


for(Int_t iz=0;iz<nz;iz++){
  BHAxis[iz] = 0.0;
  BHR16[iz] = 0.0;
  for(Int_t iphi=0;iphi<24;iphi++){
    istat = fscanf(fHigh,"%i",&iEntry);
    istat = fscanf(fHigh,"%f",&z);
    istat = fscanf(fHigh,"%f",&phi[iphi]);
    istat = fscanf(fHigh,"%f %f %f %f %f %f %f",&b1,&b2,&b3,&b4,&b5,&b5,&b6,&b7);
    //    printf("z = %f, b0 = %f, b16 = %f\n",z,b1,b5);
    if(iphi == 0){
      zHB[iz] = z;
    }
    else{
      if(fabs(z - zHB[iz]) > 1.0){
	printf("z doesn't match for different phi, high field\n");
      }
    }
    BHAxis[iz] = BHAxis[iz] + b1/24.0;
    BHR16[iz] = BHR16[iz] + b5/24.0;
  }
}

for(Int_t iz=0;iz<nz;iz++){
  BLAxis[iz] = 0.0;
  BLR16[iz] = 0.0;
  for(Int_t iphi=0;iphi<24;iphi++){
    istat = fscanf(fLow,"%i",&iEntry);
    istat = fscanf(fLow,"%f",&z);
    istat = fscanf(fLow,"%f",&phi[iphi]);
    istat = fscanf(fLow,"%f %f %f %f %f %f %f",&b1,&b2,&b3,&b4,&b5,&b5,&b6,&b7);
    if(iphi == 0){
      zLB[iz] = z;
    }
    else{
      if(fabs(z - zLB[iz]) > 1.0){
	printf("z doesn't match for different phi, low field\n");
      }
    }
    BLAxis[iz] = BLAxis[iz] + b1/24.0;
    BLR16[iz] = BLR16[iz] + b5/24.0;
  }
}

for(Int_t iz=0;iz<nz;iz++){
  BOAxis[iz] = 0.0;
  BOR16[iz] = 0.0;
  for(Int_t iphi=0;iphi<24;iphi++){
    istat = fscanf(fOPERA,"%i",&iEntry);
    istat = fscanf(fOPERA,"%f",&z);
    istat = fscanf(fOPERA,"%f",&phi[iphi]);
    istat = fscanf(fOPERA,"%f %f %f %f %f %f %f",&b1,&b2,&b3,&b4,&b5,&b5,&b6,&b7);
    if(iphi == 0){
      zOB[iz] = z;
    }
    else{
      if(fabs(z - zOB[iz]) > 1.0){
	printf("z doesn't match for different phi, OPERA\n");
      }
    }
    BOAxis[iz] = BOAxis[iz] + b1/24.0;
    BOR16[iz] = BOR16[iz] + b5/24.0;
  }
}

Float_t DiffHAxis[nz];
Float_t DiffHR16[nz];
Float_t DiffLAxis[nz];
Float_t DiffLR16[nz];

// Calculate differences
for(Int_t iz=0;iz<nz;iz++){
  if(fabs(zHB[iz] - zOB[iz]) > 1.0){
    printf("z doesn't match between high field and OPERA\n");
  }
  // High field along axis
  DiffHAxis[iz] = BHAxis[iz] - 0.98003*BOAxis[iz];
  //  printf("BH = %f,  BO = %f, BDiff = %f\n",BHAxis[iz],0.98003*BOAxis[iz],DiffHAxis[iz]);
  // High field at r = 16.5 cm
  DiffHR16[iz] = BHR16[iz] - 0.98003*BOR16[iz];
  if(fabs(zLB[iz] - zOB[iz]) > 1.0){
    printf("z doesn't match between low field and OPERA\n");
  }
  // Low field along axis
  DiffLAxis[iz] = BLAxis[iz] - 1.02002*BOAxis[iz];
  //Low field at r = 16.5 cm
  DiffLR16[iz] = BLR16[iz] - 1.02002*BOR16[iz];
}

TH2F *hHAxis;
hHAxis = new TH2F("hHAxis","",14,-700.,700.,11,-0.001,0.010);
hHAxis->SetStats(kFALSE);
hHAxis->GetYaxis()->SetTitle("#Delta B_{z} (B=1.96 T)");
hHAxis->GetYaxis()->CenterTitle();
hHAxis->GetYaxis()->SetTitleSize(Tlbl*TScale);
hHAxis->GetYaxis()->SetLabelSize(lbl*TScale);
hHAxis->GetXaxis()->SetLabelSize(0.0);
hHAxis->GetYaxis()->SetNdivisions(206,kTRUE);
hHAxis->GetXaxis()->SetNdivisions(205,kTRUE);
hHAxis->GetYaxis()->SetTitleOffset(0.9/TScale);

TH2F *hHR16;
hHR16 = new TH2F("hHR16","",14,-700.,700.,11,-0.001,0.010);
hHR16->SetStats(kFALSE);
hHR16->GetYaxis()->SetLabelSize(0.0);
hHR16->GetXaxis()->SetLabelSize(0.0);
hHR16->GetYaxis()->SetNdivisions(206,kTRUE);
hHR16->GetXaxis()->SetNdivisions(205,kTRUE);

TH2F *hLAxis;
hLAxis = new TH2F("hLAxis","",14,-700.,700.,9,-0.008,0.001);
hLAxis->SetStats(kFALSE);
hLAxis->GetYaxis()->SetTitle("#Delta B_{z} (B=2.04 T)");
hLAxis->GetYaxis()->CenterTitle();
hLAxis->GetYaxis()->SetTitleSize(Tlbl*BScale);
hLAxis->GetYaxis()->SetLabelSize(lbl*BScale);
hLAxis->GetXaxis()->SetTitle("z (mm) along z axis");
hLAxis->GetXaxis()->CenterTitle();
hLAxis->GetXaxis()->SetTitleSize(Tlbl*BScale);
hLAxis->GetXaxis()->SetLabelSize(lbl*BScale);
hLAxis->GetYaxis()->SetNdivisions(205,kTRUE);
hLAxis->GetXaxis()->SetNdivisions(205,kTRUE);
hLAxis->GetXaxis()->SetTitleOffset(0.7/BScale);
hLAxis->GetYaxis()->SetTitleOffset(0.9/BScale);

TH2F *hLR16;
hLR16 = new TH2F("hLR16","",14,-700.,700.,9,-0.008,0.001);
hLR16->SetStats(kFALSE);
hLR16->GetYaxis()->SetLabelSize(0.0);
hLR16->GetXaxis()->SetTitle("z (mm) at r = 16.5 cm");
hLR16->GetXaxis()->CenterTitle();
hLR16->GetXaxis()->SetTitleSize(Tlbl*BScale);
hLR16->GetXaxis()->SetLabelSize(lbl*BScale);
hLR16->GetYaxis()->SetNdivisions(205,kTRUE);
hLR16->GetXaxis()->SetNdivisions(205,kTRUE);
hLR16->GetXaxis()->SetTitleOffset(0.7/BScale);

TGraph *gHAxis;
gHAxis = new TGraph(nz,zHB,DiffHAxis);
gHAxis->SetMarkerStyle(8);
gHAxis->SetMarkerColor(4);

TGraph *gHR16;
gHR16 = new TGraph(nz,zHB,DiffHR16);
gHR16->SetMarkerStyle(8);
gHR16->SetMarkerColor(4);

TGraph *gLAxis;
gLAxis = new TGraph(nz,zLB,DiffLAxis);
gLAxis->SetMarkerStyle(8);
gLAxis->SetMarkerColor(4);

TGraph *gLR16;
gLR16 = new TGraph(nz,zLB,DiffLR16);
gLR16->SetMarkerStyle(8);
gLR16->SetMarkerColor(4);

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

cB->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);
hHAxis->Draw();
gHAxis->Draw("P");

cB->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);
hHR16->Draw();
gHR16->Draw("P");

cB->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);
hLAxis->Draw();
gLAxis->Draw("P");

cB->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);
hLR16->Draw();
gLR16->Draw("P");

}