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


#include "phonopolar.h"
#include <TStyle.h>
#include <TList.h>
#include <sstream>

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

   std::stringstream name;
   std::stringstream title;
   name << "hcolspectrum";
   title << "collimated photon energy spectrum";
   hcolspectrum = new TH1D(name.str().c_str(), title.str().c_str(),
                           100, 0., 13.);
   hcolspectrum->GetXaxis()->SetTitle("E_#gamma (GeV)");
   hcolspectrum->SetDirectory(0);
   GetOutputList()->Add(hcolspectrum);

   name.str("");
   title.str("");
   name << "hcolpolarization";
   title << "collimated photon linear polarization 0-90";
   hcolpolarization = new TProfile(name.str().c_str(), title.str().c_str(),
                                   100, 0., 13., -1, 1);
   hcolpolarization->GetXaxis()->SetTitle("E_#gamma (GeV)");
   hcolpolarization->GetYaxis()->SetTitle("linear polarization");
   hcolpolarization->SetDirectory(0);
   GetOutputList()->Add(hcolpolarization);

   name.str("");
   title.str("");
   name << "hkhist";
   title << "h vs k";
   hkhist = new TH2D(name.str().c_str(), title.str().c_str(),
                     220, -110, 110, 220, -110, 110);
   hkhist->GetXaxis()->SetTitle("h");
   hkhist->GetYaxis()->SetTitle("k");
   hkhist->SetDirectory(0);
   GetOutputList()->Add(hkhist);

   name.str("");
   title.str("");
   name << "lhist";
   title << "l";
   lhist = new TH1D(name.str().c_str(), title.str().c_str(),
                    30, -15, 15);
   lhist->GetXaxis()->SetTitle("l");
   lhist->SetDirectory(0);
   GetOutputList()->Add(lhist);

   for (int i=0; i < 9; ++i) {
      name.str("");
      title.str("");
      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 << "h1polar_" << i;
      title << "beam spot polarization map (0-90), "
            << 3+i << "-" << 4+i << "GeV";
      h1polar[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                            50, -10., 10., 50, -10., 10.);
      name.str("");
      title.str("");
      name << "h2polar_" << i;
      title << "beam spot w**2 p**2 map (0-90), "
            << 3+i << "-" << 4+i << "GeV";
      h2polar[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                            50, -10., 10., 50, -10., 10.);
      name.str("");
      title.str("");
      name << "h3polar_" << i;
      title << "beam spot w**3 p**2 map (0-90), "
            << 3+i << "-" << 4+i << "GeV";
      h3polar[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                            50, -10., 10., 50, -10., 10.);
      name.str("");
      title.str("");
      name << "h4polar_" << i;
      title << "beam spot w**4 p**2 map (0-90), "
            << 3+i << "-" << 4+i << "GeV";
      h4polar[i] = new TH2D(name.str().c_str(), title.str().c_str(),
                            50, -10., 10., 50, -10., 10.);
      hintens[i]->SetDirectory(0);
      h1polar[i]->SetDirectory(0);
      h2polar[i]->SetDirectory(0);
      h3polar[i]->SetDirectory(0);
      h4polar[i]->SetDirectory(0);
      GetOutputList()->Add(hintens[i]);
      GetOutputList()->Add(h1polar[i]);
      GetOutputList()->Add(h2polar[i]);
      GetOutputList()->Add(h3polar[i]);
      GetOutputList()->Add(h4polar[i]);
   }
}

Bool_t phonopolar::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[0] / k[2];
   double y = distanceRadCol * k[1] / k[2];
   double phi = atan2(y,x);
   double w = *diffXS * *wgt;
   double p0lab = *polar * 2 - 1;
   int i = int(k[2] - 3.);
   if (i >= 0 && i < 9) {
      double w2 = w * w;
      double w3 = w2 * w;
      double w4 = w3 * w;
      double p0lab2 = p0lab * p0lab;
      hintens[i]->Fill(x, y, w);
      h1polar[i]->Fill(x, y, w * p0lab);
      h2polar[i]->Fill(x, y, w2 * p0lab2);
      h3polar[i]->Fill(x, y, w3 * p0lab2);
      h4polar[i]->Fill(x, y, w4 * p0lab2);
   }
   if (sqrt(x*x + y*y) < 2.5) {
      hcolspectrum->Fill(k[2], w);
      hcolpolarization->Fill(k[2], p0lab, w);
      hkhist->Fill(hkl[0], hkl[1], w);
      lhist->Fill(hkl[2], w);
   }
   return kTRUE;
}

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