#include <TH1F.h>
#include <TH2F.h>
#include <TFile.h>
#include <TF1.h>
#include <TProfile.h>
#include <TTree.h>
#include <TGraphErrors.h>
#include <TCanvas.h>
#include <unistd.h>
void FgdRes(void){
  bool useGauss = true;
  bool divideAngle = false;
  double lowAngleCut = -1;
  double highAngleCut = 99;
  int maxHitCntPerLayer = 2;
  int hitCntPerLayer = 2;
  int maxHitCntPerEvent = 2000;
  double fitRangeHalf = 5;
  double angleConst = 0.02;
  double fgd1XInt = 120.85;
  double fgd1YInt = 130.65;
  double fgd2XInt = 1478.85;
  double fgd2YInt = 1488.65;
  double fgd1SepDist = 22.25;
  double fgd2SepDist = 52.25;
  double fgd1Center = 281.5;
  double fgd2Center = 1640.5;
  gStyle->SetOptFit(11111111);
  TChain* chain = new TChain("fgdintaligntree");
  /*
  chain->AddFile("../cmt/50kFGDMisalignedP0DReco1.root", 0);
  //chain->AddFile("../cmt/50kFGDMisalignedP0DReco2.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedP0DReco3.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedP0DReco4.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedP0DReco5.root", 0);
  */
   
  chain->AddFile("../cmt/50kFGDMisalignedRECONSBCATReco1New.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedRECONSBCATReco2New.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedRECONSBCATReco3New.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedRECONSBCATReco4New.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedRECONSBCATReco5New.root", 0);
  /*
  chain->AddFile("../cmt/50kFGDMisalignedRECONSBCATReco1.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedRECONSBCATReco2.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedRECONSBCATReco3.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedRECONSBCATReco4.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedRECONSBCATReco5.root", 0);
  */
  /*
  chain->AddFile("../cmt/50kFGDMisalignedRadonReco1.root", 0);
  //chain->AddFile("../cmt/50kFGDMisalignedRadonReco2.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedRadonReco3.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedRadonReco4.root", 0);
  chain->AddFile("../cmt/50kFGDMisalignedRadonReco5.root", 0);
  */
  //chain->AddFile("../cmt/50kFGDMisalignedReco4.root", 0);
  //chain->AddFile("../cmt/50kFGDDesignReco0.root", 0);
  //chain->AddFile("../cmt/50kFGDMisalignedReco0.root", 0);
  //chain->AddFile("../cmt/50kFGDMisalignedRECONSBCATReco5.root", 0);
  
  int layerNum;
  int evtHitCnt;
  int layerHitCnt[44];
  double layerFGD1X[15];
  double layerFGD1Y[15];
  double layerFGD2X[7];
  double layerFGD2Y[7];
  double HitCnt0[44];
  double HitCnt1[44];
  double HitCnt2[44];
  double HitCnt3[44];
  
  double layerErrFGD1[15];
  double layerErrFGD2[7];
  double layer[44];
  double layerErr[44];
  double sigmaOrRMS[44];
  double sigmaOrRMSPos[44];
  double sigmaOrRMSNeg[44];
  
  double res[44];
  double meanResDiff[44];
  //double meanResDiff[44];
  // mean residual and error
  double meanResFGD1X[15]; 
  double meanResFGD1Y[15];
  double meanResFGD2X[7];
  double meanResFGD2Y[7];
  double meanResErrFGD1X[15]; 
  double meanResErrFGD1Y[15];
  double meanResErrFGD2X[7];
  double meanResErrFGD2Y[7];

  double meanResNegAngle[44];
  double meanResNegAngleErr[44];
  double meanResPosAngle[44];
  double meanResPosAngleErr[44];
  double fgd1TrackAngleX;
  double fgd1TrackAngleY;
  double fgd2TrackAngleX;
  double fgd2TrackAngleY;

  bool fgd1XGood;
  bool fgd1YGood;
  bool fgd2XGood;
  bool fgd2YGood;
 
  double fgdIntTransConst[44];
  fgdIntTransConst[0] = -0.379351;
  fgdIntTransConst[1] = -0.115179;
  fgdIntTransConst[2] = 0.210465;
  fgdIntTransConst[3] = 0.161571;
  fgdIntTransConst[4] = -0.0505728;
  fgdIntTransConst[5] = 0.0535972; 
  fgdIntTransConst[6] = 0.0305484; 
  fgdIntTransConst[7] = -0.199944;
  fgdIntTransConst[8] = -0.0697529;
  fgdIntTransConst[9] = -0.270023;
  fgdIntTransConst[10] = -0.129371; 
  fgdIntTransConst[11] = 0.0471496; 
  fgdIntTransConst[12] = 0.255003; 
  fgdIntTransConst[13] = -0.13411; 
  fgdIntTransConst[14] = 0.214838; 
  fgdIntTransConst[15] = 0.0558979; 
  fgdIntTransConst[16] = 0.266153; 
  fgdIntTransConst[17] = 0.555704; 
  fgdIntTransConst[18] = 0.0547172; 
  fgdIntTransConst[19] = 0.351819; 
  fgdIntTransConst[20] = 0.00473694; 
  fgdIntTransConst[21] = 0.786871; 
  fgdIntTransConst[22] = 0.0380896; 
  fgdIntTransConst[23] = -0.401259; 
  fgdIntTransConst[24] = -0.245828;
  fgdIntTransConst[25] = -0.414316; 
  fgdIntTransConst[26] = -0.124361; 
  fgdIntTransConst[27] = -0.0840803; 
  fgdIntTransConst[28] = -0.0219941;
  fgdIntTransConst[29] = -0.251145; 
  fgdIntTransConst[30] = -0.0640535; 
  fgdIntTransConst[31] = -0.0783846; 
  fgdIntTransConst[32] = -0.143142; 
  fgdIntTransConst[33] = 0.408897;
  fgdIntTransConst[34] = 0.221593; 
  fgdIntTransConst[35] = -0.477095;
  fgdIntTransConst[36] = 0.20223; 
  fgdIntTransConst[37] = 0.235765;
  fgdIntTransConst[38] = -0.382516; 
  fgdIntTransConst[39] = 0.0894779; 
  fgdIntTransConst[40] = 0.401933; 
  fgdIntTransConst[41] = -0.134181;
  fgdIntTransConst[42] = -0.235573; 
  fgdIntTransConst[43] = 0.0246083; 

  double fgdIntSysError[44];
  fgdIntSysError[0] = -0.019;
  fgdIntSysError[1] = -0.01;
  fgdIntSysError[2] = -0.005;
  fgdIntSysError[3] = 0.016;
  fgdIntSysError[4] = 0.01;
  fgdIntSysError[5] = 0.024; 
  fgdIntSysError[6] = 0.016; 
  fgdIntSysError[7] = -0.008;
  fgdIntSysError[8] = -0.009;
  fgdIntSysError[9] = 0.007;
  fgdIntSysError[10] = 0.01; 
  fgdIntSysError[11] = 0.02; 
  fgdIntSysError[12] = 0.006; 
  fgdIntSysError[13] = 0.004; 
  fgdIntSysError[14] = -0.012; 
  fgdIntSysError[15] = -0.002; 
  fgdIntSysError[16] = 0.02; 
  fgdIntSysError[17] = 0.012; 
  fgdIntSysError[18] = 0.01; 
  fgdIntSysError[19] = 0.012; 
  fgdIntSysError[20] = 0.004; 
  fgdIntSysError[21] = 0.034; 
  fgdIntSysError[22] = -0.009; 
  fgdIntSysError[23] = 0.007; 
  fgdIntSysError[24] = 0.01;
  fgdIntSysError[25] = 0.01; 
  fgdIntSysError[26] = 0.007; 
  fgdIntSysError[27] = 0.0047; 
  fgdIntSysError[28] = 0.013;
  fgdIntSysError[29] = 0.02; 
  fgdIntSysError[30] = -0.028; 
  fgdIntSysError[31] = -0.0004;
  fgdIntSysError[32] = -0.009; 
  fgdIntSysError[33] = 0.02;
  fgdIntSysError[34] = 0.003; 
  fgdIntSysError[35] = 0.009;
  fgdIntSysError[36] = -0.004; 
  fgdIntSysError[37] = 0.001;
  fgdIntSysError[38] = -0.017; 
  fgdIntSysError[39] = 0.01; 
  fgdIntSysError[40] = 0.03; 
  fgdIntSysError[41] = 0.0005;
  fgdIntSysError[42] = -0.021; 
  fgdIntSysError[43] = 0.027; 

  double chisqFGD1X = 0;
  double chisqFGD1Y = 0;
  double chisqFGD2X = 0;
  double chisqFGD2Y = 0;

  for (int i = 0; i < 15; i++)
    {
      layerFGD1X[i] = i*2;
      layerFGD1Y[i] = i*2+1;
      layerErrFGD1[i] = 0;
    }
  for (int j = 0; j < 7; j++)
    {
      layerFGD2X[j] = j*2+30;
      layerFGD2Y[j] = j*2+31;
      layerErrFGD2[j] = 0;
    }
  chain->SetBranchAddress("layer", &layerNum);
  chain->SetBranchAddress("residual", res);
  chain->SetBranchAddress("layer hit count", layerHitCnt); 
  chain->SetBranchAddress("event hit count", &evtHitCnt);
  chain->SetBranchAddress("FGD1X track angle", &fgd1TrackAngleX);
  chain->SetBranchAddress("FGD1Y track angle", &fgd1TrackAngleY);
  chain->SetBranchAddress("FGD2X track angle", &fgd2TrackAngleX);
  chain->SetBranchAddress("FGD2Y track angle", &fgd2TrackAngleY);
  Long64_t nentries = chain->GetEntries();
  std::cout<<"Entries: "<<nentries<<std::endl;
  TCanvas* resCanvasFGD1X = new TCanvas("residual histograms of FGD1X layers", "residual histograms of FGD1X layers", 1000, 700);
  TCanvas* resCanvasFGD1Y = new TCanvas("residual histograms of FGD1Y layers", "residual histograms of FGD1Y layers", 1000, 700);
  TCanvas* resCanvasFGD2X = new TCanvas("residual histograms of FGD2X layers", "residual histograms of FGD2X layers", 1000, 700);
  TCanvas* resCanvasFGD2Y = new TCanvas("residual histograms of FGD2Y layers", "residual histograms of FGD2Y layers", 1000, 700);
 

  TH1F* resHist[44];
  TH1F* resHistNegAngle[44];
  TH1F* resHistPosAngle[44];
  TH1F* angHistFGD1X = new TH1F("angular distribution of FGD1X tracks", "angular distribution of FGD1X tracks", 180, -90, 90);
  TH1F* angHistFGD1Y = new TH1F("angular distribution of FGD1Y tracks", "angular distribution of FGD1Y tracks", 180, -90, 90);
  TH1F* angHistFGD2X = new TH1F("angular distribution of FGD2X tracks", "angular distribution of FGD2X tracks", 180, -90, 90);
  TH1F* angHistFGD2Y = new TH1F("angular distribution of FGD2Y tracks", "angular distribution of FGD2Y tracks", 180, -90, 90);
  TH1F* resDiffRMSHist = new TH1F("residual difference distribution", "residual difference distribution", 50, -0.1, 0.1);
  TH1F* resDiffNegRMSHist = new TH1F("negative angle residual difference distribution", "negative angle residual difference distribution", 50, -0.1, 0.1);
  TH1F* resDiffPosRMSHist = new TH1F("positive angle residual difference distribution", "positive angle residual difference distribution", 50, -0.1, 0.1);
  TCanvas* graphCanvas;
  TCanvas* resDiffGraphCanvas;
  //TH2F* hitCntHist = new TH2F("hit count distribution", "hit count distribution", 44, -0.5, 43.5, 5, -0.5, 4.5);
  TGraphErrors* HitCnt0Graph;
  TGraphErrors* HitCnt1Graph;
  TGraphErrors* HitCnt2Graph;
  TGraphErrors* HitCnt3Graph;
  TGraphErrors* resGraphFGD1X;
  TGraphErrors* resGraphFGD1Y;
  TGraphErrors* resGraphFGD2X;
  TGraphErrors* resGraphFGD2Y;
  TGraphErrors* resDiffGraph;
  TGraphErrors* sigmaOrRMSGraph;
  TGraphErrors* sigmaOrRMSPosGraph;
  TGraphErrors* sigmaOrRMSNegGraph;
  TGraphErrors* resNegAngleGraph;
  TGraphErrors* resPosAngleGraph;

  TGraphErrors* refLine;
  TGraphErrors* origTrans;
  TGraphErrors* sigmaPattern;
  TGraphErrors* fgdDivider;
  char resGName[100];
  char resGNegName[100];
  char resGPosName[100];

  for (int i = 0; i < 44; i++)
    {
      sprintf(resGName, "residuals of layer %d", i);
      sprintf(resGNegName, "residuals of layer %d from negative angled tracks", i);
      sprintf(resGPosName, "residuals of layer %d from positive angled tracks", i);
      //std::cout<<resGName<<std::endl;
      resHist[i] = new TH1F(resGName, resGName, 201, -20.25, 20.25);
      resHistNegAngle[i] = new TH1F(resGNegName, resGNegName, 201, -20.25, 20.25);
      resHistPosAngle[i] = new TH1F(resGPosName, resGPosName, 201, -20.25, 20.25);
    
      //fgdIntTransConst[i]*= -1;
      fgdIntSysError[i]*= -1;
      layer[i] = i;
      layerErr[i] = 0;
      HitCnt0[i] = 0;
      HitCnt1[i] = 0;
      HitCnt2[i] = 0;
      HitCnt3[i] = 0;
    }
  
  for (Long64_t iev = 0; iev < nentries; iev++)
    {
      fgd1XGood = true;
      fgd1YGood = true;
      fgd2XGood = true;
      fgd2YGood = true;
      chain->GetEntry(iev);
      if (evtHitCnt < maxHitCntPerEvent)
  {
    angHistFGD1X->Fill(fgd1TrackAngleX);
    angHistFGD1Y->Fill(fgd1TrackAngleY);
    angHistFGD2X->Fill(fgd2TrackAngleX);
    angHistFGD2Y->Fill(fgd2TrackAngleY);
    for (int i = 0; i < 44; i++)
      {
        if (layerHitCnt[i]==0)
    HitCnt0[i]++;
        else if (layerHitCnt[i]==1)
    HitCnt1[i]++;
        else if (layerHitCnt[i]==2)
    HitCnt2[i]++;
        else if (layerHitCnt[i]==3)
    HitCnt3[i]++;
    //hitCntHist->Fill(i, layerHitCnt[i]);
        int cat = FGDCategory(i);
        //std::cout<<cat<<std::endl;
        if (layerHitCnt[i] == 0 || layerHitCnt[i] > maxHitCntPerLayer)
    if (cat == 0)
      fgd1XGood = false;
           else if (cat == 1)
      fgd1YGood = false;
           else if (cat == 2)
      fgd2XGood = false;
           else 
      fgd2YGood = false;
      }
    if (fgd1XGood && fabs(fgd1TrackAngleX) > lowAngleCut && fabs(fgd1TrackAngleX) < highAngleCut)
      for (int i = 0; i < 30; i+=2)
        {
    if (res[i]!=9999 && layerHitCnt[i]==hitCntPerLayer)
      {
     
        //res[i]+=angleConst*tan(fgd1TrackAngleX)*(fgd1XInt + fgd1SepDist*i/2 - fgd1Center);
        if (fgd1TrackAngleX > 0)
          resHistPosAngle[i]->Fill(res[i]);
        else
          resHistNegAngle[i]->Fill(res[i]);
        resHist[i]->Fill(res[i]);
        
      }
        }
    if (fgd1YGood && fabs(fgd1TrackAngleY) > lowAngleCut && fabs(fgd1TrackAngleY) < highAngleCut)
      for (int i = 1; i < 30; i+=2)
        {
    if (res[i]!=9999 && layerHitCnt[i]==hitCntPerLayer)
      {
        
        //res[i]+=angleConst*tan(fgd1TrackAngleY)*(fgd1YInt + fgd1SepDist*(i-1)/2 - fgd1Center);
        if (fgd1TrackAngleY > 0)
          resHistPosAngle[i]->Fill(res[i]);
        else
          resHistNegAngle[i]->Fill(res[i]);
        resHist[i]->Fill(res[i]);
        
      }
        }
    if (fgd2XGood && fabs(fgd2TrackAngleX) > lowAngleCut && fabs(fgd2TrackAngleX) < highAngleCut)
      for (int i = 30; i < 44; i+=2)
        {
    if (res[i]!=9999 && layerHitCnt[i]==hitCntPerLayer)
      {
        
        //res[i]+=angleConst*tan(fgd2TrackAngleX)*(fgd2XInt + fgd2SepDist*(i-30)/2 - fgd2Center);
        if (fgd2TrackAngleX > 0)
          resHistPosAngle[i]->Fill(res[i]);
        else
          resHistNegAngle[i]->Fill(res[i]);
        resHist[i]->Fill(res[i]);
        
      }
        }
    if (fgd2YGood && fabs(fgd2TrackAngleY) > lowAngleCut && fabs(fgd2TrackAngleY) < highAngleCut)
      for (int i = 31; i < 44; i+=2)
        {
    if (res[i]!=9999 && layerHitCnt[i]==hitCntPerLayer)
      {
        
        //res[i]+=angleConst*tan(fgd2TrackAngleY)*(fgd2YInt + fgd2SepDist*(i-31)/2 - fgd2Center);
        if (fgd2TrackAngleY > 0)
          resHistPosAngle[i]->Fill(res[i]);
        else
          resHistNegAngle[i]->Fill(res[i]);
        resHist[i]->Fill(res[i]);
        
      }
        } 
  }
    }
  
  TF1 *f1 = new TF1("f1", "gaus", -fitRangeHalf, fitRangeHalf);
  TF1 *fPol1FGD1 = new TF1("f2", "pol1", 0, 29);
  TF1 *fPol1FGD2 = new TF1("f3", "pol1", 30, 43);
  resCanvasFGD1X->Divide(5, 3);
  resCanvasFGD1Y->Divide(5, 3);
  resCanvasFGD2X->Divide(4, 2);
  resCanvasFGD2X->cd(8)->Close();
  resCanvasFGD2Y->Divide(4, 2);
  resCanvasFGD2Y->cd(8)->Close();
  for (int i = 0; i < 44; i++)
    {
      int cat = FGDCategory(i);
      HitCnt0[i] = log(HitCnt0[i]);
      HitCnt1[i] = log(HitCnt1[i]);
      HitCnt2[i] = log(HitCnt2[i]);
      HitCnt3[i] = log(HitCnt3[i]);
      if (divideAngle)
  {
    if (cat == 0)
      resCanvasFGD1X->cd(i/2+1);
    else if (cat == 1)
      resCanvasFGD1Y->cd((i-1)/2+1);
    else if (cat == 2)
      resCanvasFGD2X->cd((i-30)/2+1);
    else 
      resCanvasFGD2Y->cd((i-31)/2+1);
    resHistNegAngle[i]->GetXaxis()->SetTitle("residual (mm)");
    resHistNegAngle[i]->SetLineColor(1);
    resHistNegAngle[i]->Draw();
    resHistPosAngle[i]->SetLineColor(2);
    resHistPosAngle[i]->Draw("sames");	
  }
      else
  {
    if (cat == 0)
      resCanvasFGD1X->cd(i/2+1);
    else if (cat == 1)
      resCanvasFGD1Y->cd((i-1)/2+1);
    else if (cat == 2)
      resCanvasFGD2X->cd((i-30)/2+1);
    else 
      resCanvasFGD2Y->cd((i-31)/2+1);
    resHist[i]->GetXaxis()->SetTitle("residual (mm)");
    resHist[i]->Draw();

  }
    
      //std::cout<<cat<<std::endl;
      resHist[i]->Fit("f1", "OQR");
      resHistNegAngle[i]->Fit("f1", "QR");
      resHistPosAngle[i]->Fit("f1", "QR");

      if (useGauss)
  {
    if (cat == 0)
      {
        meanResFGD1X[i/2] = resHist[i]->GetFunction("f1")->GetParameter(1);
        meanResErrFGD1X[i/2] = resHist[i]->GetFunction("f1")->GetParError(1);
      }
    else if (cat == 1)
      {
        meanResFGD1Y[(i-1)/2] = resHist[i]->GetFunction("f1")->GetParameter(1);
        meanResErrFGD1Y[(i-1)/2] = resHist[i]->GetFunction("f1")->GetParError(1);
      }
    else if (cat == 2)
      {
        meanResFGD2X[(i-30)/2] = resHist[i]->GetFunction("f1")->GetParameter(1);
        meanResErrFGD2X[(i-30)/2] = resHist[i]->GetFunction("f1")->GetParError(1);
      }
    else 
      {
        meanResFGD2Y[(i-31)/2] = resHist[i]->GetFunction("f1")->GetParameter(1);
        meanResErrFGD2Y[(i-31)/2] = resHist[i]->GetFunction("f1")->GetParError(1);
      }
    meanResDiff[i] = resHist[i]->GetFunction("f1")->GetParameter(1) - fgdIntTransConst[i];
    sigmaOrRMS[i] = resHist[i]->GetFunction("f1")->GetParameter(2);
    sigmaOrRMSPos[i] = resHistPosAngle[i]->GetFunction("f1")->GetParameter(2);
    sigmaOrRMSNeg[i] = resHistNegAngle[i]->GetFunction("f1")->GetParameter(2);
    meanResNegAngle[i] = resHistNegAngle[i]->GetFunction("f1")->GetParameter(1) - fgdIntTransConst[i];
    meanResNegAngleErr[i] = resHistNegAngle[i]->GetFunction("f1")->GetParError(1);
    meanResPosAngle[i] = resHistPosAngle[i]->GetFunction("f1")->GetParameter(1) - fgdIntTransConst[i];
    meanResPosAngleErr[i] = resHistPosAngle[i]->GetFunction("f1")->GetParError(1);
        }

      else
  { 
    if (cat == 0)
      {
        meanResFGD1X[i/2] = resHist[i]->GetMean(1);
        meanResErrFGD1X[i/2] = resHist[i]->GetMeanError(1);
      }
    else if (cat == 1)
      {
        meanResFGD1Y[(i-1)/2] = resHist[i]->GetMean(1);
        meanResErrFGD1Y[(i-1)/2] = resHist[i]->GetMeanError(1);
      }
    else if (cat == 2)
      {
        meanResFGD2X[(i-30)/2] = resHist[i]->GetMean(1);
        meanResErrFGD2X[(i-30)/2] = resHist[i]->GetMeanError(1);
      }
    else 
      {
        meanResFGD2Y[(i-31)/2] = resHist[i]->GetMean(1);
        meanResErrFGD2Y[(i-31)/2] = resHist[i]->GetMeanError(1);
      }
    sigmaOrRMS[i] = resHist[i]->GetRMS();
    sigmaOrRMSPos[i] = resHistPosAngle[i]->GetRMS();
    sigmaOrRMSNeg[i] = resHistNegAngle[i]->GetRMS();
    meanResDiff[i] = resHist[i]->GetMean(1) - fgdIntTransConst[i];
    meanResNegAngle[i] = resHistNegAngle[i]->GetMean(1) - fgdIntTransConst[i];
    meanResNegAngleErr[i] = resHistNegAngle[i]->GetMeanError(1);
    meanResPosAngle[i] = resHistPosAngle[i]->GetMean(1) - fgdIntTransConst[i];
    meanResPosAngleErr[i] = resHistPosAngle[i]->GetMeanError(1);
  }

      resDiffNegRMSHist->Fill(meanResNegAngle[i]);
      resDiffPosRMSHist->Fill(meanResPosAngle[i]);
      resDiffRMSHist->Fill(meanResDiff[i]);
               
      if (cat == 0)
  chisqFGD1X += pow((meanResFGD1X[i/2]-fgdIntTransConst[i])/meanResErrFGD1X[i/2], 2);
      else if (cat == 1){
  //if (i != 21)
  chisqFGD1Y += pow((meanResFGD1Y[(i-1)/2]-fgdIntTransConst[i])/meanResErrFGD1Y[(i-1)/2], 2);
      }
      else if (cat == 2){
  //if (i != 38 && i != 40)
  chisqFGD2X += pow((meanResFGD2X[(i-30)/2]-fgdIntTransConst[i])/meanResErrFGD2X[(i-30)/2], 2);
      }
      else{
  //if (i != 43)
  chisqFGD2Y += pow((meanResFGD2Y[(i-31)/2]-fgdIntTransConst[i])/meanResErrFGD2Y[(i-31)/2], 2);
  }
    }
  
  resGraphFGD1X = new TGraphErrors(15, layerFGD1X, meanResFGD1X, layerErrFGD1, meanResErrFGD1X);
  resGraphFGD1X->SetMarkerStyle(22);
  resGraphFGD1X->SetMarkerColor(2);

  resGraphFGD1Y = new TGraphErrors(15, layerFGD1Y, meanResFGD1Y, layerErrFGD1, meanResErrFGD1Y);
  resGraphFGD1Y->SetMarkerStyle(22);
  resGraphFGD1Y->SetMarkerColor(4);

  resGraphFGD2X = new TGraphErrors(7, layerFGD2X, meanResFGD2X, layerErrFGD2, meanResErrFGD2X);
  resGraphFGD2X->SetMarkerStyle(23);
  resGraphFGD2X->SetMarkerColor(2);

  resGraphFGD2Y = new TGraphErrors(7, layerFGD2Y, meanResFGD2Y, layerErrFGD2, meanResErrFGD2Y);
  resGraphFGD2Y->SetMarkerStyle(23);
  resGraphFGD2Y->SetMarkerColor(4);

  resDiffGraph = new TGraphErrors(44, layer, meanResDiff, layerErr, layerErr);
  resDiffGraph->SetMarkerStyle(20);
  resDiffGraph->SetMarkerColor(3);
  
  sigmaOrRMSGraph = new TGraphErrors(44, layer, sigmaOrRMS, layerErr, layerErr);
  sigmaOrRMSGraph->SetMarkerStyle(20);
  sigmaOrRMSGraph->SetMarkerColor(1);

  sigmaOrRMSPosGraph = new TGraphErrors(44, layer, sigmaOrRMSPos, layerErr, layerErr);
  sigmaOrRMSPosGraph->SetMarkerStyle(20);
  
  sigmaOrRMSNegGraph = new TGraphErrors(44, layer, sigmaOrRMSNeg, layerErr, layerErr);
  sigmaOrRMSNegGraph->SetMarkerStyle(20);

  resNegAngleGraph = new TGraphErrors(44, layer, meanResNegAngle, layerErr, meanResNegAngleErr);
  resNegAngleGraph->SetMarkerStyle(20);
  resNegAngleGraph->SetMarkerColor(2);

 

  resPosAngleGraph = new TGraphErrors(44, layer, meanResPosAngle, layerErr, meanResPosAngleErr);
  resPosAngleGraph->SetMarkerStyle(20);
  resPosAngleGraph->SetMarkerColor(4);
  

  origTrans = new TGraphErrors(44, layer, fgdIntTransConst, layerErr, layerErr);
  origTrans->SetMarkerStyle(29);
  origTrans->SetMarkerColor(1);

  refLine1 = new TGraphErrors(44, layer, layerErr, layerErr, layerErr);
  refLine1->GetXaxis()->SetTitle("layer number");
  refLine1->GetYaxis()->SetTitle("mean residual (mm)");
  //refLine->SetLineColor(1);
 

  double fgdDivideX[2];
  double fgdDivideY[2];
  fgdDivideX[0] = 29.5;
  fgdDivideX[1] = 29.5;
  fgdDivideY[0] = -1;
  fgdDivideY[1] = 1;
  fgdDivider = new TGraphErrors(2, fgdDivideX, fgdDivideY, layerErr, layerErr);
  
  graphCanvas = new TCanvas("residual graph", "residual graph", 1000, 700);

  legend = new TLegend(0.1,0.7,0.3,0.9);
  legend->AddEntry(refLine, "0 line, FGD divider", "l");
  legend->AddEntry(origTrans, "misalignment constants", "p");
  refLine1->SetMinimum(-1);
  refLine1->SetMaximum(1);
  refLine1->Draw("al");
  fgdDivider->Draw("l");
  origTrans->Draw("p");
 
  legend->AddEntry(resGraphFGD1X, "calculated mean residual of FGD1X", "p");
  legend->AddEntry(resGraphFGD1Y, "calculated mean residual of FGD1Y", "p");
  legend->AddEntry(resGraphFGD2X, "calculated mean residual of FGD2X", "p");
  legend->AddEntry(resGraphFGD2Y, "calculated mean residual of FGD2Y", "p");
  legend->AddEntry(resDiffGraph, "calculated mean residual - original constant", "p");
  resGraphFGD1X->Draw("p");
  resGraphFGD1Y->Draw("p");
  resGraphFGD2X->Draw("p");
  resGraphFGD2Y->Draw("p");
  resDiffGraph->Draw("p");
  legend->Draw(); 

  legend1 = new TLegend(0.1,0.7,0.3,0.9);
  refLine2 = new TGraphErrors(44, layer, layerErr, layerErr, layerErr);
  refLine2->GetXaxis()->SetTitle("layer number");
  refLine2->GetYaxis()->SetTitle("mean residual difference to original (mm)");
  refLine2->SetMinimum(-0.5);
  refLine2->SetMaximum(0.5);
  resDiffGraphCanvas = new TCanvas("difference in mean residuals with track angle", "difference in mean residuals with track angle", 1000, 700);
  refLine2->Draw("al");
  fgdDivider->Draw("l");
  legend1->AddEntry(resNegAngleGraph, "from tracks with angle < 0", "p");
  legend1->AddEntry(resPosAngleGraph, "from tracks with angle > 0", "p");
  resNegAngleGraph->Draw("p");
  resPosAngleGraph->Draw("p");
  resNegAngleGraph->Fit("f2", "OQR");
  resNegAngleGraph->GetFunction("f2")->DrawCopy("LSAMES");
  resNegAngleGraph->Fit("f3", "OQR");
  resNegAngleGraph->GetFunction("f3")->DrawCopy("LSAMES");
  resPosAngleGraph->Fit("f2", "OQR");
  resPosAngleGraph->GetFunction("f2")->DrawCopy("LSAMES");
  resPosAngleGraph->Fit("f3", "OQR");
  resPosAngleGraph->GetFunction("f3")->DrawCopy("LSAMES");
  legend1->Draw();

 
  
  std::cout<<"Chi2/ndof for FGD1X layers: "<<chisqFGD1X<<"/15"<<std::endl;
  std::cout<<"Chi2/ndof for FGD1Y layers: "<<chisqFGD1Y<<"/15"<<std::endl;
  std::cout<<"Chi2/ndof for FGD2X layers: "<<chisqFGD2X<<"/7"<<std::endl;
  std::cout<<"Chi2/ndof for FGD2Y layers: "<<chisqFGD2Y<<"/7"<<std::endl;
  std::cout<<"Total chi2/ndof: "<<chisqFGD1X+chisqFGD1Y+chisqFGD2X+chisqFGD2Y<<"/44"<<std::endl;
  std::cout<<"Probability: "<<TMath::Prob(chisqFGD1X+chisqFGD1Y+chisqFGD2X+chisqFGD2Y, 44)<<std::endl;



  double FGD1XTransConst[15];
  double FGD1YTransConst[15];
  double FGD2XTransConst[7];
  double FGD2YTransConst[7];	
  
  double FGD1XResDiff[15];
  double FGD1YResDiff[15];
  double FGD2XResDiff[7];
  double FGD2YResDiff[7];
  
  double FGD1XLayer[15];
  double FGD1YLayer[15];
  double FGD2XLayer[7];
  double FGD2YLayer[7];
    
 for (int i = 0; i < 15; i++)
   {

     FGD1XTransConst[i] = fgdIntTransConst[i*2];
     FGD1YTransConst[i] = fgdIntTransConst[i*2+1];

     FGD1XResDiff[i] = meanResFGD1X[i] - fgdIntTransConst[i*2];
     FGD1YResDiff[i] = meanResFGD1Y[i] - fgdIntTransConst[i*2+1];
     FGD1XLayer[i] = i;
     FGD1YLayer[i] = i;
   }
  for (int i = 0; i < 7; i++)
   {

     FGD2XTransConst[i] = fgdIntTransConst[i*2+30];
     FGD2YTransConst[i] = fgdIntTransConst[i*2+31];

     FGD2XResDiff[i] = meanResFGD2X[i] - fgdIntTransConst[i*2+30];
     FGD2YResDiff[i] = meanResFGD2Y[i] - fgdIntTransConst[i*2+31];
     FGD2XLayer[i] = i;
     FGD2YLayer[i] = i;

   }
 
 
 TGraphErrors* corrPlotFGD1X = new TGraphErrors(15, FGD1XTransConst, meanResFGD1X, layerErr, meanResErrFGD1X);
 corrPlotFGD1X->GetXaxis()->SetTitle("FGD1X Original misalignment (mm)");
 corrPlotFGD1X->GetYaxis()->SetTitle("FGD1X Calculated misalignment (mm)");
 TGraphErrors* corrPlotFGD1Y = new TGraphErrors(15, FGD1YTransConst, meanResFGD1Y, layerErr, meanResErrFGD1Y);
 corrPlotFGD1Y->GetXaxis()->SetTitle("FGD1Y Original misalignment (mm)");
 corrPlotFGD1Y->GetYaxis()->SetTitle("FGD1Y Calculated misalignment (mm)");
 TGraphErrors* corrPlotFGD2X = new TGraphErrors(7, FGD2XTransConst, meanResFGD2X, layerErr, meanResErrFGD2X);
 corrPlotFGD2X->GetXaxis()->SetTitle("FGD2X Original misalignment (mm)");
 corrPlotFGD2X->GetYaxis()->SetTitle("FGD2X Calculated misalignment (mm)");
 TGraphErrors* corrPlotFGD2Y = new TGraphErrors(7, FGD2YTransConst, meanResFGD2Y, layerErr, meanResErrFGD2Y);
 corrPlotFGD2Y->GetXaxis()->SetTitle("FGD2Y Original misalignment (mm)");
 corrPlotFGD2Y->GetYaxis()->SetTitle("FGD2Y Calculated misalignment (mm)");
 TCanvas* corrPlotCanvas = new TCanvas("correlation canvas", "correlation canvas", 1000, 700);
 corrPlotCanvas->Divide(2, 2);
 corrPlotCanvas->cd(1);
 corrPlotFGD1X->Fit("pol1", "Q");
 corrPlotFGD1X->Draw("ap");
 corrPlotCanvas->cd(2);
 corrPlotFGD1Y->Fit("pol1", "Q");
 corrPlotFGD1Y->Draw("ap");
 corrPlotCanvas->cd(3);
 corrPlotFGD2X->Fit("pol1", "Q");
 corrPlotFGD2X->Draw("ap");
 corrPlotCanvas->cd(4);
 corrPlotFGD2Y->Fit("pol1", "Q");
 corrPlotFGD2Y->Draw("ap");
 
 TGraphErrors* geomPlotFGD1X = new TGraphErrors(15, FGD1XLayer, FGD1XResDiff, layerErr, meanResErrFGD1X);
 geomPlotFGD1X->SetMarkerStyle(4);
 geomPlotFGD1X->GetXaxis()->SetTitle("FGD1X Original misalignment (mm)");
 geomPlotFGD1X->GetYaxis()->SetTitle("FGD1X Calculated misalignment (mm)");
 TGraphErrors* geomPlotFGD1Y = new TGraphErrors(15, FGD1YLayer, FGD1YResDiff, layerErr, meanResErrFGD1Y);
 geomPlotFGD1Y->SetMarkerStyle(4);
 geomPlotFGD1Y->GetXaxis()->SetTitle("FGD1Y Original misalignment (mm)");
 geomPlotFGD1Y->GetYaxis()->SetTitle("FGD1Y Calculated misalignment (mm)");
 TGraphErrors* geomPlotFGD2X = new TGraphErrors(7, FGD2XLayer, FGD2XResDiff, layerErr, meanResErrFGD2X);
 geomPlotFGD2X->SetMarkerStyle(4);
 geomPlotFGD2X->GetXaxis()->SetTitle("FGD2X Original misalignment (mm)");
 geomPlotFGD2X->GetYaxis()->SetTitle("FGD2X Calculated misalignment (mm)");
 TGraphErrors* geomPlotFGD2Y = new TGraphErrors(7, FGD2YLayer, FGD2YResDiff, layerErr, meanResErrFGD2Y);
 geomPlotFGD2Y->SetMarkerStyle(4);
 geomPlotFGD2Y->GetXaxis()->SetTitle("FGD2Y Original misalignment (mm)");
 geomPlotFGD2Y->GetYaxis()->SetTitle("FGD2Y Calculated misalignment (mm)");
 TCanvas* geomPlotCanvas = new TCanvas("geometry correlation canvas", "geometry correlation canvas", 1000, 700);
 geomPlotCanvas->Divide(2, 2);
 geomPlotCanvas->cd(1);
 //geomPlotFGD1X->Fit("pol1", "Q");
 geomPlotFGD1X->Draw("ap");
 geomPlotCanvas->cd(2);
 //geomPlotFGD1Y->Fit("pol1", "Q");
 geomPlotFGD1Y->Draw("ap");
 geomPlotCanvas->cd(3);
 //geomPlotFGD2X->Fit("pol1", "Q");
 geomPlotFGD2X->Draw("ap");
 geomPlotCanvas->cd(4);
 //geomPlotFGD2Y->Fit("pol1", "Q");
 geomPlotFGD2Y->Draw("ap");

 TCanvas* angCanvas = new TCanvas("angular distribution", "angular distribution", 1000, 700);
 angCanvas->Divide(2, 2);
 angCanvas->cd(1);
 angHistFGD1X->Draw();
 angCanvas->cd(2);
 angHistFGD1Y->Draw();
 angCanvas->cd(3);
 angHistFGD2X->Draw();
 angCanvas->cd(4);
 angHistFGD2Y->Draw();

 TCanvas* sigmaOrRMSCanvas = new TCanvas("sigma or RMS distribution", "sigma or RMS distribution", 1000, 700);
 sigmaOrRMSCanvas->cd(1);
 if (divideAngle)
   {
     sigmaOrRMSPosGraph->SetMarkerColor(4);
     sigmaOrRMSPosGraph->Draw("ap");
     sigmaOrRMSNegGraph->SetMarkerColor(2);
     sigmaOrRMSNegGraph->Draw("p");
   }
 else
 sigmaOrRMSGraph->Draw("ap");

 TCanvas* resDiffRMSCanvas = new TCanvas("mean residual difference distribution", "mean residual difference distribution", 1000, 700);
 resDiffRMSCanvas->cd(1);
 if (divideAngle)
   {
     resDiffRMSHist->GetXaxis()->SetTitle("calculated mean residual - misalignment constant (mm)");
     resDiffRMSHist->Draw();
     //resDiffNegRMSHist->GetXaxis()->SetTitle("calculated mean residual - misalignment constant (mm)");
     resDiffNegRMSHist->SetLineColor(2);
     resDiffNegRMSHist->Draw("sames");
     resDiffPosRMSHist->SetLineColor(4);
     resDiffPosRMSHist->Draw("sames");
   }
 else
   {
     resDiffRMSHist->GetXaxis()->SetTitle("calculated mean residual - misalignment constant (mm)");
     resDiffRMSHist->Draw();
   }
 for (int i = 0; i < 44; i+=2)
   {
     std::cout<<"Module "<<i/2<<" mean residuals: "<<-meanResDiff[i] - fgdIntTransConst[i]<<", "<<-meanResDiff[i+1] - fgdIntTransConst[i+1]<<" mm"<<std::endl;
   }

 HitCnt0Graph = new TGraphErrors(44, layer, HitCnt0, layerErr, layerErr);
 HitCnt0Graph->SetMarkerStyle(20);
 HitCnt1Graph = new TGraphErrors(44, layer, HitCnt1, layerErr, layerErr);
 HitCnt1Graph->SetMarkerStyle(20);
 HitCnt1Graph->SetMarkerColor(2);
 HitCnt1Graph->SetMinimum(0);
 HitCnt1Graph->SetMaximum(12);
 HitCnt2Graph = new TGraphErrors(44, layer, HitCnt2, layerErr, layerErr);
 HitCnt2Graph->SetMarkerStyle(20);
 HitCnt2Graph->SetMarkerColor(3);
 HitCnt3Graph = new TGraphErrors(44, layer, HitCnt3, layerErr, layerErr);
 HitCnt3Graph->SetMarkerStyle(20);
 HitCnt3Graph->SetMarkerColor(4);
 TCanvas* hitCntCanvas = new TCanvas("hit count distribution", "hit count distribution", 1000, 700);
 hitCntCanvas->cd(1);
 //HitCnt0Graph->Draw("alp");
 HitCnt1Graph->Draw("alp");
 HitCnt2Graph->Draw("lp");
 HitCnt3Graph->Draw("lp");

} 
int FGDCategory(int layer)
{
  if (layer >= 0 && layer < 44)
    {
      if (layer < 30 && layer%2 == 0)
  return 0;
      else if (layer < 30)
  return 1;
      else if (layer%2 == 0)
  return 2;
      else
  return 3;
    }
}