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;
  ULong_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();
}