#define Analiza_cxx
#include "Analiza.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <iostream>
#include <cmath>

using namespace std;

void Analiza::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L Analiza.C+ (lepsze niz .L Analiza.C, bo program 
//                            jest kompilowany)
//      Root > Analiza 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;


  fChain->SetBranchStatus("*",0);  // disable all branches

  // activate chosen branchnames
  fChain->SetBranchStatus("fHeader.*",1);
  fChain->SetBranchStatus("evthdr.*",1);
  fChain->SetBranchStatus("evt.*",1);
  fChain->SetBranchStatus("stp.*",1);
  fChain->SetBranchStatus("trk.*",1);
  fChain->SetBranchStatus("mc.*",1);
  fChain->SetBranchStatus("stdhep.*",1);
  fChain->SetBranchStatus("thstp.*",1);
  fChain->SetBranchStatus("thslc.*",1);
  fChain->SetBranchStatus("thtrk.*",1);
  fChain->SetBranchStatus("thevt.*",1);
  
  // Nie dziala dla TChain 
  //Long64_t nentries = fChain->GetEntriesFast();
  //cout << "number of entries in fChain= " << nentries << endl;
  
  Long64_t nentries = fChain->GetEntries(); 
  //number of all entries in the chain
  cout << "number of entries in fChain= " << nentries << endl;  
  
  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;
	fChain->GetEntry(jentry);
	
	Wypelnij_mc();
	Wypelnij_stdhep();
	
   }	 // End of loop
  
  // Wydruk do pliku .ps (PostScript)
  DrawAll(); 
  // Zapisywanie histogramow do pliku .root ,do dalszej analizy
  Write("hist.root");
}	

//--------------------------------------------------------------
void Analiza::BookAll()
{
  Book(Pxmu1,"Pxmu1","Muon px (truth)",30,-5,5);  
  Book(Pymu1,"Pymu1","Muon py (truth)",30,-5,5);  
  Book(Pzmu1,"Pzmu1","Muon pz (truth)",50,-20,20); 
  Book(Emu1,"Emu1","Energy (truth)",50,0,40); 
  Book(Massmu1,"Massmu1","Muon mass",50,0,0.2); 

  Book(Pxstdhep,"Pxstdhep","Stdhep px (truth)",30,-5,5);  
  Book(Pystdhep,"Pystdhep","Stdhep py (truth)",30,-5,5);  
  Book(Pzstdhep,"Pzstdhep","Stdhep pz (truth)",50,-20,20); 
  Book(Estdhep,"Estdhep","Energy stdhep (truth)",50,0,40); 
 
}
//---------------------------------------------------------------

void Analiza::Book(TH1D* &hist, const string& name, const string& title, Int_t nbins,const Double_t *binarr){
// Zmienna szerokosc przedzialow w histogramach
//
  const char* hisname = name.c_str();
  const char* histitle = title.c_str();
  hist = new TH1D(hisname,histitle,nbins,binarr);
}

//---------------------------------------------------------------

void Analiza::Book(TH1D* &hist, const string& name, const string& title, Int_t nbinsx, Axis_t xlow, Axis_t xup){
// Stala szerokosc przedzialow w histogramach
//
  const char* hisname = name.c_str();
  const char* histitle = title.c_str();
  hist = new TH1D(hisname,histitle,nbinsx,xlow,xup);
}

//---------------------------------------------------------------

void Analiza::Book(vector<TH1D*> &hist, Int_t howMany, const string& name, const string& title, Int_t nbinsx, Axis_t xlow, Axis_t xup){
//
 for(Int_t ii=0; ii<howMany; ++ii){
   const char* hisname = (name+FourMom[ii]).c_str();
   const char* histitle = (title+FourMom[ii]).c_str();
   hist.push_back(new TH1D(hisname,histitle,nbinsx,xlow,xup));
 }
}

//--------------------------------------------------------------
void Analiza::DrawAll()
// Tworzenie gotowych do wydruku rysunkow zapisywanych w formacie PostScript
{

TCanvas *canvas1 = new TCanvas("canvas1","canvas1",150,10,700,700);

  canvas1->Clear();
  canvas1->Divide(2,2,0.005,0.005);
  canvas1->cd(1);
  Pxmu1->SetXTitle("px mionu [GeV]");
  Pxmu1->SetYTitle("Liczba przypadkow");
  Pxmu1->Draw();
  canvas1->cd(2);
  Pymu1->SetXTitle("py mionu [GeV]");
  Pymu1->SetYTitle("Liczba przypadkow");
  Pymu1->Draw();
  canvas1->cd(3);
  Pzmu1->SetXTitle("pz mionu [GeV]");
  Pzmu1->SetYTitle("Liczba przypadkow");
  Pzmu1->Draw();
  canvas1->cd(4);
  Emu1->SetXTitle("Energia mionu [GeV]");
  Emu1->SetYTitle("Liczba przypadkow");
  Emu1->Draw();

  canvas1->Update();
  canvas1->Print("All.ps(");

  canvas1->Print("mion1.eps");

  canvas1->Clear();
  canvas1->Divide(2,2,0.005,0.005);
  canvas1->cd(1);
  Massmu1->SetXTitle("masa mionu [GeV]");
  Massmu1->SetYTitle("Liczba przypadkow");
  Massmu1->Draw();
  canvas1->Update();

  canvas1->Print("All.ps");

  canvas1->Clear();
  canvas1->Divide(2,2,0.005,0.005);
  canvas1->cd(1);
  Pxstdhep->SetXTitle("px stdhep [GeV]");
  Pxstdhep->SetYTitle("Liczba przypadkow");
  Pxstdhep->Draw();
  canvas1->cd(2);
  Pystdhep->SetXTitle("py stdhep [GeV]");
  Pystdhep->SetYTitle("Liczba przypadkow");
  Pystdhep->Draw();
  canvas1->cd(3);
  Pzstdhep->SetXTitle("pz stdhep [GeV]");
  Pzstdhep->SetYTitle("Liczba przypadkow");
  Pzstdhep->Draw();
  canvas1->cd(4);
  Estdhep->SetXTitle("Energia stdhep [GeV]");
  Estdhep->SetYTitle("Liczba przypadkow");
  Estdhep->Draw();

  canvas1->Update();
 
  canvas1->Print("All.ps)");

}

//--------------------------------------------------------------
void Analiza::Wypelnij_stdhep()
{

  Double_t pxstdhep;
  Double_t pystdhep;
  Double_t pzstdhep;
  Double_t estdhep;

  for(Int_t is = 0; is<stdhep_; ++is){
	if(stdhep_IstHEP[is]==1 && stdhep_IdHEP[is]!=28 ){
	  pxstdhep=stdhep_p4[is][0];
	  pystdhep=stdhep_p4[is][1];
	  pzstdhep=stdhep_p4[is][2];
	  estdhep=fabs(stdhep_p4[is][3]);
	  if(estdhep>0 ){
		Pxstdhep->Fill(pxstdhep);
		Pystdhep->Fill(pystdhep);
		Pzstdhep->Fill(pzstdhep);
		Estdhep->Fill(estdhep);
	  }            		 
	}
  } 
}	

//--------------------------------------------------------------
void Analiza::Wypelnij_mc()
{
     for(Int_t ie = 0; ie<mc_; ++ie){      

	   Double_t px=mc_p4mu1[ie][0];
	   Double_t py=mc_p4mu1[ie][1];
	   Double_t pz=mc_p4mu1[ie][2];
	   Double_t e=fabs(mc_p4mu1[ie][3]);     

	   if(e>0){
		 Pxmu1->Fill(px);
		 Pymu1->Fill(py);
		 Pzmu1->Fill(pz);
		 Emu1->Fill(e);
		 Massmu1->Fill(sqrt(e*e-px*px-py*py-pz*pz));
	   } 

	 }

}


//--------------------------------------------------------------
void Analiza::Write(char* filename){
// Zapisywanie gotowych histogramow do pliku 
TFile *file1 = new TFile(filename,"RECREATE");
if ( file1->IsOpen() ) cout << "File opened successfully" << endl;

// Jesli plik otwarty jest zanim histogramy sa tworzone, 
// to polecenie    
// file1->Write();
// zapisuje je wszystkie na raz do pliku

  Pxmu1->Write();
  Pymu1->Write();
  Pzmu1->Write();
  Emu1->Write();

  file1->Print();
  file1->Close();
}