#define DData_cxx
#include "DData.h"

#include <TCanvas.h>
#include <TH1F.h>
#include <TLorentzVector.h>

#include <iostream>

// #if !defined(__CINT__)
// ClassImp(DData);
// #endif


//==========================================================================
DData::DData(TString name, TString title, TTree *tree) : TNamed(name, title), fChain(nullptr), 
                                                         fHistoList(nullptr), fHistosOutFile(nullptr),
                                                         fVerbose(0) 
{
	// ctor
	// if parameter tree is not specified (or zero), connect the file
	// used to generate this class and read the Tree.
	if (tree == 0) {
      TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("event.root");
      if (!f || !f->IsOpen()) {
         f = new TFile("event.root");
      }
      f->GetObject("h100",tree);
   }
   Init(tree);
}

//==========================================================================
DData::~DData()
{
	// dtor
	if (fChain)
   	delete fChain->GetCurrentFile();
   if (fHistosOutFile) {
      fHistosOutFile->Write(); 
      // fHistosOutFile->Close();
   //    fHistoList->Delete(); 
   //    delete fHistoList; 
   }
}

//==========================================================================
Double_t DData::ChargedEnergy()
{
   // calculates the charged energy assuming all charged particles
   TLorentzVector mom4; 
   Double_t ene = 0.0; 
   for (int index = 0; index < Npac; index++) {
      mom4.SetXYZM(Papx[index], Papy[index], Papz[index], Mass(Paid[index]));
      ene += mom4.E();
   }
   return ene / Ecm;
}

//==========================================================================
void DData::CreateHistograms()
{
   // open output file and create histograms

   fHistosOutFile = new TFile("histos.root", "RECREATE", "Delphi basic histograms");
   fHistoList = new TList(); 
   fHistoList->Add(new TH1F("ChargedMultiplicity", "Charged multiplicity",              100,   0.0, 100.0)); 
   fHistoList->Add(new TH1F("ChargedEnergy",       "Charged energy / ECM",              100,   0.0,   1.2));
   fHistoList->Add(new TH1F("Thrust",              "Event Thrust",                      100,   0.5,   1.0));
   fHistoList->Add(new TH1F("Oblateness",          "Event Oblateness",                  100,   0.0,   1.0));
   fHistoList->Add(new TH1F("Major",               "Event Major",                       100,   0.0,   1.0));
   fHistoList->Add(new TH1F("Sphericity",          "Event Sphericity",                  100,   0.0,   1.0));
   fHistoList->Add(new TH1F("Aplanarity",          "Event Aplanarity",                  100,   0.0,   1.0));
   fHistoList->Add(new TH1F("RPhi",                "Impact R/Phi(pT)",                  100,  -5.0,   5.0));
   fHistoList->Add(new TH1F("Z",                   "Impact Z(pT)",                      110, -11.0,  11.0));
   fHistoList->Add(new TH1F("Momentum",            "p[GeV]",                            500,   0.0,   5.0));
   fHistoList->Add(new TH1F("Theta08",             "Theta p < 0.8 GeV",                 180,   0.0, 180.0));
   fHistoList->Add(new TH1F("Theta082",            "Theta 0.8 < p < 2 GeV",             180,   0.0, 180.0));
   fHistoList->Add(new TH1F("Theta25",             "Theta 2 < p < 5 GeV",               180,   0.0, 180.0));
   fHistoList->Add(new TH1F("Theta5",              "Theta p > 5 GeV",                   180,   0.0, 180.0));
   fHistoList->Add(new TH1F("Phi",                 "Phi",                               360,-180.0, 180.0));
   fHistoList->Add(new TH1F("pLong",               "p Long [GeV]",                      200,   0.0,  50.0));
   fHistoList->Add(new TH1F("pIn",                 "pt In [GeV]",                       200,   0.0,  20.0));
   fHistoList->Add(new TH1F("pOut",                "pt Out [GeV]",                      200,   0.0,   5.0));
}

//==========================================================================
Int_t DData::Cut(Long64_t /*entry*/)
{
	// This function may be called from Loop.
	// returns  1 if entry is accepted.
	// returns -1 otherwise.
	return 1;
}

//==========================================================================
void  DData::FillHistos()
{
   // Fill the histograms
   if (!fHistoList)
      CreateHistograms(); 
   (dynamic_cast<TH1F*>(fHistoList->FindObject("ChargedMultiplicity")))->Fill(Npac);  
   (dynamic_cast<TH1F*>(fHistoList->FindObject("ChargedEnergy")))->Fill(ChargedEnergy());  
   (dynamic_cast<TH1F*>(fHistoList->FindObject("Thrust")))->Fill(Thrval1);  
   (dynamic_cast<TH1F*>(fHistoList->FindObject("Oblateness")))->Fill(Thrval2 - Thrval3);  
   (dynamic_cast<TH1F*>(fHistoList->FindObject("Major")))->Fill(Thrval2);  
   (dynamic_cast<TH1F*>(fHistoList->FindObject("Sphericity")))->Fill((Sphval2 + Sphval3) * 3 / 2);  
   (dynamic_cast<TH1F*>(fHistoList->FindObject("Aplanarity")))->Fill(Sphval3 * 3 / 2);  
   TLorentzVector mom4;
   for (Int_t part = 0; part < Npac; part++) {
      (dynamic_cast<TH1F*>(fHistoList->FindObject("RPhi")))->Fill(Rvtx[part]);  
      (dynamic_cast<TH1F*>(fHistoList->FindObject("Z")))->Fill(Zvtx[part]);  
      mom4.SetXYZM(Papx[part], Papy[part], Papz[part], Mass(Paid[part]));
      (dynamic_cast<TH1F*>(fHistoList->FindObject("Momentum")))->Fill(mom4.P());  
      if (mom4.P() < 0.8)
         (dynamic_cast<TH1F*>(fHistoList->FindObject("Theta08")))->Fill(TMath::RadToDeg() *mom4.Theta());  
      else if (mom4.P() >= 0.8 && mom4.P() < 2.0)
         (dynamic_cast<TH1F*>(fHistoList->FindObject("Theta082")))->Fill(TMath::RadToDeg() *mom4.Theta());  
      else if (mom4.P() >= 2.0 && mom4.P() < 5.0)
         (dynamic_cast<TH1F*>(fHistoList->FindObject("Theta25")))->Fill(TMath::RadToDeg() *mom4.Theta());  
      else 
         (dynamic_cast<TH1F*>(fHistoList->FindObject("Theta5")))->Fill(TMath::RadToDeg() *mom4.Theta());  
      (dynamic_cast<TH1F*>(fHistoList->FindObject("Phi")))->Fill(TMath::RadToDeg() *mom4.Phi());  
      TVector3 p = mom4.Vect(); 
      TVector3 vThrust(Thrvec1x, Thrvec1y, Thrvec1z);
      (dynamic_cast<TH1F*>(fHistoList->FindObject("pLong")))->Fill(p.Dot(vThrust));
      TVector3 vMajor(Thrvec2x, Thrvec2y, Thrvec2z);
      (dynamic_cast<TH1F*>(fHistoList->FindObject("pIn")))->Fill(p.Dot(vMajor));
      TVector3 vMinor(Thrvec3x, Thrvec3y, Thrvec3z);
      (dynamic_cast<TH1F*>(fHistoList->FindObject("pOut")))->Fill(p.Dot(vMinor));
   }
}


//==========================================================================
Int_t DData::GetEntry(Long64_t entry)
{
	// Read contents of event #entry [0,fChain->GetEntriesFast()].
   if (!fChain) return 0;
   return fChain->GetEntry(entry);
}

//==========================================================================
Double_t DData::Mass(Int_t pdg)
{
	// returns mass of particle with pdg code pdg
	// for non pdg code see myjob.car
   	Double_t mass = 0.0; 
   	if (fMassTable.size() == 0) {
   		fMassTable[11]    = 0.0005109989461; // electron
      	fMassTable[-11]   = 0.0005109989461; // positron
      	fMassTable[13]    = 0.1056583745;    // mu+
      	fMassTable[-13]   = 0.1056583745;    // mu-
      	fMassTable[22]    = 0.0;             // photon
      	fMassTable[111]   = 0.1349770;       // pi0
      	fMassTable[130]   = 0.497611;        // K0L
      	fMassTable[211]   = 0.13957061;      // pi+
      	fMassTable[-211]  = 0.13957061;      // pi-
      	fMassTable[310]   = 0.497611;        // K0S
      	fMassTable[311]   = 0.497611;        // K0
      	fMassTable[321]   = 0.493677;        // K+
      	fMassTable[-321]  = 0.493677;        // K-
      	fMassTable[2121]  = 0.9395654133;    // neutron
      	fMassTable[2212]  = 0.9382720813;    // proton
      	fMassTable[-2212] = 0.9382720813;    // anti proton
      	fMassTable[3122]  = 1.115683;        // lambda
   	}
   	mass = fMassTable[pdg]; 
   	if (mass == 0.0) {
         if (fVerbose == 1)
            std::cout << "Warning: GetMass --> mass for  pdg = " << pdg << " not found! mass set to pion mass" << std::endl; 
      	mass = 0.13957061;
   	}
   	return mass; 
}

//==========================================================================
void DData::Init(TTree *tree)
{
   // The Init() function is called when the selector needs to initialize
   // a new tree or chain. Typically here the branch addresses and branch
   // pointers of the tree will be set.
   // It is normally not necessary to make changes to the generated
   // code, but the routine can be extended by the user if needed.
   // Init() will be called many times when running on PROOF
   // (once per file to be processed).

   // Set branch addresses and branch pointers
   if (!tree) return;
   fChain = tree;
   fCurrent = -1;
   fChain->SetMakeClass(1);

   fChain->SetBranchAddress("Ecm", &Ecm, &b_Ecm);
   fChain->SetBranchAddress("Sphvec1x", &Sphvec1x, &b_Sphvec1x);
   fChain->SetBranchAddress("Sphvec1y", &Sphvec1y, &b_Sphvec1y);
   fChain->SetBranchAddress("Sphvec1z", &Sphvec1z, &b_Sphvec1z);
   fChain->SetBranchAddress("Sphvec2x", &Sphvec2x, &b_Sphvec2x);
   fChain->SetBranchAddress("Sphvec2y", &Sphvec2y, &b_Sphvec2y);
   fChain->SetBranchAddress("Sphvec2z", &Sphvec2z, &b_Sphvec2z);
   fChain->SetBranchAddress("Sphvec3x", &Sphvec3x, &b_Sphvec3x);
   fChain->SetBranchAddress("Sphvec3y", &Sphvec3y, &b_Sphvec3y);
   fChain->SetBranchAddress("Sphvec3z", &Sphvec3z, &b_Sphvec3z);
   fChain->SetBranchAddress("Sphval1", &Sphval1, &b_Sphval1);
   fChain->SetBranchAddress("Sphval2", &Sphval2, &b_Sphval2);
   fChain->SetBranchAddress("Sphval3", &Sphval3, &b_Sphval3);
   fChain->SetBranchAddress("Thrvec1x", &Thrvec1x, &b_Thrvec1x);
   fChain->SetBranchAddress("Thrvec1y", &Thrvec1y, &b_Thrvec1y);
   fChain->SetBranchAddress("Thrvec1z", &Thrvec1z, &b_Thrvec1z);
   fChain->SetBranchAddress("Thrvec2x", &Thrvec2x, &b_Thrvec2x);
   fChain->SetBranchAddress("Thrvec2y", &Thrvec2y, &b_Thrvec2y);
   fChain->SetBranchAddress("Thrvec2z", &Thrvec2z, &b_Thrvec2z);
   fChain->SetBranchAddress("Thrvec3x", &Thrvec3x, &b_Thrvec3x);
   fChain->SetBranchAddress("Thrvec3y", &Thrvec3y, &b_Thrvec3y);
   fChain->SetBranchAddress("Thrvec3z", &Thrvec3z, &b_Thrvec3z);
   fChain->SetBranchAddress("Thrval1", &Thrval1, &b_Thrval1);
   fChain->SetBranchAddress("Thrval2", &Thrval2, &b_Thrval2);
   fChain->SetBranchAddress("Thrval3", &Thrval3, &b_Thrval3);
   fChain->SetBranchAddress("Npa", &Npa, &b_Npa);
   fChain->SetBranchAddress("Npac", &Npac, &b_Npac);
   fChain->SetBranchAddress("Paid", Paid, &b_Paid);
   fChain->SetBranchAddress("Pafl", Pafl, &b_Pafl);
   fChain->SetBranchAddress("Pav0d", Pav0d, &b_Pav0d);
   fChain->SetBranchAddress("Pav0m", Pav0m, &b_Pav0m);
   fChain->SetBranchAddress("Rvtx", Rvtx, &b_Rvtx);
   fChain->SetBranchAddress("Zvtx", Zvtx, &b_Zvtx);
   fChain->SetBranchAddress("Papx", Papx, &b_Papx);
   fChain->SetBranchAddress("Papy", Papy, &b_Papy);
   fChain->SetBranchAddress("Papz", Papz, &b_Papz);
   fChain->SetBranchAddress("Paje", Paje, &b_Paje);
   fChain->SetBranchAddress("Pani", Pani, &b_Pani);
   fChain->SetBranchAddress("Njer", &Njer, &b_Njer);
   fChain->SetBranchAddress("Tgenr", &Tgenr, &b_Tgenr);
   fChain->SetBranchAddress("Dminr", &Dminr, &b_Dminr);
   fChain->SetBranchAddress("Jepx", Jepx, &b_Jepx);
   fChain->SetBranchAddress("Jepy", Jepy, &b_Jepy);
   fChain->SetBranchAddress("Jepz", Jepz, &b_Jepz);
   fChain->SetBranchAddress("Jee", Jee, &b_Jee);
   fChain->SetBranchAddress("Jem", Jem, &b_Jem);
   fChain->SetBranchAddress("Jep", Jep, &b_Jep);
   Notify();
}

//==========================================================================
Long64_t DData::LoadTree(Long64_t entry)
{
	// Set the environment to read one entry
   	if (!fChain) return -5;
   	Long64_t centry = fChain->LoadTree(entry);
   	if (centry < 0) return centry;
   	if (fChain->GetTreeNumber() != fCurrent) {
    	 fCurrent = fChain->GetTreeNumber();
     	 Notify();
   	}
   	return centry;
}

//==========================================================================
void DData::Loop()
{
//   In a ROOT session, you can do:
//      root> .L DData.C
//      root> DData data
//      root> data.GetEntry(12); // Fill t data members with entry number 12
//      root> data.Show();       // Show values of entry 12
//      root> data.Show(16);     // Read and show values of entry 16
//      root> data.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();

   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;
      FillHistos();  
   }
}

//==========================================================================
Bool_t DData::Notify()
{
   // The Notify() function is called when a new file is opened. This
   // can be either for a new TTree in a TChain or when when a new TTree
   // is started when using PROOF. It is normally not necessary to make changes
   // to the generated code, but the routine can be extended by the
   // user if needed. The return value is currently not used.

   return kTRUE;
}

 //==========================================================================
void DData::Plot() const
{
   // plot histograms 
   Int_t nHistos = TMath::Sqrt(fHistoList->GetSize()) + 1; 
   TCanvas *can = new TCanvas("can","Delphi control histos");
   can->SetCanvasSize(1000, 1000);
   can->SetWindowSize(500, 500);
   can->Divide (nHistos, nHistos); 

   TIter next(fHistoList); 
   Int_t npad = 1;    
   while ( TObject *h = next() ) {
      can->cd(npad++)->SetLogy();
      dynamic_cast<TH1F*>(h)->Draw();

      KolmogorovTest()
   } 
}

//==========================================================================
void DData::PrintMomentum() const
{
	// print ntuple 3-momentum components
	for (int index = 0; index < Npa; index++)
    	std:cout << "Info: " << index << ": Px = " << Papx[index] << " Py = " << Papy[index] << " Pz = " << Papz[index] << std::endl;
}

//==========================================================================
void DData::PrintPID() const 
{
	// print PID of all particles in the event
	for (int index = 0; index < Npa; index++)
		std::cout << "Info: PID(" << index << ") = " << Paid[index] << std::endl; 
}

//==========================================================================
void DData::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
   if (!fChain) return;
   fChain->Show(entry);
}