#include <fstream>
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <TObject.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TStyle.h>
#include <TPad.h>
#include <TLegend.h>
#include <TF1.h>

using namespace std;

int main(int argc, char* argv[]){

TFile* f0 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.390GeV.root");
TFile* f1 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.395GeV.root");
TFile* f2 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.400GeV.root");
TFile* f3 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.405GeV.root");
TFile* f4 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.410GeV.root");
TFile* f5 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.420GeV.root");
TFile* f6 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.430GeV.root");
TFile* f7 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.440GeV.root");
TFile* f8 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.445GeV.root");
TFile* f9 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.460GeV.root");
TFile* f10 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.480GeV.root");
TFile* f11 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.500GeV.root");
TFile* f12 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.520GeV.root");
TFile* f13 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.540GeV.root");
TFile* f14 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.551GeV.root");
TFile* f15 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.560GeV.root");
TFile* f16 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.570GeV.root");
TFile* f17 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.580GeV.root");
TFile* f18 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.590GeV.root");
TFile* f19 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.599GeV.root");
TFile* f20 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.600GeV.root");
TFile* f21 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.610GeV.root");
TFile* f22 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.620GeV.root");
TFile* f23 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.630GeV.root");
TFile* f24 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.644GeV.root");
TFile* f25 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.729GeV.root");
TFile* f26 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.808GeV.root");
TFile* f27 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.881GeV.root");
TFile* f28 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_0.954GeV.root");
TFile* f29 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_1.100GeV.root");
TFile* f30 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_1.219GeV.root");
TFile* f31 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_1.300GeV.root");
TFile* f32 = new TFile("/home/malakhov/CheckGeom/Protons/ProtonSeparation2/ParDirCalc_Protons_1.400GeV.root");
TH1F* h0 = (TH1F*)f0->Get("Proton");
TH1F* h1 = (TH1F*)f1->Get("Proton");
TH1F* h2 = (TH1F*)f2->Get("Proton");
TH1F* h3 = (TH1F*)f3->Get("Proton");
TH1F* h4 = (TH1F*)f4->Get("Proton");
TH1F* h5 = (TH1F*)f5->Get("Proton");
TH1F* h6 = (TH1F*)f6->Get("Proton");
TH1F* h7 = (TH1F*)f7->Get("Proton");
TH1F* h8 = (TH1F*)f8->Get("Proton");
TH1F* h9 = (TH1F*)f9->Get("Proton");
TH1F* h10 = (TH1F*)f10->Get("Proton");
TH1F* h11 = (TH1F*)f11->Get("Proton");
TH1F* h12 = (TH1F*)f12->Get("Proton");
TH1F* h13 = (TH1F*)f13->Get("Proton");
TH1F* h14 = (TH1F*)f14->Get("Proton");
TH1F* h15 = (TH1F*)f15->Get("Proton");
TH1F* h16 = (TH1F*)f16->Get("Proton");
TH1F* h17 = (TH1F*)f17->Get("Proton");
TH1F* h18 = (TH1F*)f18->Get("Proton");
TH1F* h19 = (TH1F*)f19->Get("Proton");
TH1F* h20 = (TH1F*)f20->Get("Proton");
TH1F* h21 = (TH1F*)f21->Get("Proton");
TH1F* h22 = (TH1F*)f22->Get("Proton");
TH1F* h23 = (TH1F*)f23->Get("Proton");
TH1F* h24 = (TH1F*)f24->Get("Proton");
TH1F* h25 = (TH1F*)f25->Get("Proton");
TH1F* h26 = (TH1F*)f26->Get("Proton");
TH1F* h27 = (TH1F*)f27->Get("Proton");
TH1F* h28 = (TH1F*)f28->Get("Proton");
TH1F* h29 = (TH1F*)f29->Get("Proton");
TH1F* h30 = (TH1F*)f30->Get("Proton");
TH1F* h31 = (TH1F*)f31->Get("Proton");
TH1F* h32 = (TH1F*)f32->Get("Proton");
Float_t x[33]; Float_t xe[33]; Float_t xr[33]; Float_t xre[33]; Float_t G[33]; Float_t Ge[33]; Float_t xd[33];
for(Int_t i=0;i<33;i++)xe[i]=0;
x[0]=0.390; x[1]=0.395; x[2]=0.400; x[3]=0.405; x[4]=0.410; x[5]=0.420; x[6]=0.430; x[7]=0.440; x[8]=0.444; x[9]=0.460;
x[10]=0.480; x[11]=0.500; x[12]=0.520; x[13]=0.540; x[14]=0.551; x[15]=0.560; x[16]=0.570; x[17]=0.580; x[18]=0.590;
x[19]=0.599; x[20]=0.600; x[21]=0.610; x[22]=0.620; x[23]=0.630; x[24]=0.644; x[25]=0.729; x[26]=0.808; x[27]=0.883;
x[28]=0.954; x[29]=1.100; x[30]=1.218; x[31]=1.300; x[32]=1.400;
xr[0]=h0->GetMean(); xre[0]=h0->GetMeanError();
xr[1]=h1->GetMean(); xre[1]=h1->GetMeanError();
xr[2]=h2->GetMean(); xre[2]=h2->GetMeanError();
xr[3]=h3->GetMean(); xre[3]=h3->GetMeanError();
xr[4]=h4->GetMean(); xre[4]=h4->GetMeanError();
xr[5]=h5->GetMean(); xre[5]=h5->GetMeanError();
xr[6]=h6->GetMean(); xre[6]=h6->GetMeanError();
xr[7]=h7->GetMean(); xre[7]=h7->GetMeanError();
xr[8]=h8->GetMean(); xre[8]=h8->GetMeanError();
xr[9]=h9->GetMean(); xre[9]=h9->GetMeanError();
xr[10]=h10->GetMean(); xre[10]=h10->GetMeanError();
xr[11]=h11->GetMean(); xre[11]=h11->GetMeanError();
xr[12]=h12->GetMean(); xre[12]=h12->GetMeanError();
xr[13]=h13->GetMean(); xre[13]=h13->GetMeanError();
xr[14]=h14->GetMean(); xre[14]=h14->GetMeanError();
xr[15]=h15->GetMean(); xre[15]=h15->GetMeanError();
xr[16]=h16->GetMean(); xre[16]=h16->GetMeanError();
xr[17]=h17->GetMean(); xre[17]=h17->GetMeanError();
xr[18]=h18->GetMean(); xre[18]=h18->GetMeanError();
xr[19]=h19->GetMean(); xre[19]=h19->GetMeanError();
xr[20]=h20->GetMean(); xre[20]=h20->GetMeanError();
xr[21]=h21->GetMean(); xre[21]=h21->GetMeanError();
xr[22]=h22->GetMean(); xre[22]=h22->GetMeanError();
xr[23]=h23->GetMean(); xre[23]=h23->GetMeanError();
xr[24]=h24->GetMean(); xre[24]=h24->GetMeanError();
xr[25]=h25->GetMean(); xre[25]=h25->GetMeanError();
xr[26]=h26->GetMean(); xre[26]=h26->GetMeanError();
xr[27]=h27->GetMean(); xre[27]=h27->GetMeanError();
xr[28]=h28->GetMean(); xre[28]=h28->GetMeanError();
xr[29]=h29->GetMean(); xre[29]=h29->GetMeanError();
xr[30]=h30->GetMean(); xre[30]=h30->GetMeanError();
xr[31]=h31->GetMean(); xre[31]=h31->GetMeanError();
xr[32]=h32->GetMean(); xre[32]=h32->GetMeanError();
G[0]=3.14*3600*h0->GetEntries()/400000; Ge[0]=3.14*3600*sqrt(h0->GetEntries())/400000;
G[1]=3.14*3600*h1->GetEntries()/400000; Ge[1]=3.14*3600*sqrt(h1->GetEntries())/400000;
G[2]=3.14*3600*h2->GetEntries()/400000; Ge[2]=3.14*3600*sqrt(h2->GetEntries())/400000;
G[3]=3.14*3600*h3->GetEntries()/400000; Ge[3]=3.14*3600*sqrt(h3->GetEntries())/400000;
G[4]=3.14*3600*h4->GetEntries()/400000; Ge[4]=3.14*3600*sqrt(h4->GetEntries())/400000;
G[5]=3.14*3600*h5->GetEntries()/400000; Ge[5]=3.14*3600*sqrt(h5->GetEntries())/400000;
G[6]=3.14*3600*h6->GetEntries()/400000; Ge[6]=3.14*3600*sqrt(h6->GetEntries())/400000;
G[7]=3.14*3600*h7->GetEntries()/400000; Ge[7]=3.14*3600*sqrt(h7->GetEntries())/400000;
G[8]=3.14*3600*h8->GetEntries()/400000; Ge[8]=3.14*3600*sqrt(h8->GetEntries())/400000;
G[9]=3.14*3600*h9->GetEntries()/400000; Ge[9]=3.14*3600*sqrt(h9->GetEntries())/400000;
G[10]=3.14*3600*h10->GetEntries()/400000; Ge[10]=3.14*3600*sqrt(h10->GetEntries())/400000;
G[11]=3.14*3600*h11->GetEntries()/400000; Ge[11]=3.14*3600*sqrt(h11->GetEntries())/400000;
G[12]=3.14*3600*h12->GetEntries()/400000; Ge[12]=3.14*3600*sqrt(h12->GetEntries())/400000;
G[13]=3.14*3600*h13->GetEntries()/400000; Ge[13]=3.14*3600*sqrt(h13->GetEntries())/400000;
G[14]=3.14*3600*h14->GetEntries()/400000; Ge[14]=3.14*3600*sqrt(h14->GetEntries())/400000;
G[15]=3.14*3600*h15->GetEntries()/400000; Ge[15]=3.14*3600*sqrt(h15->GetEntries())/400000;
G[16]=3.14*3600*h16->GetEntries()/400000; Ge[16]=3.14*3600*sqrt(h16->GetEntries())/400000;
G[17]=3.14*3600*h17->GetEntries()/400000; Ge[17]=3.14*3600*sqrt(h17->GetEntries())/400000;
G[18]=3.14*3600*h18->GetEntries()/400000; Ge[18]=3.14*3600*sqrt(h18->GetEntries())/400000;
G[19]=3.14*3600*h19->GetEntries()/400000; Ge[19]=3.14*3600*sqrt(h19->GetEntries())/400000;
G[20]=3.14*3600*h20->GetEntries()/400000; Ge[20]=3.14*3600*sqrt(h20->GetEntries())/400000;
G[21]=3.14*3600*h21->GetEntries()/400000; Ge[21]=3.14*3600*sqrt(h21->GetEntries())/400000;
G[22]=3.14*3600*h22->GetEntries()/400000; Ge[22]=3.14*3600*sqrt(h22->GetEntries())/400000;
G[23]=3.14*3600*h23->GetEntries()/400000; Ge[23]=3.14*3600*sqrt(h23->GetEntries())/400000;
G[24]=3.14*3600*h24->GetEntries()/400000; Ge[24]=3.14*3600*sqrt(h24->GetEntries())/400000;
G[25]=3.14*3600*h25->GetEntries()/400000; Ge[25]=3.14*3600*sqrt(h25->GetEntries())/400000;
G[26]=3.14*3600*h26->GetEntries()/400000; Ge[26]=3.14*3600*sqrt(h26->GetEntries())/400000;
G[27]=3.14*3600*h27->GetEntries()/400000; Ge[27]=3.14*3600*sqrt(h27->GetEntries())/400000;
G[28]=3.14*3600*h28->GetEntries()/400000; Ge[28]=3.14*3600*sqrt(h28->GetEntries())/400000;
G[29]=3.14*3600*h29->GetEntries()/400000; Ge[29]=3.14*3600*sqrt(h29->GetEntries())/400000;
G[30]=3.14*3600*h30->GetEntries()/400000; Ge[30]=3.14*3600*sqrt(h30->GetEntries())/400000;
G[31]=3.14*3600*h31->GetEntries()/400000; Ge[31]=3.14*3600*sqrt(h31->GetEntries())/400000;
G[32]=3.14*3600*h32->GetEntries()/400000; Ge[32]=3.14*3600*sqrt(h32->GetEntries())/400000;
for(Int_t i=0;i<33;i++)xd[i]=x[i]-xr[i];
TGraphErrors* Geom = new TGraphErrors(33,x,G,xe,Ge);
TGraphErrors* dif1 = new TGraphErrors(33,x,xd,xe,xre);
TGraphErrors* dif2 = new TGraphErrors(33,xr,xd,xre,xre);
TGraphErrors* dif3 = new TGraphErrors(33,xr,x,xre,xe);
TCanvas* c1 = new TCanvas("c1","c1",1200,800);
Geom->Draw("AP");
Geom->SetName("Geom");
c1->SaveAs((TString)argv[1]+"Geom_33points.png");
delete c1;
TCanvas* c2 = new TCanvas("c2","c2",1200,800);
dif1->Draw("AP");
dif1->SetName("dif1");
c2->SaveAs((TString)argv[1]+"Rshift1_33points.png");
delete c2;
TCanvas* c3 = new TCanvas("c3","c3",1200,800);
dif2->Draw("AP");
dif2->SetName("dif2");
c3->SaveAs((TString)argv[1]+"Rshift2_33points.png");
delete c3;
TCanvas* c4 = new TCanvas("c4","c4",1200,800);
dif3->Draw("AP");
dif3->SetName("dif3");
c4->SaveAs((TString)argv[1]+"Rshift2_33points.png");
delete c4;
TFile* out = new TFile((TString)argv[1]+"CheckGeom33points.root","recreate");
Geom->Write();
dif1->Write();
dif2->Write();
dif3->Write();
out->Close();

}
