#include <fstream>
#include <iostream>
#include <TString.h>
#include <TH1D.h>
#include <TH1F.h>
#include <TF1.h>
#include <TFile.h>
#include <TLatex.h>
#include <TStyle.h>
#include <TProfile.h>

#include <stdio.h>
#include <cmath>

 #define Debug
 #define ERROR_OUT

using namespace std;

const int NStations = 2; // first and last

const int MaxNtr = 20000;

const char* fres_name = "../Output/track_param.dat";

float chi2[MaxNtr];

float Par[NStations][MaxNtr][6][3];		// [0] - value;  [1] - error; [2] - MC
TString Par_name[5]={"X","Y","Tx","Ty","qp"};
float Par_max[5][2];  // [0] - max value;  [1] - max (value - MCvalue)

float xboundary = 60;
float yboundary = 400;
float txboundary = 2.;
float tyboundary = 4;//3.;
float pboundary = 10;//4.;
const int nbins = 100;

void AsciiToRoot(TString nSt=""){

  int Ntr = 0; // number of tracks
  int NgoodTr = 0;     // number of tracks with normal parameters
#ifdef Debug 
  cout << nSt<<endl;
#endif
  if(!nSt.IsNull()) 
  {
    TString nSt_temp = nSt;
    nSt = "_";
    nSt+=nSt_temp;
  }

    //        fres_name += nSt.Data();
#ifdef Debug 
  cout << fres_name <<"111111111" << endl;
#endif
    // 	Style2m();
    // --------------------- Read results -----------------------------
  { 
    FILE * fpres; 
    fpres = fopen(fres_name, "r");

    int n = -1; 
#ifdef Debug 
    cout << "sds !!!!!!"<< feof(fpres)<<endl;
#endif
    for (int i = 0; (i < MaxNtr) && !feof(fpres); i++, NgoodTr++){
#ifdef Debug 
      // cout << i << endl;
#endif
      float Tmc[NStations][7],T[NStations][6],C[NStations][15];
      fscanf (fpres, "%d", &(n)); // nHits
#ifdef Debug 
    // cout << "NHits "<< n <<endl;
#endif
//       if (n!=i+Ntr-NgoodTr) {
// // #ifdef ERROR_OUT
//         cout << "Error read file. iTrack is wrong." << " n= " << n << " i+Ntr-NgoodTr= " << i+Ntr-NgoodTr << endl;
// // #endif // ERROR_OUT
//         return;
        //       }
      
      for(int is = 0; is<NStations; is++)
      {
        for(int k = 0; k<7; k++) fscanf (fpres, "%f", &(Tmc[0][k]));
        for(int k = 0; k<6; k++) fscanf (fpres, "%f", &(T[0][k]));
        for(int k = 0; k<15; k++) fscanf (fpres, "%f", &(C[0][k]));
          // there Tmc:  x,y,z,px,py,pz,q
        Par[is][i][0][2] = Tmc[0][0];
        Par[is][i][1][2] = Tmc[0][1];
        if(Tmc[0][5] == 0 ){std::cout << is << std::endl; continue;}
        Par[is][i][2][2] = Tmc[0][3]/Tmc[0][5];
        Par[is][i][3][2] = Tmc[0][4]/Tmc[0][5];
        Par[is][i][4][2] = Tmc[0][6]/sqrt(Tmc[0][3]*Tmc[0][3]+Tmc[0][4]*Tmc[0][4]+Tmc[0][5]*Tmc[0][5]);
        Par[is][i][5][2] = Tmc[0][2];

        for(int k = 0; k<5; k++){
          if (Par_max[k][0] < abs(T[0][k])) Par_max[k][0] = abs(T[0][k]);
          if (Par_max[k][1] < abs(Par[is][i][k][2]-T[0][k])) Par_max[k][0] = abs(Par[is][i][k][2]-T[0][k]);
        }

        for(int k = 0,kk=0; k<5; k++,kk+=k+1){
          Par[is][i][k][0] = T[0][k];
          if ( C[0][kk] > 0 ){
            Par[is][i][k][1] = sqrt(C[0][kk]);
          }
          else{
#ifdef ERROR_OUT
            cout << "Error: In track " << i+Ntr-NgoodTr << " er^2 < 0. Entry will be skiped." << endl;
#endif // ERROR_OUT
            Par[is][i][k][1] = -sqrt(-C[0][kk]); // mark as bad
            // i--;
            // NgoodTr--;
            // break;
          }
        }
        Par[is][i][5][0] = T[0][5];
        Par[is][i][5][1] = 1;//error

        // float tmp;
        // for(int k = 0; k<7+6+15; k++) fscanf (fpres, "%f", &tmp);
      }
    }
#ifdef Debug 
    cout << "sds !!!!!!"<< feof(fpres)<<endl;
#endif
    if (feof(fpres)) cout << " Data from file " << fres_name << " read succesfully " << endl; 
    fclose(fpres);
  } // read results
	

    // -------------- Fill histo ---------------------------

    //     TCanvas *c = new TCanvas("asdf","asdf",0,0,400,500);
    //     c ->Divide(2,3);
#ifdef Debug 
  cout << "tut" << endl;
  cout << fres_name <<"Fill histo"<< endl;
  cout << "N Read Tracks = " << NgoodTr << endl;
#endif
  TH1D *Pull[5],*X[5];//,*dX[5];

  TString FRootName;
  if(!nSt.IsNull()) 
  {
    TString name_temp = ".root";
    FRootName = "track_param";
    FRootName += nSt;
    FRootName += name_temp;
  }
  else FRootName = "track_param.root";
  TFile f(FRootName.Data(),"RECREATE");
#ifdef Debug 
  cout << FRootName << endl;
#endif

  
  for (int is = 0; is < NStations; is++) {
    TString Sta_name = "Station";
    Sta_name += is;
    Sta_name += ".";
    
    gStyle->SetOptStat(110010);
    for (int i = 0; i < NgoodTr; i++)
        // cout << Par[is][i][5][0] - Par[is][i][5][2] << endl;
      if ( (Par[is][i][5][0] - Par[is][i][5][2]) / Par[is][i][5][1] > 1e-2 ) {
        cout << " ERROR: Track: " << i << " DIFERENT Z OF RECO AND MC!! dz > 1e-2" << endl;
        cout << " Reco: " << Par[is][i][5][0] << " MC: " << Par[is][i][5][2] << endl;
        continue;
      }
  
    for (int k = 0; k < 5; k++){

        // ----- Pulls
      Pull[k] = new TH1D(Sta_name+"Pull "+Par_name[k],Sta_name+"Pull "+Par_name[k],nbins,-6.,6.);
      for (int i = 0; i < NgoodTr; i++){
        if( Par[is][i][k][1] > 0 ) Pull[k]->Fill( (Par[is][i][k][0] - Par[is][i][k][2]) / Par[is][i][k][1] );
      }

      TF1  *pullfit = new TF1(Sta_name+"gausFitPull "+Par_name[k],"gaus");
      pullfit->SetLineColor(2);
      Pull[k]->Fit(Sta_name+"gausFitPull "+Par_name[k],"","",Pull[k]->GetXaxis()->GetXmin(),Pull[k]->GetXaxis()->GetXmax());

      Pull[k]->GetXaxis()->SetTitle(Sta_name+"Pull "+Par_name[k]);
      Pull[k]->GetXaxis()->SetTitleSize(0.04);

        // ----- Values
      X[k] = new TH1D(Sta_name+Par_name[k],Sta_name+Par_name[k],nbins,-Par_max[k][0],Par_max[k][0]);
      for (int i = 0; i < NgoodTr; i++){
        // if ( Par[is][i][5][0] > 1. )
        X[k]->Fill( Par[is][i][k][0] );
      }

        // ----- Residual
        // 		dX[k] = new TH1D(Par_name[k]+" Residual",Par_name[k]+" Residual",nbins,-Par_max[k][1]/5,Par_max[k][1]/5);
        // 		for (int i = 0; i < NgoodTr; i++){
        // 			dX[k]->Fill( Par[is][i][k][0] - Par[is][i][k][2] );
        // 		}
        // 
        //   	TF1  *dXfit = new TF1("Residual "+ Par_name[k],"gaus");
        // 		dXfit->SetLineColor(2);
        //   	dX[k]->Fit("Residual "+Par_name[k],"","",dX[k]->GetXaxis()->GetXmin(),dX[k]->GetXaxis()->GetXmax());
        // 
        // 		dX[k]->GetXaxis()->SetTitle("Residual "+Par_name[k]);
        // 		dX[k]->GetXaxis()->SetTitleSize(0.04);

        //     c->cd(k);
        // 		Pull[k]->Draw("colz");
        // Save histo
      Pull[k]->Write();
      X[k]->Write();
        //  		dX[k]->Write();
    }

    TLatex Sss(1,2," Residual1");

      // ----- Residual x
    TH1D *dx = new TH1D(Sta_name+" Residualx",Sta_name+" Residualx",nbins,-xboundary,xboundary);
      // 	TH1D *dx = new TH1D(Sss,Sss,nbins,-200,200);
    for (int i = 0; i < NgoodTr; i++){
      dx->Fill( (Par[is][i][0][0] - Par[is][i][0][2])*10000. ); // from cm to um
    }

    TF1  *dxfit = new TF1(Sta_name+" Residualx","gaus");
    dxfit->SetLineColor(2);
    dx->Fit(Sta_name+" Residualx","","",dx->GetXaxis()->GetXmin(),dx->GetXaxis()->GetXmax());

    dx->GetXaxis()->SetTitle("Residual (x^{reco}-x^{mc}) [um]");
    dx->GetXaxis()->SetTitleSize(0.04);

      // ----- Residual tx
    TH1D *dtx = new TH1D(Sta_name+" Residualtx",Sta_name+" Residualtx",nbins,-txboundary,txboundary);
    for (int i = 0; i < NgoodTr; i++){
      dtx->Fill( (Par[is][i][2][0] - Par[is][i][2][2])*1000. );
    }

    TF1  *dtxfit = new TF1(Sta_name+" Residualtx","gaus");
    dtxfit->SetLineColor(2);
    dtx->Fit(Sta_name+" Residualtx","","",dtx->GetXaxis()->GetXmin(),dtx->GetXaxis()->GetXmax());

    dtx->GetXaxis()->SetTitle("Residual (t_{x}^{reco}-t_{x}^{mc}) [mrad]");
    dtx->GetXaxis()->SetTitleSize(0.04);

      // ----- Residual y
    TH1D *dy = new TH1D(Sta_name+" Residualy",Sta_name+" Residualy",nbins,-yboundary,yboundary);
      //  TH1D *dy = new TH1D(Sss,Sss,nbins,-200,200);
    for (int i = 0; i < NgoodTr; i++){
      dy->Fill( (Par[is][i][1][0] - Par[is][i][1][2])*10000. ); // from cm to um
    }

    TF1  *dyfit = new TF1(Sta_name+" Residualy","gaus");
    dyfit->SetLineColor(2);
    dy->Fit(Sta_name+" Residualy","","",dy->GetXaxis()->GetXmin(),dy->GetXaxis()->GetXmax());

    dy->GetXaxis()->SetTitle("Residual (y^{reco}-y^{mc}) [um]");
    dy->GetXaxis()->SetTitleSize(0.04);

      // ----- Residual ty
    TH1D *dty = new TH1D(Sta_name+" Residualty",Sta_name+" Residualty",nbins,-tyboundary,tyboundary);
    for (int i = 0; i < NgoodTr; i++){
      dty->Fill( (Par[is][i][3][0] - Par[is][i][3][2])*1000. );
    }

    TF1  *dtyfit = new TF1(Sta_name+" Residualty","gaus");
    dtyfit->SetLineColor(2);
    dty->Fit(Sta_name+" Residualty","","",dty->GetXaxis()->GetXmin(),dty->GetXaxis()->GetXmax());

    dty->GetXaxis()->SetTitle("Residual (t_{y}^{reco}-t_{y}^{mc}) [mrad]");
    dty->GetXaxis()->SetTitleSize(0.04);

  
      // ---- Resolution
    TH1D *resp = new TH1D(Sta_name+"p Resolution",Sta_name+"p Resolution",nbins,-pboundary,pboundary);
    for (int i = 0; i < NgoodTr; i++){
      if (Par[is][i][4][2] != 0.) resp->Fill( (fabs(Par[is][i][4][0])/fabs(Par[is][i][4][2]) - 1.0)*100 );
        // 		cout << ( Par[is][i][4][0]/Par[is][i][4][2] - 1.0 ) << endl;
    }


    TF1  *resfit = new TF1(Sta_name+"p Resolution","gaus");
    resfit->SetLineColor(2);
    resp->Fit(Sta_name+"p Resolution","","",resp->GetXaxis()->GetXmin(),resp->GetXaxis()->GetXmax());

    resp->GetXaxis()->SetTitle("Resolution (p^{reco}-p^{mc})/p^{mc}");
    resp->GetXaxis()->SetTitleSize(0.04);


      // ---- Resolution dependence of p
    double p[MaxNtr],rez[MaxNtr]; 
    for (int i = 0; i < NgoodTr; i++){
      if (Par[is][i][4][2] != 0.){
        rez[i] = fabs(Par[is][i][4][0] / Par[is][i][4][2] - 1.0) ;
        p[i] = fabs(Par[is][i][4][2]);
      }
    }
      //   TGraph *rezGr = new TGraph(NgoodTr,p,rez);
      //   rezGr->GetXaxis()->SetTitle("p\^\{mc\} [GeV/c]");
      //   rezGr->GetYaxis()->SetTitle("Resolution (p\^\{reco\}-p\^\{mc\})/p\^\{mc\}");
      //   rezGr->SetMarkerColor(2);
      //   rezGr->Draw();

    gStyle->SetOptStat(00000000);
    TProfile *rezGr = new TProfile(Sta_name+"Momentum resolution",Sta_name+"Momentum resolution",100, 0.0, 2.5,0,3.0);
    rezGr->SetBinContent(1,0);
    for (int i=0; i<NgoodTr; i++){
      rezGr->Fill(p[i],rez[i]*100);
        //     cout << p[i] << " " << rez[i] << endl;
    }
    rezGr->GetXaxis()->SetTitle("p^{mc} [GeV/c]");
    rezGr->GetYaxis()->SetTitle("Resolution (p^{reco}-p^{mc})/p^{mc},  %");
    rezGr->GetYaxis()->SetTitleOffset(1.3);
    rezGr->SetMarkerColor(2);
    rezGr->SetLineColor(2);
    rezGr->SetLineWidth(1);


    rezGr->Write();
    resp->Write();
    dtx->Write();
    dx->Write();
    dty->Write();
    dy->Write();

  }
}