JDIS.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <map>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TRandom3.h"

#include "JPhysics/JDIS.hh"

#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * Example application to display DIS of muon.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string         outputFile;
  int            numberOfEvents;
  UInt_t         seed;
  int            debug;

  try { 

    JParser<> zap("Example application to display DIS of muon.");

    zap['o'] = make_field(outputFile)       = "dis.root";
    zap['n'] = make_field(numberOfEvents)   = 0;
    zap['S'] = make_field(seed)             = 0;
    zap['d'] = make_field(debug)            = 0;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  gRandom->SetSeed(seed);
 

  const JDIS dis;


  TH1D h0("sigma", NULL, 1000, -1.0, 9.0);

  JManager<double, TH1D> HA(new TH1D("v  [% GeV]", NULL, 1000, -5.0, 0.0));
  JManager<double, TH1D> H1(new TH1D("ran[% GeV]", NULL, 1000, -5.0, 0.0));
  JManager<double, TH1D> HB(new TH1D("RMS[% GeV]", NULL, 1000, -5.0, 0.0));

  H1->Sumw2();

  const map<double, double> setups = { { 1.0e2, 0.7 },
                                       { 1.0e3, 0.6 },
                                       { 1.0e4, 0.5 } };


  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
      
    const double x = h0.GetBinCenter(i);
    const double E = pow(10.0, x);
    const double y = dis.getCrossSection(E);
    
    h0.SetBinContent(i, y * 1.0e+30);  // ub
  }
  
  for (const auto& setup : setups) {

    const double E = setup.first;

    TH1D* ha = HA[E];
    TH1D* h1 = H1[E];
    TH1D* hb = HB[E];

    for (int i = 1; i <= ha->GetNbinsX(); ++i) {
      
      const double x = ha->GetBinCenter(i);
      const double v = pow(10.0, x);

      double y = dis.getP(E, v);
      
      if (numberOfEvents == 0) {
        y *= setup.second / dis.getP(E, dis.getV(E));
      }

      ha->SetBinContent(i, y);
    }

    if (numberOfEvents != 0) {

      for (int i = 0; i != numberOfEvents; ++i) {

        const double Es = dis.getE(E);
        const double v  = Es/E;

        h1->Fill(log10(v));
      }

      for (int i = 1; i <= h1->GetNbinsX(); ++i) {

        const double y =  h1->GetBinContent(i);
        const double z =  h1->GetBinError  (i);

        const double xmin = h1->GetXaxis()->GetBinLowEdge(i);
        const double xmax = h1->GetXaxis()->GetBinUpEdge (i);

        const double w = (pow(10.0,xmax) - pow(10.0,xmin)) * numberOfEvents;

        h1->SetBinContent(i, y / w);
        h1->SetBinError  (i, z / w);
      }
    }

    for (int i = 1; i <= hb->GetNbinsX(); ++i) {
      
      const double x = hb->GetBinCenter(i);
      const double v = pow(10.0, x);
      const double y = dis.getThetaRMS(E, v);

      hb->SetBinContent(i, y);
    }
  }


  TFile out(outputFile.c_str(), "recreate");

  out << h0 << HA << H1 << HB;

  out.Write();
  out.Close();
}