#define ReadLHE_cxx
#include "ReadLHE.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

#include "TLorentzVector.h"

#include <cstdlib>
#include <cmath>
#include <iostream>

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

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      if (jentry > 1) break; // just the first 2 entries
#if 0 /* 0 or 1 */
      nb = fChain->GetEntry(jentry); nbytes += nb;
#else /* 0 or 1 */
      nb = 0;
      nb += b_Particle_->GetEntry(ientry);
      nb += b_Particle_PID->GetEntry(ientry);
      nb += b_Particle_Px->GetEntry(ientry);
      nb += b_Particle_Py->GetEntry(ientry);
      nb += b_Particle_Pz->GetEntry(ientry);
      nb += b_Particle_E->GetEntry(ientry);
      nbytes += nb;
#endif /* 0 or 1 */
      // if (Cut(ientry) < 0) continue;
      std::cout << "nb = " << nb << std::endl;
      int quakr_count = 0;
      for (int im = 0; im < Particle_; im++) {
        std::cout << im << "\tParticle_Py = " << Particle_E[im] << std::endl;
        std::cout << im << "\tParticle_PID = " << Particle_PID[im] << std::endl;
        if ((std::abs(Particle_PID[im]) < 9) || (std::abs(Particle_PID[im]) == 21)) {
          quakr_count++;
          std::cout << "quakr_count = " << quakr_count << std::endl;
          if (quakr_count==1) {
            std::cout << "px, py, pz, E " << Particle_Px[im] << "\t" << Particle_Py[im] << "\t" << Particle_Pz[im] << "\t"<< Particle_E[im] << std::endl;
            std::cout << im << "\tParticle_ BEFORE Iqrk0 = " << Particle_ << std::endl;
#if 1 /* 0 or 1 */
            TLorentzVector Iqrk0(0, 0, 0, 0);
            // TLorentzVector Iqrk0(Particle_Px[im], Particle_Py[im], Particle_Pz[im], Particle_E[im]);
#else /* 0 or 1 */
            TLorentzVector Iqrk0;
            Iqrk0.SetPxPyPzE(0, 0, 0, 0);
            // Iqrk0.SetPxPyPzE(Particle_Px[im], Particle_Py[im], Particle_Pz[im], Particle_E[im]);
#endif /* 0 or 1 */
            std::cout << im << "\tParticle_ AFTER Iqrk0 = " << Particle_ << std::endl;
          }
          if (quakr_count==2)
            std::cout << "px, py, pz, E " << Particle_Px[im] << "\t" << Particle_Py[im] << "\t" << Particle_Pz[im] << "\t" << Particle_E[im] << std::endl;
          if (quakr_count==3)
            std::cout << "px, py, pz, E " << Particle_Px[im] << "\t" << Particle_Py[im] << "\t" << Particle_Pz[im] << "\t" << Particle_E[im] << std::endl;
          if (quakr_count==4)
            std::cout << "px, py, pz, E " << Particle_Px[im] << "\t" << Particle_Py[im] << "\t" << Particle_Pz[im] << "\t" << Particle_E[im] << std::endl;
        }
      }
   }
}