#define HP_Ge_cxx
#include "HP_Ge.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

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

   Long64_t nentries = fChain->GetEntriesFast();
 
      Double_t Raw   = 0 ;
      Double_t LaBr3 = 0 ;
      
      Double_t nEnergyBins = 3000;
        
   TH1F *Energy_raw = new TH1F("Energy_raw","Energy_raw",nEnergyBins,-100,2900);   
   TH1F *Energy_Calib = new TH1F("Energy_Calib","Energy_Calib",nEnergyBins,-100,2900);  
 
        TFile outtracker("outtracker.root","recreate");
 	TTree *EnergyTree = new TTree("E_Calib","E_Calib") ;
	EnergyTree->Branch("Energy_raw", &Raw) ;
	EnergyTree->Branch("Energy_Calib", &LaBr3) ;  
	

   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;
      

      Raw   = LaBr3Raw[3];
      LaBr3 = channel_Energy(Raw);
      
      //-78.61 + (0.6816*LaBr3Raw[3]) + (1.043*TMath::Power(10,-6)*TMath::Power(LaBr3Raw[3],2));
       
       
      Energy_raw->Fill(Raw);  
      Double_t dE = Energy_Calib->GetBinWidth(1) / 2.0; // half of the "fixed" bin width   
      Energy_Calib->Fill(LaBr3+dE);
      
       EnergyTree  -> Fill();   
      // if (Cut(ientry) < 0) continue;
   }
  // Energy_Calib->Draw();
  
//cout<<Energy_Calib->GetBinWidth(1);
 TCanvas *c1 = new TCanvas("c1","Energy spectra of HP Ge",800,600);;
 c1->Divide(1,2);  
 
 gStyle->SetOptStat(2211);
 gStyle->SetOptFit(1111);
 
 c1->cd(1);
 Energy_raw->Draw();
 c1->cd(2); 
 Energy_Calib->Draw();
 
    // Saving the hitograms in root file

   Energy_raw->Write();
   Energy_Calib->Write();
   
   
   EnergyTree  ->  Write();
  
   outtracker.Close(); 
   

}



Double_t HP_Ge::channel_Energy(Double_t x)
{
  Double_t v;
  //v = -78.61 + (0.6816*LaBr3Raw[3]) + (1.043*TMath::Power(10,-6)*TMath::Power(LaBr3Raw[3],2));     
  v = 0.6861 * x - 82.4698; // "linear scaling" function example
  return v;
}

/*

void ScaleXaxis(TH1 *h, Double_t (*Scale)(Double_t))
{
  if (!h) return; // just a precaution
  ScaleAxis(h->GetXaxis(), Scale);
  return;
}

void ScaleAxis(TAxis *a, Double_t (*Scale)(Double_t))
{
  if (!a) return; // just a precaution
  if (a->GetXbins()->GetSize())
    {
      // an axis with variable bins
      // note: bins must remain in increasing order, hence the "Scale"
      // function must be strictly (monotonically) increasing
      TArrayD X(*(a->GetXbins()));
      for(Int_t i = 0; i < X.GetSize(); i++) X[i] = Scale(X[i]);
      a->Set((X.GetSize() - 1), X.GetArray()); // new Xbins
    }
  else
    {
      // an axis with fix bins
      // note: we modify Xmin and Xmax only, hence the "Scale" function
      // must be linear (and Xmax must remain greater than Xmin)
      a->Set( a->GetNbins(),
              Scale(a->GetXmin()), // new Xmin
              Scale(a->GetXmax()) ); // new Xmax
    }
  return;
}*/