#define EsempioAnalisi_cxx
#include "EsempioAnalisi.h"
#include <TH2.h>
#include <TF1.h>
#include <TStyle.h>
#include <TCanvas.h>

void EsempioAnalisi::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L EsempioAnalisi.C
//      Root > EsempioAnalisi 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

  gStyle->SetOptStat(1111111);
  gStyle->SetOptFit(11111);

  if (fChain == 0) return;
  
  TH1F *hTDCa0 = new TH1F("hTDCa0","hTDCa0",100,0,1170);
  TH1F *hTDCb0 = new TH1F("hTDCb0","hTDCb0",100,0,1170);

  TH1F *hTDCaZoom = new TH1F("hTDCaZoom","hTDCaZoom",100,0,1100);
  TH1F *hTDCbZoom = new TH1F("hTDCbZoom","hTDCbZoom",100,0,1100);

  TH1F *hTDCaCut = new TH1F("hTDCaCut","hTDCaCut",1000,0,1170);
  TH1F *hTDCbCut = new TH1F("hTDCbCut","hTDCbCut",1000,0,1170);

  TH1F *hTDCa = new TH1F("hTDCa","hTDCa",100,0,1170);
  TH1F *hTDCb = new TH1F("hTDCb","hTDCb",100,0,1170);

  TH2F *hTDCab = new TH2F("hTDCab","hTDaCb",100,0,1170,100,0,1170);
  TH2F *hTDCabCut = new TH2F("hTDCabCut","hTDaCb",100,0,1170,100,0,1170);

  TF1 *fCut = new TF1("fCut","[0]*x+[1]",0,1000);
  fCut->SetParameter(0,-4./3.);
  fCut->SetParameter(1,400);

  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;

    hTDCa0 ->Fill(tdca);
    hTDCb0 ->Fill(tdcb);
    
    hTDCaZoom ->Fill(tdca);
    hTDCbZoom ->Fill(tdcb);
    
    if(tdca > 1100) continue;
    if(tdca < 0) continue;

    if(tdcb > 1100) continue;
    if(tdcb < 0) continue;

    hTDCaCut ->Fill(tdca);
    hTDCbCut ->Fill(tdcb);

    hTDCa ->Fill(tdca);
    hTDCb ->Fill(tdcb);

    hTDCab ->Fill(tdca,tdcb);
	
    if(tdcb < (-4./3. * tdca + 400)) continue;
    hTDCabCut ->Fill(tdca,tdcb);
  

  }

  TCanvas *c0 = new TCanvas("c0","c0",800,400);
  c0->Divide(2,1);
  c0->cd(1);
  hTDCa0 ->DrawCopy();
  c0->cd(2);
  hTDCb0 ->DrawCopy();

  TCanvas *c1 = new TCanvas("c1","c1",800,400);
  c1->Divide(2,1);
  c1->cd(1);
  hTDCaZoom ->DrawCopy();
  c1->cd(2);
  hTDCbZoom ->DrawCopy();
    
  TCanvas *c2 = new TCanvas("c2","c2",800,400);
  c2->Divide(2,1);
  c2->cd(1);
  hTDCaCut ->DrawCopy();
  c2->cd(2);
  hTDCbCut ->DrawCopy();

  TCanvas *c3 = new TCanvas("c3","c3",800,400);
  c3->Divide(2,1);
  c3->cd(1);
  hTDCa ->DrawCopy();
  c3->cd(2);
  hTDCb ->DrawCopy();

  TCanvas *c4 = new TCanvas("c4","c4",800,400);
  c4->Divide(2,1);
  c4->cd(1);
  hTDCab ->DrawCopy("colz");
  fCut->Draw("same");
  c4->cd(2);
  hTDCabCut ->DrawCopy("colz");
  
  // Analisi sulla proiezione 

  TH1D *hTDCaP = (TH1D*)hTDCabCut->ProjectionX("hTDCaP");
  TH1D *hTDCbP = (TH1D*)hTDCabCut->ProjectionY("hTDCbP");

  TCanvas *cPro = new TCanvas("cPro","cPro",800,400);
  cPro->Divide(2,1);
  cPro->cd(1);
  hTDCaP->DrawCopy();
  cPro->cd(2);
  hTDCbP->DrawCopy();
  
  // fit e controllo bonta fit

  TF1 *fFit1 = new TF1("fFit1","[0]*TMath::Landau(x,[1],[2])");
  fFit1->SetParameters(1000,580,180);
  fFit1->SetParLimits(0,0,10000);
  fFit1->SetParLimits(2,0,500);
  
  TF1 *fFit2 = new TF1("fFit2","[0]*TMath::Landau(x,[1],[2])");
  fFit2->SetParameters(1000,580,180);
  fFit2->SetParLimits(0,0,10000);
  fFit2->SetParLimits(2,0,500);
  
 
  
  TCanvas *cFitPro = new TCanvas("cFitPro","cFitPro",800,1200);
  cFitPro->Divide(2,3);
  cFitPro->cd(1);
  hTDCaP->DrawCopy("E");
  hTDCaP->Fit("fFit1");
  cFitPro->cd(2);
  hTDCbP->DrawCopy("E");
  hTDCbP->Fit("fFit2");
 
  TH1F *hRes1 = new TH1F("hRes1","hRes1",100,0,1170);
  TH1F *hRes2 = new TH1F("hRes2","hRes2",100,0,1170);

  TH1F *hRes01 = new TH1F("hRes01","hRes01",20,-5,5);
  TH1F *hRes02 = new TH1F("hRes02","hRes02",20,-5,5);

  Float_t err = -1;
  Float_t resid = -1;

  for (int i=1; i <= hTDCaP->GetXaxis()->GetNbins(); i++) {

    hRes1->SetBinContent(i,-999);
    
    err = hTDCaP->GetBinError(i); 
    
   if(err>0){ 
      resid = (hTDCaP->GetBinContent(i)-fFit1->Eval(hTDCaP->GetBinCenter(i)))/hTDCaP->GetBinError(i);
      
      hRes1->SetBinContent(i,resid);
      hRes01->Fill(resid);
      //     cout<<resid<<endl;
    }          
    else 
      continue;

  }
    
  resid = -1;

  for (int i=1; i <= hTDCbP->GetXaxis()->GetNbins(); i++) {
    hRes2->SetBinContent(i,-999);
    
    err = hTDCbP->GetBinError(i); 
          
    if(err>0)  {
      
      resid = (hTDCbP->GetBinContent(i)-fFit2->Eval(hTDCbP->GetBinCenter(i)))/hTDCbP->GetBinError(i);
      
      hRes2->SetBinContent(i,resid);
      hRes02->Fill(resid);

      // cout<<"2: "<<resid<<" err: "<<err<<endl;
      // cout<<hTDCbP->GetBinContent(i)<<" "<<fFit2->Eval(hTDCbP->GetBinCenter(i))<<" "<<hTDCbP->GetBinCenter(i)<<endl;

    }
    else 
      continue;
  }

  hRes1->GetYaxis()->SetRangeUser(-10,10);
  hRes2->GetYaxis()->SetRangeUser(-10,10);

  cFitPro->cd(3);
  hRes1->DrawCopy("E");
  cFitPro->cd(4);
  hRes2->DrawCopy("E");


  cFitPro->cd(5);
  hRes01->DrawCopy("E");
  hRes01->Fit("gaus");
  cFitPro->cd(6);
  hRes02->DrawCopy("E");
  hRes02->Fit("gaus");

}