#include "Data_Dp.C"
#include "Data_D0.C"

const static double beam_angle = 0.011;
const static double Ecms = 4.180;
const static double M_Dp = 1.86962;
const static double M_D0 = 1.86483;
const static double M_Dstar = 2.01026;

TLorentzVector ecms(Ecms*sin(beam_angle), 0, 0, Ecms);

void make_CrBg(TF1* fun_in, TF1* fun_out1, TF1* fun_out2){
  for( int i = 0; i < 3; i++ ){
    fun_out1->SetParameter(i, fun_in->GetParameter(i));
    fun_out1->SetParError(i, fun_in->GetParError(i));
    }
  for(int i = 3; i < fun_in->GetNpar(); i++){
    fun_out2->SetParameter(i-3, fun_in->GetParameter(i));
    fun_out2->SetParError(i-3, fun_in->GetParError(i));
    }
}

void setpar(TF1* fun, vector<double> list_of_values, vector<string> list_of_names){
  for ( int i = 0; i < list_of_values.size(); i++){
    fun->SetParameter(i, list_of_values.at(i));
    fun->SetParName(i, list_of_names.at(i).c_str());
    }
}

void make_grid(TH1* his){
  gStyle->SetOptStat(10);
  gStyle->SetOptFit(1);
  gPad->SetGrid();
  his->SetMinimum(0);
}

void make_legend(TH1* his, vector<TF1*> list_of_function, vector<string> list_of_title){
  TLegend* legend = new TLegend(0.1, 0.7, 0.4, 0.9);
  legend->SetHeader("Data","C");
  legend->AddEntry(his, his->GetName());
  for ( int i = 0; i < list_of_title.size(); i++ )
    legend->AddEntry(list_of_function.at(i), list_of_title.at(i).c_str());
  legend->Draw();
}

void DrwFun(vector<TF1*> list_of_function){
  list_of_function.at(0)->SetLineColor(kRed);
  list_of_function.at(1)->SetLineColor(kGreen+2);
  list_of_function.at(2)->SetLineColor(kGray+2);
  
  list_of_function.at(0)->Draw("same");
  list_of_function.at(1)->Draw("same");
  list_of_function.at(2)->Draw("same");
}

void ShowInt(TF1* core, TF1* bg, TMatrixD& mat_core, TMatrixD& mat_bg)
{
  cout << "В распределнии Гаусса: " << core->Integral(core->GetParameter(1) - 3 * core->GetParameter(2), core->GetParameter(1) + 3 * core->GetParameter(2))
       << "\nВ линейном распределении: " << bg->Integral(core->GetParameter(1) - 3 * core->GetParameter(2), core->GetParameter(1) + 3 * core->GetParameter(2)) << endl;
  
  cout << "Вычисление ошибок инетгралов:"
       << "\nОшибка интеграла Гаусса: " << core->IntegralError(core->GetParameter(1) - 3 * core->GetParameter(2), core->GetParameter(1) + 3 * core->GetParameter(2), core->GetParameters(), mat_core.GetMatrixArray() )
       << "\nОшибка линейной функции: " << bg->IntegralError(core->GetParameter(1) - 3 * core->GetParameter(2), core->GetParameter(1) + 3 * core->GetParameter(2), bg->GetParameters(), mat_bg.GetMatrixArray() ) << endl;
}

void ShowInt(TF1* core, TF1* bg, TMatrixD& mat_core, TMatrixD& mat_bg, double main, double width)
{
  cout << "В распределнии Гаусса: " << core->Integral(main - width, main + width)
       << "\nВ линейном распределении: " << bg->Integral(main - width, main + width) << endl;
  
  cout << "Вычисление ошибок инетгралов:"
       << "\nОшибка интеграла Гаусса: " << core->IntegralError(main - width, main + width, core->GetParameters(), mat_core.GetMatrixArray() )
       << "\nОшибка линейной функции: " << bg->IntegralError(main - width, main + width, bg->GetParameters(), mat_bg.GetMatrixArray() ) << endl;
}

void build_his(){ 

////////////////////////////////////////////////////////////
////// Создание и заполнение гистограмм для Dp
////////////////////////////////////////////////////////////

  Data_Dp* Dp = new Data_Dp;

  // histogams for Dp
  TH1D* mass_D          = new TH1D("m_{D}","Mass D^{#pm};m_{D^{#pm}}, GeV", 80, 1.86966 - 0.040, 1.86966 + 0.040);
  TH1D* mass_BC         = new TH1D("m_{BC}","Mass BC;m_{BC}, GeV", 160, 1.83, 1.99);
  TH1D* num_ntrk        = new TH1D("ntrk","Number other tracks;N_{other}", 18, 0, 18);
  TH1D* num_ntrk_charge = new TH1D("ntrk_charge","Number other charge tracks;N_{other__charge}", 8, 0, 8);
  TH2D* mass_BC_D       = new TH2D("m_BC_D","Masses D and BC;m_{BC}, GeV;m_{D}, GeV", 160, 1.83, 1.99, 80, 1.86966 - 0.040, 1.86966 + 0.040);
  TH2D* mBC_eD          = new TH2D("m_BC_eD","Momentum D and masses BC;m_{BC}, GeV;p_{D}, GeV", 160, 1.83, 1.99, 400, 0.6, 1.);
  TH1D* M_BC            = new TH1D("Ь_ВС","Mass BC;m_{BC}, GeV", 330, 1.79, 2.12);
  TH1D* x_p             = new TH1D("x_P", "", 10, 0.65, 1.15);
  TH1D* num_DTag_Dp     = new TH1D("num_DTag_Dp","Number DTag for Dp", 10, 0, 10);
  
  TH1D* his_Dp[7] = {mass_D, mass_BC, num_ntrk, num_ntrk_charge, M_BC, x_p, num_DTag_Dp};
  
  TH1D* M_DD[9];
  for (int i = 0; i < 9; i++)
    M_DD[i] = new TH1D( ("M_D" + to_string(i)).c_str(), "Mass D from DTag;M_{D^{#pm}}, GeV", 70, 1.86966 - 0.016, 1.86966 + 0.016);
  
  //auto& mass_Dp = Dp->mass_false;
  int run_no_Dp = 0;
  int num_dtag_Dp = 0;
  
  for (int i = 0; i < Dp->fChain->GetEntriesFast(); i++){
    Dp->GetEntry(i);
    
    mass_D->Fill(Dp->mass_false);
    mass_BC->Fill(Dp->massBC_false);
    num_ntrk->Fill(Dp->ntrk);
    
    TLorentzVector mD;
    mD.SetXYZM(Dp->p_D * TMath::Cos(Dp->phi_D) * TMath::Sin( TMath::ACos(Dp->cosTheta_D) ), Dp->p_D * TMath::Sin(Dp->phi_D) * TMath::Sin( TMath::ACos(Dp->cosTheta_D) ), Dp->p_D * Dp->cosTheta_D, 1.86965);
    M_BC->Fill((mD-ecms).M());
    
    int k = 0;
    for ( int j = 0; j < Dp->ntrk; j++ )
      if ( Dp->charge_other[j] != 10 )
        k++;
    num_ntrk_charge->Fill(k);
    
    mass_BC_D->Fill(Dp->massBC_false, Dp->mass_false);
    mBC_eD->Fill(Dp->massBC_false, Dp->p_D);
    
    //if ( (mD-ecms).M() > 1.89 )
    x_p->Fill( Dp->p_D / sqrt( Ecms * Ecms/4 - M_Dp * M_Dp) );
    for ( int i = 0; i < 9; i++)
       if ( Dp->p_D / sqrt( Ecms * Ecms/4 - M_Dp * M_Dp) > 0.65 + i * 0.05 && Dp->p_D/sqrt( Ecms * Ecms/4 - M_Dp * M_Dp) < 0.70 + i * 0.05)
         M_DD[i]->Fill(Dp->mass_false);
    //if ( p_D/sqrt(Ecms*Ecms/4 - M_Dp * M_Dp) > 0.75 && p_D/sqrt(Ecms*Ecms/4 - M_Dp * M_Dp) < 0.8 )
    //M_D1->Fill(mass);
    
    if ( Dp->run_no == run_no_Dp )
      num_dtag_Dp++;
    else
      {
      run_no_Dp = Dp->run_no;
      num_DTag_Dp->Fill(num_dtag_Dp);
      num_dtag_Dp = 1;
      }
  }
  
////////////////////////////////////////////////////////////
////// Создание и заполнение гистограмм для D0
////////////////////////////////////////////////////////////

  Data_D0* D0 = new Data_D0;
  
  // hisograms for D0
  TH1D* mass_D0     = new TH1D("m_D0","Mass D^{0};m_{D^{0}}, GeV", 80, 1.865 - 0.040, 1.865 + 0.040);
  TH1D* mass_BC0    = new TH1D("m_BC0","Mass BC;m_{BC}, GeV", 160, 1.83, 1.99);
  TH1D* num_ntrk0   = new TH1D("ntrk0","Number other tracks;N_{other}", 18, 0, 18);
  TH1D* num_DTag_D0 = new TH1D("num_DTag_D0","Number DTag for D0", 10, 0, 10);
  
  TH1D* his_D0[4] = {mass_D0, mass_BC0, num_ntrk0, num_DTag_D0};
  
  int run_no_D0 = 0;
  int num_dtag_D0 = 0;
  
  for ( int i = 0; i < D0->fChain->GetEntriesFast(); i++ ){
    D0->GetEntry(i);
    //cout << D0->run_no0 << ": " << D0->mass_D0 << endl;
    mass_D0->Fill(D0->mass_D0);
    mass_BC0->Fill(D0->massBC_D0);
    num_ntrk0->Fill(D0->ntrk0);
    
    if ( D0->run_no0 == run_no_D0 )
      num_dtag_D0++;
    else
      {
      run_no_D0 = D0->run_no0;
      num_DTag_D0->Fill(num_dtag_D0);
      num_dtag_D0 = 1;
      }
  }
  
  //mass_D->Sumw2();
  //mass_D->Scale(1., "width");
  
  //mass_BC->Sumw2();
  //mass_BC->Scale(1., "width");

////////////////////////////////////////////////////////////
////// Фитирование и подсчет числа частиц
////////////////////////////////////////////////////////////

  TF1* f_total = new TF1("f_total","gausn(0) + pol1(3)", 1.86966 - 0.050, 1.86966 + 0.050);
  
  vector<double> val = {10., 1.86966, 0.004, 120., 0.};
  vector<string> na = {"N_{#infty}", "#mu", "#sigma", "p_{0}", "p_{1}"};
  
  setpar(f_total, val, na);
  
  TFitResultPtr res_fit_mD = mass_D->Fit(f_total, "RSNQ");
  
  TMatrixDSym covTot = res_fit_mD->GetCovarianceMatrix();
  
  TMatrixD covGauss = covTot.GetSub(0, 2, 0, 2);
  TMatrixD covLine = covTot.GetSub(3, 4, 3, 4);
  
  TF1* f_core = new TF1 ("f_core","gausn(0)", 1.86966 - 0.050, 1.86966 + 0.050);
  TF1* f_bg = new TF1 ("f_bg","pol1(0)", 1.86966 - 0.050, 1.86966 + 0.050);
  
  make_CrBg(f_total, f_core, f_bg);
  
  // Нахождение количества частиц в Гауссе и линейном распределении для Dp
  //ShowInt(f_core, f_bg, covGauss, covLine, M_Dp, 0.016);
       
  TF1* f_total_BC = new TF1("f_total_BC","gausn(0) + pol1(3)", 1.86966 - 0.016, 1.86966 + 0.016);
  
  vector<double> val1 = {10., 1.86966, 0.004, 120., 0.};
  vector<string> na1 = {"N_{#infty}", "#mu", "#sigma", "p_{0}", "p_{1}"};
  
  setpar(f_total_BC, val1, na1);
  
  TFitResultPtr res_fit_mBC = mass_BC->Fit(f_total_BC, "RSNQ");
  
  TMatrixDSym covTot_BC = res_fit_mBC->GetCovarianceMatrix();
  
  TMatrixD covGauss_BC = covTot_BC.GetSub(0, 2, 0, 2);
  TMatrixD covLine_BC = covTot_BC.GetSub(3, 4, 3, 4);
  
  TF1* f_core_BC = new TF1 ("f_core_BC","gausn(0)", 1.86966 - 0.016, 1.86966 + 0.016);
  
  TF1* f_bg_BC = new TF1 ("f_bg_BC","pol1(0)", 1.86966 - 0.016, 1.86966 + 0.016);
  
  make_CrBg(f_total_BC, f_core_BC, f_bg_BC);
  
  TF1* f_total_BC_2 = new TF1("f_total_BC_2","gausn(0) + pol1(3)", 2.01 - 0.016, 2.01 + 0.016);
  
  vector<double> val2 = {10., 2.01, 0.004, 120., 0.};
  vector<string> na2 = {"N_{#infty}", "#mu", "#sigma", "p_{0}", "p_{1}"};
  
  setpar(f_total_BC_2, val2, na2);
  
  TFitResultPtr res_fit_mBC_2 = M_BC->Fit(f_total_BC_2, "RSNQ");
  
  TMatrixDSym covTot_BC_2 = res_fit_mBC_2->GetCovarianceMatrix();
  
  TMatrixD covGauss_BC_2 = covTot_BC_2.GetSub(0, 2, 0, 2);
  TMatrixD covLine_BC_2 = covTot_BC_2.GetSub(3, 4, 3, 4);
  
  TF1* f_core_BC_2 = new TF1 ("f_core_BC_2","gausn(0)", 2.01 - 0.016, 2.01 + 0.016);
  
  TF1* f_bg_BC_2 = new TF1 ("f_bg_BC_2","pol1(0)", 2.01 - 0.016, 2.01 + 0.016);
  
  make_CrBg(f_total_BC_2, f_core_BC_2, f_bg_BC_2);
  
  // Нахождение количества частиц в Гауссе и линейном распределении для Dm из BC для Dp
  //ShowInt(f_core_BC, f_bg_BC, covGauss_BC, covLine_BC, M_Dp, 0.016);
  
  // Нахождение количества частиц в Гауссе и линейном распределении для второго пика Dstar из BC для Dp
  //ShowInt(f_core_BC_2, f_bg_BC_2, covGauss_BC_2, covLine_BC_2);
  
  //cout << "Отношение интеграла D мезона ко второму пику:" << f_core_BC->Integral(1.86966 - 0.016, 1.86966 + 0.016)/f_core_BC_2->Integral(2.012 - 0.016, 2.012 + 0.016) << endl;
  
  vector<TF1*> fun_Dp = {f_total, f_core, f_bg};
  vector<TF1*> fun_BC_Dp_Dm = {f_total_BC, f_core_BC, f_bg_BC};
  vector<TF1*> fun_BC_Dp_Dstar = {f_total_BC_2, f_core_BC_2, f_bg_BC_2};
  
  vector<string> names_Dp = {"Gauss + Bg", "Gauss", "Bg"};
  vector<string> names_BC_Dp_Dm = {"Gauss + Bg", "Gauss", "Bg"};
  vector<string> names_BC_Dp_Dstar = {"Gauss + Bg", "Gauss", "Bg"};
  
  vector<TF1*> list_fun[7] = {fun_Dp, fun_BC_Dp_Dm, {}, {}, {}, {}, {}};
  vector<string> list_names[7] = {names_Dp, names_BC_Dp_Dm, {}, {}, {}, {}, {}};
  
  TF1* f_total_D0 = new TF1("f_total_D0","gausn(0) + pol1(3)", M_D0 - 0.050, M_D0 + 0.050);
  
  vector<double> val_D0 = {10., M_D0, 0.004, 120., 0.};
  vector<string> na_val_D0 = {"N_{#infty}", "#mu", "#sigma", "p_{0}", "p_{1}"};
  
  setpar(f_total_D0, val_D0, na_val_D0);
  
  TFitResultPtr res_fit_D0 = mass_D0->Fit(f_total_D0, "RSNQ");
  
  TMatrixDSym covTot_D0 = res_fit_D0->GetCovarianceMatrix();
  
  TMatrixD covGauss_D0 = covTot_D0.GetSub(0, 2, 0, 2);
  TMatrixD covLine_D0 = covTot_D0.GetSub(3, 4, 3, 4);
  
  TF1* f_core_D0 = new TF1 ("f_core_D0","gausn(0)", M_D0 - 0.050, M_D0 + 0.050);
  
  TF1* f_bg_D0 = new TF1 ("f_bg_D0","pol1(0)", M_D0 - 0.050, M_D0 + 0.050);
  
  make_CrBg(f_total_D0, f_core_D0, f_bg_D0);
  
  vector<TF1*> fun_D0 = {f_total_D0, f_core_D0, f_bg_D0};
  vector<string> names_D0 = {"Gauss + Bg", "Gauss", "Bg"};
  
  vector<TF1*> list_fun_D0[4] = {fun_D0, {}, {}, {}};
  vector<string> list_names_D0[4] = {names_D0, {}, {}, {}};

////////////////////////////////////////////////////////////
////// Печатанье гистограмм
////////////////////////////////////////////////////////////
  
  string str_windows_Dp = "1000000"; // всего 7 цифр
  TCanvas* c_Dp[str_windows_Dp.size()];
  
  for ( int i = 0; i < str_windows_Dp.size(); i++ )
    if ( str_windows_Dp[i] == '1' )
    {
      c_Dp[i] = new TCanvas( ("Window " + to_string(i) + " for Dp").c_str(), "", 20 + i * 40, 20, 600, 600 );
      c_Dp[i]->cd();
      his_Dp[i]->Draw();
      make_grid(his_Dp[i]);
      if ( list_fun[i].empty() != 1 )
      {
        DrwFun(list_fun[i]);
        make_legend(his_Dp[i], list_fun[i], list_names[i]);
      }
    }
  
  string str_windows_D0 = "0000"; // всего 4 цифры
  TCanvas* c_D0[str_windows_D0.size()];
  
  for( int i = 0; i < str_windows_D0.size(); i++)
    if ( str_windows_D0[i] == '1' )
    {
      c_D0[i] = new TCanvas( ("Window " + to_string(i) + " for D0").c_str(),"", 20, 20, 600, 600 );
      c_D0[i]->cd();
      his_D0[i]->Draw();
      make_grid(his_D0[i]);
      if ( list_fun_D0[i].empty() != 1 )
      {
        DrwFun(list_fun_D0[i]);
        make_legend(his_D0[i], list_fun_D0[i], list_names_D0[i]);
      }
    }
  
  string str_windows = "000000000"; // всего 9 цифр
  TCanvas* c[str_windows.size()];
  
  double main_value[str_windows.size()];
  double error_value[str_windows.size()];
  double error_x[str_windows.size()];
  for ( int i = 0; i < 9; i++)
    error_x[i]=0;
  double x_value[str_windows.size()];
  for ( int i = 0; i < 9; i++)
    x_value[i] = 0.65 + 0.005 + i * 0.05;
  
  for ( int i = 0; i < str_windows.size(); i++)
    if ( str_windows[i] == '1' ){
      c[i] = new TCanvas( ("Windiw" + to_string(i)).c_str(), "", 20, 20, 600, 600);
      c[i]->cd();
      TF1* fTot = new TF1("fTot", "gausn(0) + pol2(3)", 1.86966 - 0.016, 1.86966 + 0.016);
      TF1* fGausn = new TF1("fGausn", "gausn(0)", 1.86966 - 0.016, 1.86966 + 0.016);
      TF1* fBg = new TF1("fBg", "pol2(0)", 1.86966 - 0.016, 1.86966 + 0.016);
      vector<double> v = {24., 1.86966, 0.004, 120., 0., 400.};
      vector<string> n = {"N_{#infty}", "#mu", "#sigma", "p_{0}", "p_{1}", "p_{2}"};
      setpar(fTot, v, n);
      TFitResultPtr res_fit_mD = M_DD[i]->Fit(fTot, "RSN");
      make_CrBg(fTot, fGausn, fBg);
      make_grid(M_DD[i]);
      vector<TF1*> funs = {fTot, fGausn, fBg};
      M_DD[i]->Draw("E");
      DrwFun(funs);
      fBg->Draw("same");
      main_value[i] = fTot->GetParameter(0);
      error_value[i] = fTot->GetParError(0);
      c[i]->Update();
      }
  //cout << error_value[7] << endl;
  TGraphErrors* xp_D = new TGraphErrors(9, x_value, main_value, error_x, error_value);
  
}