#define Decay_cxx
#include "Decay.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

void Decay::Loop()
{
//   In a ROOT session, you can do:
//      root> .L Decay.C
//      root> Decay t
//      root> t.GetEntry(12); // Fill t data members with entry number 12
//      root> t.Show();       // Show values of entry 12
//      root> t.Show(16);     // Read and show values of entry 16
//      root> t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   TFile *histofile = new TFile("2011_MU_B_Mass_PID-BDTCuts_FrozenPion.root", "RECREATE");

   TH1D* h[50];
   TF1* f[50];

   gStyle->SetOptStat(1110);
   TH1D* h[1] = new TH1D("h_Reconstructed_B_Mass_cfg1", " ", 100, 6350, 6500);
   TH1D* h[2] = new TH1D("h_Reconstructed_B_Mass_BDTCut1_cfg1", " ", 100, 6350, 6500);
   TH1D* h[3] = new TH1D("h_Reconstructed_B_Mass_BDTCut2_cfg1", " ", 100, 6350, 6500);
   TH1D* h[4] = new TH1D("h_Reconstructed_B_Mass_BDTCut3_cfg1", " ", 100, 6350, 6500);
   TH1D* h[5] = new TH1D("h_Reconstructed_B_Mass_BDTCut4_cfg1", " ", 100, 6350, 6500);
   TH1D* h[6] = new TH1D("h_Reconstructed_B_Mass_BDTCut5_cfg1", " ", 100, 6350, 6500);
   TH1D* h[7] = new TH1D("h_Reconstructed_B_Mass_BDTCut6_cfg1", " ", 100, 6350, 6500);
   TH1D* h[8] = new TH1D("h_Reconstructed_B_Mass_BDTCut7_cfg1", " ", 100, 6350, 6500);
   TH1D* h[9] = new TH1D("h_Reconstructed_B_Mass_BDTCut8_cfg1", " ", 100, 6350, 6500);
   TH1D* h[10] = new TH1D("h_Reconstructed_B_Mass_BDTCut9_cfg1", " ", 100, 6350, 6500);
   TH1D* h[11]= new TH1D("h_Reconstructed_B_Mass_BDTCut10_cfg1", " ", 100, 6350, 6500);
   TH1D* h[12] = new TH1D("h_Reconstructed_B_Mass_BDTCut11_cfg1", " ", 100, 6350, 6500);

   TH1D* h[13] = new TH1D("h_Reconstructed_B_Mass_cfg2", " ", 100, 6350, 6500);
   TH1D* h[14]= new TH1D("h_Reconstructed_B_Mass_BDTCut1_cfg2", " ", 100, 6350, 6500);
   TH1D* h[15] = new TH1D("h_Reconstructed_B_Mass_BDTCut2_cfg2", " ", 100, 6350, 6500);
   TH1D* h[16] = new TH1D("h_Reconstructed_B_Mass_BDTCut3_cfg2", " ", 100, 6350, 6500);
   TH1D* h[17] = new TH1D("h_Reconstructed_B_Mass_BDTCut4_cfg2", " ", 100, 6350, 6500);
   TH1D* h[18] = new TH1D("h_Reconstructed_B_Mass_BDTCut5_cfg2", " ", 100, 6350, 6500);
   TH1D* h[19] = new TH1D("h_Reconstructed_B_Mass_BDTCut6_cfg2", " ", 100, 6350, 6500);
   TH1D* h[20] = new TH1D("h_Reconstructed_B_Mass_BDTCut7_cfg2", " ", 100, 6350, 6500);
   TH1D* h[21] = new TH1D("h_Reconstructed_B_Mass_BDTCut8_cfg2", " ", 100, 6350, 6500);
   TH1D* h[22] = new TH1D("h_Reconstructed_B_Mass_BDTCut9_cfg2", " ", 100, 6350, 6500);
   TH1D* h[23] = new TH1D("h_Reconstructed_B_Mass_BDTCut10_cfg2", " ", 100, 6350, 6500);
   TH1D* h[24] = new TH1D("h_Reconstructed_B_Mass_BDTCut11_cfg2", " ", 100, 6350, 6500);

   TH1D* h[25] = new TH1D("h_Reconstructed_B_Mass_cfg3", " ", 100, 6350, 6500);
   TH1D* h[26] = new TH1D("h_Reconstructed_B_Mass_BDTCut1_cfg3", " ", 100, 6350, 6500);
   TH1D* h[27] = new TH1D("h_Reconstructed_B_Mass_BDTCut2_cfg3", " ", 100, 6350, 6500);
   TH1D* h[28] = new TH1D("h_Reconstructed_B_Mass_BDTCut3_cfg3", " ", 100, 6350, 6500);
   TH1D* h[29] = new TH1D("h_Reconstructed_B_Mass_BDTCut4_cfg3", " ", 100, 6350, 6500);
   TH1D* h[30] = new TH1D("h_Reconstructed_B_Mass_BDTCut5_cfg3", " ", 100, 6350, 6500);
   TH1D* h[31] = new TH1D("h_Reconstructed_B_Mass_BDTCut6_cfg3", " ", 100, 6350, 6500);
   TH1D* h[32] = new TH1D("h_Reconstructed_B_Mass_BDTCut7_cfg3", " ", 100, 6350, 6500);
   TH1D* h[33] = new TH1D("h_Reconstructed_B_Mass_BDTCut8_cfg3", " ", 100, 6350, 6500);
   TH1D* h[34] = new TH1D("h_Reconstructed_B_Mass_BDTCut9_cfg3", " ", 100, 6350, 6500);
   TH1D* h[35] = new TH1D("h_Reconstructed_B_Mass_BDTCut10_cfg3", " ", 100, 6350, 6500);
   TH1D* h[36] = new TH1D("h_Reconstructed_B_Mass_BDTCut11_cfg3", " ", 100, 6350, 6500);

   TH1D* h[37] = new TH1D("h_Reconstructed_B_Mass_cfg4", " ", 100, 6350, 6500);
   TH1D* h[38] = new TH1D("h_Reconstructed_B_Mass_BDTCut1_cfg4", " ", 100, 6350, 6500);
   TH1D* h[39] = new TH1D("h_Reconstructed_B_Mass_BDTCut2_cfg4", " ", 100, 6350, 6500);
   TH1D* h[40] = new TH1D("h_Reconstructed_B_Mass_BDTCut3_cfg4", " ", 100, 6350, 6500);
   TH1D* h[41] = new TH1D("h_Reconstructed_B_Mass_BDTCut4_cfg4", " ", 100, 6350, 6500);
   TH1D* h[42] = new TH1D("h_Reconstructed_B_Mass_BDTCut5_cfg4", " ", 100, 6350, 6500);
   TH1D* h[43] = new TH1D("h_Reconstructed_B_Mass_BDTCut6_cfg4", " ", 100, 6350, 6500);
   TH1D* h[44] = new TH1D("h_Reconstructed_B_Mass_BDTCut7_cfg4", " ", 100, 6350, 6500);
   TH1D* h[45] = new TH1D("h_Reconstructed_B_Mass_BDTCut8_cfg4", " ", 100, 6350, 6500);
   TH1D* h[46] = new TH1D("h_Reconstructed_B_Mass_BDTCut9_cfg4", " ", 100, 6350, 6500);
   TH1D* h[47] = new TH1D("h_Reconstructed_B_Mass_BDTCut10_cfg4", " ", 100, 6350, 6500);
   TH1D* h[48] = new TH1D("h_Reconstructed_B_Mass_BDTCut11_cfg4", " ", 100, 6350, 6500);

   cout << "Histogramas declarados." << endl;

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;

     if ((d1_MC12TuneV2_ProbNNk > 0.2) && (d2_MC12TuneV2_ProbNNk > 0.2) && (d3_PIDK < -5)) {

       h[1]-> Fill (B_M);
       if (BDT_KKpi2011_total > 0) {
         h[2] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.02) {
         h[3] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.04) {
         h[4] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.06) {
         h[5] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.08) {
         h[6] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.1) {
         h[7] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.12) {
         h[8] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.14) {
         h[9] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.16) {
         h[10] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.18) {
         h[11] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.2) {
         h[12] -> Fill (B_M);
       }

     } // fim cfg1

     if ((d1_MC12TuneV2_ProbNNk > 0.3) && (d2_MC12TuneV2_ProbNNk > 0.3) && (d3_PIDK < -5)) {

       h[13] -> Fill (B_M);
       if (BDT_KKpi2011_total > 0) {
         h[14] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.02) {
         h[15] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.04) {
         h[16] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.06) {
         h[17] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.08) {
         h[18] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.1) {
         h[19] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.12) {
         h[20] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.14) {
         h[21] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.16) {
         h[22] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.18) {
         h[23] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.2) {
         h[24] -> Fill (B_M);
       }

     } // fim cfg2

     if ((d1_MC12TuneV2_ProbNNk > 0.4) && (d2_MC12TuneV2_ProbNNk > 0.4) && (d3_PIDK < -5)) {

       h[25] -> Fill (B_M);
       if (BDT_KKpi2011_total > 0) {
         h[26] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.02) {
         h[27] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.04) {
         h[28] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.06) {
         h[29] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.08) {
         h[30] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.1) {
         h[31] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.12) {
         h[32] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.14) {
         h[33] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.16) {
         h[34] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.18) {
         h[35] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.2) {
         h[36] -> Fill (B_M);
       }

     } // fim cfg3

     if ((d1_MC12TuneV2_ProbNNk > 0.5) && (d2_MC12TuneV2_ProbNNk > 0.5) && (d3_PIDK < -5)) {

       h[37] -> Fill (B_M);
       if (BDT_KKpi2011_total > 0) {
         h[38] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.02) {
         h[39] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.04) {
         h[40] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.06) {
         h[41] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.08) {
         h[42] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.1) {
         h[43] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.12) {
         h[44] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.14) {
         h[45] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.16) {
         h[46] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.18) {
         h[47] -> Fill (B_M);
       }
       if (BDT_KKpi2011_total > 0.2) {
         h[48] -> Fill (B_M);
       }

     } // fim cfg4
 
   } // fim do Long64

   cout << "Histogramas preenchidos." << endl;

   Int_t i;

   for (i=1; i<49; i++) {
     h[i]-> Sumw2();
   }

   cout << "Sumw2 calculado." << endl;

   histofile->Write();
   histofile->Close();

   Double_t UpEvents[100] = {0};
   Double_t B[100] = {0};
   Double_t SigEvents[100] = {0};
   Double_t integral = 0;

   cout << "Vetores declarados." << endl;

   for (i=1; i<49; i++) {
     integral = h[i]->Integral();
     UpEvents[i] = integral;
   }

   cout << "Vetor de número de eventos no upper sideband preenchido." << endl;

   // cfg1 

   for (i=1; i<7; i++) { // No cut to 0.08
     TF1* f[i] = new TF1("fit","pol1",6350,6490);
     h[i].Fit(f[i],"l","",6350,6490);
     B[i] = f[i]->Integral(6200,6350);
     SigEvents[i] = (double)B[i]/(double)1.5;
   }

   TF1* f[7] = new TF1("fit7","pol1",6350,6490); //0.10
   h[7].Fit(f[7],"l","",6350,6490);
   B[7] = f[7]->Integral(6200,6350);
   SigEvents[7] = B[7]/1.5;
   Double_t Cfg1Ratio = (double)SigEvents[7]/(double)UpEvents[7];

   for (i=8; i<10; i++) { //0.12 and 0.14
     TF1* f[i] = new TF1("fit","pol1",6350,6490);
     h[i].Fit(fi,"l","",6350,6490);
     B[i] = fi->Integral(6200,6350);
     SigEvents[i] = (double)B[i]/(double)1.5;
   }

   for (i=10; i<13; i++) { // 0.16 to 0.2
     SigEvents[i] = (double)UpEvents[i]*(double)Cfg1Ratio;
   }

   cout << "Variáveis referentes à primeira configuração calculadas." << endl;

   // cfg2

   for (i=13; i<19; i++) { // No cut to 0.08
     TF1* f[i] = new TF1("fit","pol1",6350,6490);
     h[i].Fit(f[i],"l","",6350,6490);
     B[i] = f[i]->Integral(6200,6350);
     SigEvents[i] = (double)B[i]/(double)1.5;
   }

   TF1* f[19] = new TF1("fit19","pol1",6350,6490); //0.10
   h[19].Fit(f[19],"l","",6350,6490);
   B[19] = f[19]->Integral(6200,6350);
   SigEvents[19] = B[19]/1.5;
   Double_t Cfg2Ratio = (double)SigEvents[19]/(double)UpEvents[19];

   for (i=20; i<22; i++) { //0.12 and 0.14
     TF1* f[i] = new TF1("fit","pol1",6350,6490);
     h[i].Fit(f[i],"l","",6350,6490);
     B[i] = f[i]->Integral(6200,6350);
     SigEvents[i] = (double)B[i]/(double)1.5;
   }

   for (i=22; i<25; i++) { // 0.16 to 0.2
     SigEvents[i] = (double)UpEvents[i]*(double)Cfg2Ratio;
   }

   cout << "Variáveis referentes à segunda configuração calculadas." << endl;

   // cfg3

   for (i=25; i<31; i++) { // No cut to 0.08
     TF1* f[i] = new TF1("fit","pol1",6350,6490);
     h[i].Fit(f[i],"l","",6350,6490);
     B[i] = f[i]->Integral(6200,6350);
     SigEvents[i] = (double)B[i]/(double)1.5;
   }

   TF1* f[31] = new TF1("fit31","pol1",6350,6490); //0.10
   h[31].Fit(f[31],"l","",6350,6490);
   B[31] = f[31]->Integral(6200,6350);
   SigEvents[31] = B[31]/1.5;
   Double_t Cfg3Ratio = (double)SigEvents[31]/(double)UpEvents[31];

   for (i=32; i<34; i++) { //0.12 and 0.14
     TF1* f[i] = new TF1("fit","pol1",6350,6490);
     h[i].Fit(f[i],"l","",6350,6490);
     B[i] = f[i]->Integral(6200,6350);
     SigEvents[i] = (double)B[i]/(double)1.5;
   }

   for (i=34; i<37; i++) { // 0.16 to 0.2
     SigEvents[i] = (double)UpEvents[i]*(double)Cfg3Ratio;
   }

   cout << "Variáveis referentes à terceira configuração calculadas." << endl;

   // cfg4

   for (i=37; i<43; i++) { // No cut to 0.08
     TF1* f[i] = new TF1("fit","pol1",6350,6490);
     h[i].Fit(f[i],"l","",6350,6490);
     B[i] = f[i]->Integral(6200,6350);
     SigEvents[i] = (double)B[i]/(double)1.5;
   }

   TF1* f[43] = new TF1("fit43","pol1",6350,6490); //0.10
   h[43].Fit(f[43],"l","",6350,6490);
   B[43] = f[43]->Integral(6200,6350);
   SigEvents[43] = B[43]/1.5;
   Double_t Cfg4Ratio = (double)SigEvents[43]/(double)UpEvents[43];

   for (i=44; i<46; i++) { //0.12 and 0.14
     TF1* f[i] = new TF1("fit","pol1",6350,6490);
     h[i].Fit(f[i],"l","",6350,6490);
     B[i] = f[i]->Integral(6200,6350);
     SigEvents[i] = (double)B[i]/(double)1.5;
   }

   for (i=46; i<49; i++) { // 0.16 to 0.2
     SigEvents[i] = (double)UpEvents[i]*(double)Cfg4Ratio;
   }

   cout << "Variáveis referentes à quarta configuração calculadas." << endl;

   ofstream out ("Expected_Background_Events_2011_MU.txt");
   out << "Número de eventos no sideband e número esperado de eventos na região do sinal: " << endl;
   out << " " << endl;

   out << "Primeira configuração do PID" << endl;
   for (i=1; i<13; i++) {
   out << "Corte de número: " << i << " Eventos no Upper Sideband: " << UpEvents[i] << " Eventos esperados na região do sinal: " << SigEvents[i];
   }

   out << "Segunda configuração do PID" << endl;
   for (i=13; i<25; i++) {
   out << "Corte de número: " << i << " Eventos no Upper Sideband: " << UpEvents[i] << " Eventos esperados na região do sinal: " << SigEvents[i];
   }

   out << "Terceira configuração do PID" << endl;
   for (i=25; i<37; i++) {
   out << "Corte de número: " << i << " Eventos no Upper Sideband: " << UpEvents[i] << " Eventos esperados na região do sinal: " << SigEvents[i];
   }

   out << "Quarta configuração do PID" << endl;
   for (i=37; i<49; i++) {
   out << "Corte de número: " << i << " Eventos no Upper Sideband: " << UpEvents[i] << " Eventos esperados na região do sinal: " << SigEvents[i];
   }
 
   ifstream in("Expected_Background_Events_2011_MU.txt");
   out.close();

   cout << "Arquivos .root e .txt salvos." << endl;


} // fim da macro