#define h1dose_cxx
// The class definition in h1dose.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("h1dose.C")
// Root > T->Process("h1dose.C","some options")
// Root > T->Process("h1dose.C+")
//

#include "h1dose.h"
#include <TH2.h>
#include <TStyle.h>
#include <string>
#include <iostream>

#define NEUTRON_MASS_GEV 0.939565378

TGraph *deq_ambient(const char *input_file);

h1dose::h1dose(TTree *tree, int nxbins, double xlow_cm, double xhigh_cm,
                            int nybins, double ylow_cm, double yhigh_cm,
                            int nzbins, double zlow_cm, double zhigh_cm,
               double duration_s)
{
   fChain = tree;
   fNbins[0] = nxbins;
   fNbins[1] = nybins;
   fNbins[2] = nzbins;
   fXlow[0] = xlow_cm;
   fXlow[1] = ylow_cm;
   fXlow[2] = zlow_cm;
   fXhigh[0] = xhigh_cm;
   fXhigh[1] = yhigh_cm;
   fXhigh[2] = zhigh_cm;
   fDuration = duration_s;

   fGraph_deq = deq_ambient("deq_ambient.dat");

   std::string name(fChain->GetName());
   std::string title(fChain->GetTitle());
   if (nxbins > 1) {
      name += "_dose";
      title += " dose (mrem/hr) as function of x (cm)";
      fH1dose = new TH1D(name.c_str(),title.c_str(),nxbins,xlow_cm,xhigh_cm);
   }
   else if (nybins > 1) {
      name += "_dose";
      title += " dose (mrem/hr) as function of y (cm)";
      fH1dose = new TH1D(name.c_str(),title.c_str(),nybins,ylow_cm,yhigh_cm);
   }
   else {
      name += "_dose";
      title += " dose (mrem/hr) as function of z (cm)";
      fH1dose = new TH1D(name.c_str(),title.c_str(),nzbins,zlow_cm,zhigh_cm);
   }
   fH1dose->Sumw2();
}

h1dose::~h1dose()
{
   delete fGraph_deq;
}

void h1dose::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 h1dose::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();

}

Bool_t h1dose::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.
   // It can be passed to either h1dose::GetEntry() or TBranch::GetEntry()
   // to read either all or the required parts of the data. 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.

   GetEntry(entry);
   if (kind != 13)
      return kFALSE;

   double bin_area_cm2 = 1;
   for (int dim=0; dim < 3; ++dim) {
      if (fNbins[dim] > 0) {
         bin_area_cm2 *= fabs(fXhigh[dim] - fXlow[dim])/fNbins[dim];
      }
      if (xin[dim] < fXlow[dim] || xin[dim] > fXhigh[dim]) {
         return kFALSE;
      }
   }  
   double u_cm = (fNbins[0] > 1)? xin[0] :
                 (fNbins[1] > 1)? xin[1] :
                 (fNbins[2] > 1)? xin[2] : xin[2];
   
   //if (pin[3]-NEUTRON_MASS_GEV < 0.001) return kFALSE;
   double w = pow(10.,fGraph_deq->Eval(log10(pin[3] - NEUTRON_MASS_GEV)));
   fH1dose->Fill(u_cm, w / bin_area_cm2 * 1e-7 * 3600/fDuration);

   return kTRUE;
}

void h1dose::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 h1dose::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.

}