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


#include "hpolar.h"
#include <TH2.h>
#include <TStyle.h>
#include <TList.h>
#include <sstream>

void hpolar::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 hpolar::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();
   for (int i=0; i < 9; ++i) {
      std::stringstream name;
      std::stringstream title;
      name << "hintens_" << i;
      title << "beam spot intensity map, "
            << 3+i << "-" << 4+i << "GeV";
      hintens[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                            50, -10., 10., 50, -10., 10.);
      name.str("");
      title.str("");
      name << "hpolar0_" << i;
      title << "beam spot polarization map (0-90), "
            << 3+i << "-" << 4+i << "GeV";
      hpolar0[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                            50, -10., 10., 50, -10., 10.);
      name.str("");
      title.str("");
      name << "hpolar45_" << i;
      title << "beam spot polarization map (45-135), "
            << 3+i << "-" << 4+i << "GeV";
      hpolar45[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                             50, -10., 10., 50, -10., 10.);
      name.str("");
      title.str("");
      name << "h2polar0_" << i;
      title << "beam spot w**2 p**2 map (0-90), "
            << 3+i << "-" << 4+i << "GeV";
      h2polar0[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                             50, -10., 10., 50, -10., 10.);
      name.str("");
      title.str("");
      name << "h2polar45_" << i;
      title << "beam spot w**2 p**2 map (45-135), "
            << 3+i << "-" << 4+i << "GeV";
      h2polar45[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                              50, -10., 10., 50, -10., 10.);
      name.str("");
      title.str("");
      name << "h3polar0_" << i;
      title << "beam spot w**3 p**2 map (0-90), "
            << 3+i << "-" << 4+i << "GeV";
      h3polar0[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                             50, -10., 10., 50, -10., 10.);
      name.str("");
      title.str("");
      name << "h3polar45_" << i;
      title << "beam spot w**3 p**2 map (45-135), "
            << 3+i << "-" << 4+i << "GeV";
      h3polar45[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                              50, -10., 10., 50, -10., 10.);
      name.str("");
      title.str("");
      name << "h4polar0_" << i;
      title << "beam spot w**4 p**2 map (0-90), "
            << 3+i << "-" << 4+i << "GeV";
      h4polar0[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                             50, -10., 10., 50, -10., 10.);
      name.str("");
      title.str("");
      name << "h4polar45_" << i;
      title << "beam spot w**4 p**2 map (45-135), "
            << 3+i << "-" << 4+i << "GeV";
      h4polar45[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                              50, -10., 10., 50, -10., 10.);
      hintens[i]->SetDirectory(0);
      hpolar0[i]->SetDirectory(0);
      hpolar45[i]->SetDirectory(0);
      h2polar0[i]->SetDirectory(0);
      h2polar45[i]->SetDirectory(0);
      h3polar0[i]->SetDirectory(0);
      h3polar45[i]->SetDirectory(0);
      h4polar0[i]->SetDirectory(0);
      h4polar45[i]->SetDirectory(0);
      GetOutputList()->Add(hintens[i]);
      GetOutputList()->Add(hpolar0[i]);
      GetOutputList()->Add(hpolar45[i]);
      GetOutputList()->Add(h2polar0[i]);
      GetOutputList()->Add(h2polar45[i]);
      GetOutputList()->Add(h3polar0[i]);
      GetOutputList()->Add(h3polar45[i]);
      GetOutputList()->Add(h4polar0[i]);
      GetOutputList()->Add(h4polar45[i]);
   }
}

Bool_t hpolar::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.

   fReader.SetLocalEntry(entry);
   static double distanceRadCol = 76e3; // units are mm
   double x = distanceRadCol * k[1] / k[3];
   double y = distanceRadCol * k[2] / k[3];
   double phi = atan2(y,x);
   double w = *diffXS * *weight;
   double p0lab = (*polar_0_90 * 2 - 1) * cos(2 * phi) -
                  (*polar_45_135 * 2 - 1) * sin(2 * phi);
   double p45lab = (*polar_45_135 * 2 - 1) * cos(2 * phi) +
                   (*polar_0_90 * 2 - 1) * sin(2 * phi);
   int i = int(k[0] - 3.);
   if (i >= 0 && i < 9) {
      double w2 = w * w;
      double w3 = w2 * w;
      double w4 = w3 * w;
      double p0lab2 = p0lab * p0lab;
      double p45lab2 = p45lab * p45lab;
      hintens[i]->Fill(x, y, w);
      hpolar0[i]->Fill(x, y, w * p0lab);
      hpolar45[i]->Fill(x, y, w * p45lab);
      h2polar0[i]->Fill(x, y, w2 * p0lab2);
      h2polar45[i]->Fill(x, y, w2 * p45lab2);
      h3polar0[i]->Fill(x, y, w3 * p0lab2);
      h3polar45[i]->Fill(x, y, w3 * p45lab2);
      h4polar0[i]->Fill(x, y, w4 * p0lab2);
      h4polar45[i]->Fill(x, y, w4 * p45lab2);
   }
   return kTRUE;
}

void hpolar::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 hpolar::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("hpolar.root", "recreate");
   TIter next(GetOutputList());
   while (TObject *obj = next()) {
      TH2D *h = (TH2D*)obj;
      h->Write();
   }
}