es.cpp

/**
 * Plot ES scattering angle for 8B solar neutrinos vs. recoil energy,
 * weighted by the neutrino spectrum.
 *
 * This uses the 8B spectrum and ES cross section code in RAT.
 *
 * Build:
 *
 *     g++ -g -o es es.cpp -L$RATROOT/lib -I$RATROOT/include \
 *     -lrat_Linux `root-config --libs --cflags` \
 *     `geant4-config --libs --cflags` -lcurl -lbz2
 *
 * Run:
 *
 *     $ ./es
 */

#include <iostream>
#include <RAT/DB.hh>
#include <RAT/RunManager.hh>
#include <RAT/ESCrossSec.hh>
#include <RAT/ESgen.hh>
#include <TGraph.h>
#include <TCanvas.h>
#include <TAxis.h>
#include <G4RunManager.hh>
#include <G4UImanager.hh>
#include <CLHEP/Units/PhysicalConstants.h>

#define DEBUG false  // Enable/disable detailed output

int main(int argc, char* argv[]) {
  // Minimal initialization of RAT
  RAT::DB::Get()->LoadDefaults();
  RAT::gTheRunManager = new RAT::RunManager;
  G4UImanager::GetUIpointer()->ApplyCommand("/PhysicsList/OmitHadronicProcesses true");
  G4RunManager::GetRunManager()->Initialize();
  float me = electron_mass_c2;

  // Get solar 8B spectrum
  RAT::DBLinkPtr lsolar = RAT::DB::Get()->GetLink("SOLAR", "b8");
  std::vector<double> x = lsolar->GetDArray("spec_e");
  std::vector<double> y = lsolar->GetDArray("spec_flux");
  double total = 0;
  for (size_t i=0; i<y.size(); i++) {
    total += y[i];
  }
  for (size_t i=0; i<y.size(); i++) {
    y[i] /= total;
  }

  RAT::ESCrossSec xs("nue");

  TCanvas c;
  TGraph g;

  for (size_t i=0; i<x.size(); i++) {
    double ee = x[i] * MeV;
    double thetasum = 0;
    double count = 0;
    for (size_t j=i; j<x.size(); j++) {
      double enu = x[j] * MeV;
      double t = y[j] * (x[1]-x[0]) * acos(sqrt((ee*(me + enu)*(me + enu)) / (2*me*enu*enu + enu*enu*ee)));
      if (t == t) {
        thetasum += t;
        count += y[j] * (x[1]-x[0]);
      }
      if (DEBUG) {
        std::cout << ee << " " << enu << " " << y[j] << " " << me << " " << t << std::endl;
      }
    }

    if (count > 0) {
      thetasum /= count;
    }
    else {
      thetasum = 0;
    }

    g.SetPoint(i, ee, thetasum * 180 / 3.14159);
    std::cout << i << ": " << ee << " / " << thetasum << std::endl;
  }

  // Plot
  g.Draw("al");
  g.SetName("esangles");
  g.GetXaxis()->SetTitle("#beta kinetic energy (MeV)");
  g.GetYaxis()->SetTitle("#theta");
  c.Update();
  c.SaveAs("es_angles.root");
  c.SaveAs("es_angles.C");
  c.SaveAs("es_angles.pdf");

  return 0;
}