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

using namespace std;

Long64_t Analiza::Laduj()
{
//   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 1;


  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("shw.*",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;

  return nentries;

}

void Analiza::Loop()
{

  Long64_t nentries = Laduj(); 

  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);
	// Monte Carlo
	Wypelnij_mc();
	Wypelnij_stdhep();
        if(mc_iaction[0]==1){
	  Wypelnij_energie();
	}
    }	 // 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(Pmu1, "Pmu1","Muon p",50,0,20); 
  Book(Emu1,"Emu1","Energy (truth)",50,0,20); 
  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); 

  Book(Erec1,"Erec1","Reconstructed energy z evthdr",25,0,500000);  
  Book(Erec2,"Erec2","Reconstructed energy z evt (gev)",25,0,40);  
  Book(Erec3,"Erec3","Reconstructed energy z evt (mip)",25,0,450);  
  Book(Erec4,"Erec4","Reconstructed energy z shw (gev)",25,0,20); 
  Book(Erec4a,"Erec4a","Reconstructed energy, linCC (gev)",25,0,20); 
  Book(Erec4b,"Erec4b","Reconstructed energy, wtCC (gev)",25,0,20);
  Book(Erec4c,"Erec4c","Reconstructed energy, linNC (gev)",25,0,20); 
  Book(Erec4d,"Erec4d","Reconstructed energy, wtNC (gev)",25,0,20);
  Book(Erec5,"Erec5","Reconstructed energy z shw (mip)",25,0,450); 
  Book(Erec6,"Erec6","Reconstructed energy z trk (gev)",25,0,20); 
  Book(Erec7,"Erec7","Reconstructed energy z trk (mip)",25,0,450); 
  Book(Erec8,"Erec8","Reconstructed energy z shw+trk (gev)",25,0,20); 
  Book(Erec9,"Erec9","Reconstructed energy z shw+trk (mip)",25,0,450); 
 
  Book(MomBest,"MomBest","Best momentum",25,0,20); 
  Book(EnBest,"EnBest","Best energy",25,0,20); 


}
//---------------------------------------------------------------

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->SetNameTitle("Px mu1","");
  Pxmu1->Draw();
  canvas1->cd(2);
  Pymu1->SetXTitle("py mionu [GeV]");
  Pymu1->SetYTitle("Liczba przypadkow");
  //Pymu1->SetNameTitle("Py mu1","");
  Pymu1->Draw();
  canvas1->cd(3);
  Pzmu1->SetXTitle("pz mionu [GeV]");
  Pzmu1->SetYTitle("Liczba przypadkow");
  //Pzmu1->SetNameTitle("Pz mu1","");
  Pzmu1->Draw();
  canvas1->cd(4);
  Emu1->SetXTitle("Energia mionu [GeV]");
  Emu1->SetYTitle("Liczba przypadkow");
  //Emu1->SetNameTitle("Energia mu1","");
  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->SetNameTitle("Px mu1","");
  Massmu1->Draw();
  canvas1->cd(2);
  Pmu1->SetXTitle("ped mionu [GeV]");
  Pmu1->SetYTitle("Liczba przypadkow");
  //Pmu1->SetNameTitle("P mu1","");
  Pmu1->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->SetNameTitle("Px stdhep","");
  Pxstdhep->Draw();
  canvas1->cd(2);
  Pystdhep->SetXTitle("py stdhep [GeV]");
  Pystdhep->SetYTitle("Liczba przypadkow");
  //Pystdhep->SetNameTitle("Py stdhep","");
  Pystdhep->Draw();
  canvas1->cd(3);
  Pzstdhep->SetXTitle("pz stdhep [GeV]");
  Pzstdhep->SetYTitle("Liczba przypadkow");
  //Pzstdhep->SetNameTitle("Pz stdhep","");
  Pzstdhep->Draw();
  canvas1->cd(4);
  Estdhep->SetXTitle("Energia stdhep [GeV]");
  Estdhep->SetYTitle("Liczba przypadkow");
  //Estdhep->SetNameTitle("Energia stdhep","");
  Estdhep->Draw();

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

  canvas1->Clear();
  canvas1->Divide(2,2,0.005,0.005);
  canvas1->cd(1);
  Erec1->SetXTitle("Zrekonstruowana energia, evthdr, sigcor");
  Erec1->SetYTitle("Liczba przypadkow");
  Erec1->Draw(); 

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

  canvas1->Clear();
  canvas1->Divide(2,2,0.005,0.005);
  canvas1->cd(1);
  Erec2->SetXTitle("Zrekonstruowana energia event [GeV]");
  Erec2->SetYTitle("Liczba przypadkow");
  Erec2->Draw();
  canvas1->cd(2);
  Erec3->SetXTitle("Zrekonstruowana energia event [mip]");
  Erec3->SetYTitle("Liczba przypadkow");
  Erec3->Draw();
  canvas1->cd(3);
  Erec4->SetXTitle("Zrekonstruowana energia shower [GeV]");
  Erec4->SetYTitle("Liczba przypadkow");
  Erec4->Draw();
  canvas1->cd(4);
  Erec5->SetXTitle("Zrekonstruowana energia shower [mip]");
  Erec5->SetYTitle("Liczba przypadkow");
  Erec5->Draw();
  canvas1->Update();
 
  canvas1->Print("All.ps");

  canvas1->Clear();
  canvas1->Divide(2,2,0.005,0.005);
  canvas1->cd(1);
  Erec4a->SetXTitle("Zrekonstruowana energia, shower linCC [GeV]");
  Erec4a->SetYTitle("Liczba przypadkow");
  Erec4a->Draw();
  canvas1->cd(2);
  Erec4b->SetXTitle("Zrekonstruowana energia, shower wtCC [GeV]");
  Erec4b->SetYTitle("Liczba przypadkow");
  Erec4b->Draw();
  canvas1->cd(3);
  Erec4c->SetXTitle("Zrekonstruowana energia, shower linNC [GeV]");
  Erec4c->SetYTitle("Liczba przypadkow");
  Erec4c->Draw();
  canvas1->cd(4);
  Erec4d->SetXTitle("Zrekonstruowana energia, shower wtNC[GeV]");
  Erec4d->SetYTitle("Liczba przypadkow");
  Erec4d->Draw();
  canvas1->Update();
 
  canvas1->Print("All.ps");


  canvas1->Clear();
  canvas1->Divide(2,2,0.005,0.005);
  canvas1->cd(1);
  Erec6->SetXTitle("Zrekonstruowana energia track [GeV]");
  Erec6->SetYTitle("Liczba przypadkow");
  Erec6->Draw();
  canvas1->cd(2);
  Erec7->SetXTitle("Zrekonstruowana energia track [mip]");
  Erec7->SetYTitle("Liczba przypadkow");
  Erec7->Draw();
  canvas1->cd(3);
  MomBest->SetXTitle("Ped mionu (best: zasieg lub krzywizna) [GeV]");
  MomBest->SetYTitle("Liczba przypadkow");
  MomBest->Draw();
  canvas1->cd(4);
  EnBest->SetXTitle("Energia mionu (best) [GeV]");
  EnBest->SetYTitle("Liczba przypadkow");
  EnBest->Draw();
 
  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);
                 Pmu1->Fill(sqrt(px*px+py*py+pz*pz));
		 Massmu1->Fill(sqrt(e*e-px*px-py*py-pz*pz));
	   } 

	 }

}

//--------------------------------------------------------------
void Analiza::Wypelnij_energie()
// Histogramy: energie kaskad hadronowych, toru itp.
// Petla po przypadkach (events !)
{
 
  

  Erec1->Fill(evthdr_ph_sigcor);


  for(Int_t ie=0;ie<evt_;ie++){ 

     Erec2->Fill(evt_ph_gev[ie]);
     Erec3->Fill(evt_ph_mip[ie]);
    
        for(Int_t it=0;it<evt_ntrack[ie];it++){
           Int_t* wsk=evt_trk[ie];
              if(it!=0){
		 wsk++;
	      }
	      Erec6->Fill(trk_ph_gev[*wsk]); 
	      Erec7->Fill(trk_ph_mip[*wsk]); 
              MomBest->Fill(trk_momentum_best[*wsk]);
              Double_t muen=
		sqrt((trk_momentum_best[*wsk]*trk_momentum_best[*wsk])
		     +(0.1057*0.1057));
		EnBest->Fill(muen);  

	}
        for(Int_t it=0;it<evt_nshower[ie];it++){
           Int_t* wsk=evt_shw[ie];
              if(it!=0){
		 wsk++;
	      }
	      Erec4->Fill(shw_ph_gev[*wsk]);
	      Erec4a->Fill(shw_shwph_linCCgev[*wsk]);
	      Erec4b->Fill(shw_shwph_wtCCgev[*wsk]);
	      Erec4c->Fill(shw_shwph_linNCgev[*wsk]);
	      Erec4d->Fill(shw_shwph_wtNCgev[*wsk]);
	      Erec5->Fill(shw_ph_mip[*wsk]); 
	}
    
  }
}

//--------------------------------------------------------------
void Analiza::Rysuj()
// Prosty event display. Na zielono zaznaczone niskie sygnaly
// (ponizej 2 fotoelektronow) 
{

  Long64_t nentries = Laduj(); 
  
   gStyle->SetMarkerStyle(20);
   gStyle->SetCanvasColor(kWhite);

   TCanvas *c1 = new TCanvas("Canva1","events",150,10,1000,500);
   c1->Divide(2,1,0.005,0.005);  


  for (Long64_t jentry=0; jentry<nentries;jentry++) {
      fChain->GetEntry(jentry);

      // Jezeli przypadek typu CC
      if(mc_iaction[0]==1 && evt_>0){
        if(evt_ph_gev[0]<12){
	  c1->cd(1); 
	  fChain->SetMarkerColor(1); 
	  fChain->Draw("stp.tpos:stp.z","stp.planeview==2","",1,jentry);
	  fChain->SetMarkerColor(4);
	  fChain->Draw("stp.tpos:stp.z","stp.planeview==2 && (stp.ph0.pe+stp.ph1.pe)>2","same",1,jentry);
	  fChain->SetMarkerColor(3);
	  fChain->Draw("stp.tpos:stp.z","stp.planeview==2 && (stp.ph0.pe+stp.ph1.pe)<=2","same",1,jentry);
	  c1->cd(2);
	  fChain->SetMarkerColor(1); 
	  fChain->Draw("stp.tpos:stp.z","stp.planeview==3","",1,jentry);
	  fChain->SetMarkerColor(4);
	  fChain->Draw("stp.tpos:stp.z","stp.planeview==3 && (stp.ph0.pe+stp.ph1.pe)>2","same",1,jentry);
	  fChain->SetMarkerColor(3);
	  fChain->Draw("stp.tpos:stp.z","stp.planeview==3 && (stp.ph0.pe+stp.ph1.pe)<=2","same",1,jentry);
	  c1->Update();
	  cout << "Next(y/n) ?" << endl;
	  string next;
	  cin >> next;      
	  if(next=="n"){
	    return;
	  }   
	}  
      }
  }	 // End of loop
  

}	


//--------------------------------------------------------------
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();
  Pmu1->Write();

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