#define BHgen_cxx
// The class definition in BHgen.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.


// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// root> T->Process("BHgen.C")
// root> T->Process("BHgen.C","some options")
// root> T->Process("BHgen.C+")
//

#include <iostream>
#include <sstream>
#include <fstream>
#include <HDDM/hddm_s.hpp>

#include "BHgen.h"
#include <TH1D.h>
#include <TStyle.h>

void BHgen::Begin(TTree * /*tree*/)
{
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();
}

void BHgen::SlaveBegin(TTree * /*tree*/)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();
   h0 = new TH1D("h0", "", 200, 0, 12);
   h1 = new TH1D("h1", "", 200, 0, 12);
   h1->GetXaxis()->SetTitle("beam photon energy (GeV)");
   h1->GetYaxis()->SetTitle("Bethe Heitler cross section (ub)");
   h1->Sumw2();
   GetOutputList()->Add(h0);
   GetOutputList()->Add(h1);
}

Bool_t BHgen::Process(Long64_t entry)
{
   // The Process() function is called for each entry in the tree (or possibly
   // keyed object in the case of PROOF) to be processed. The entry argument
   // specifies which entry in the currently loaded tree is to be processed.
   // When processing keyed objects with PROOF, the object is already loaded
   // and is available via the fObject pointer.
   //
   // This function should contain the \"body\" of the analysis. It can contain
   // simple or elaborate selection criteria, run algorithms on the data
   // of the event and typically fill histograms.
   //
   // The processing can be stopped by calling Abort().
   //
   // Use fStatus to set the return value of TTree::Process().
   //
   // The return value is currently not used.

   fChain->GetEntry(entry);
   std::stringstream rootpath;
   rootpath << "root://nod28.phys.uconn.edu/Gluex/simulation/BHgen-5-2020"
            << "/BHgen_" << run << ".hddm";
   std::ifstream infile(rootpath.str().c_str());
   hddm_s::istream inhddm(infile);
   hddm_s::HDDM rec;
   while ((inhddm >> rec)) {
      hddm_s::ReactionList rea = rec.getReactions();
      hddm_s::MomentumList mom = rec.getMomenta();
      h0->Fill(mom(0).getE());
      if (mom.size() == 5) {
         double pairE = mom(2).getE() + mom(3).getE();
         double pairPx = mom(2).getPx() + mom(3).getPx();
         double pairPy = mom(2).getPy() + mom(3).getPy();
         double pairPz = mom(2).getPz() + mom(3).getPz();
         double invmass2 = pairE*pairE - pairPx*pairPx - pairPy*pairPy - pairPz*pairPz;
         h1->Fill(mom(0).getE(), rea(0).getWeight());
      }
   }
   return kTRUE;
}

void BHgen::SlaveTerminate()
{
   // The SlaveTerminate() function is called after all entries or objects
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server.

}

void BHgen::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.

   TFile fout("BHgen.root", "update");
   TIter next(GetOutputList());
   while (TObject *myobj = next()) {
      TH1 *h = dynamic_cast<TH1*>(myobj);
      if (h && h->GetEntries() > 0)
         h->Write();
   }
}