GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gtestXSecCEvNS.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestXSecCEvNS
5 
6 \brief test program used for testing the CEvNS cross-section calculation
7 
8 \author Costas Andreopoulos <c.andreopoulos \at cern.ch>
9  University of Liverpool
10 
11 \created July 15, 2019
12 
13 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
14  For the full text of the license visit http://copyright.genie-mc.org
15 */
16 //____________________________________________________________________________
17 
18 #include <TFile.h>
19 #include <TGraph.h>
20 
27 #include "Framework/Utils/RunOpt.h"
28 
29 using namespace genie;
30 //__________________________________________________________________________
31 int main(int argc, char ** argv)
32 {
33  // <temp/>
35  if ( ! RunOpt::Instance()->Tune() ) {
36  LOG("test", pFATAL) << " No TuneId in RunOption";
37  exit(-1);
38  }
40  // </temp>
41 
42  LOG("test",pINFO) << "Testing the CEvNS cross-section calculation";
43 
44  // Instantiate the coherent elastic cross-section model
46  AlgId id("genie::PattonCEvNSPXSec","Default");
47  const Algorithm * alg = algf->GetAlgorithm(id);
48  LOG("test", pINFO) << *alg;
49  const XSecAlgorithmI * xsecalg = dynamic_cast<const XSecAlgorithmI*>(alg);
50 
51  const int target = 1000180400;
52  const int probe = 12;
53  double E = 0.010; // GeV
54  double Q2 = 0.0001; //GeV^2
55 
56  Interaction * interaction = Interaction::CEvNS(target, probe, E);
57  interaction->KinePtr()->SetQ2(Q2);
58 
59 /*
60  // test differential cross-section calculation at a single point
61  double dxsec = xsecalg->XSec(interaction,kPSQ2fE);
62 
63  LOG("test", pNOTICE)
64  << "dxsec[CEvNS; target = " << target << "]/dQ2"
65  << "(E = " << E << " GeV, Q2 = " << Q2 << " GeV^2) = "
66  << dxsec/(units::cm2) << " cm^2/GeV^2";
67 */
68 
69 /*
70  // test integrated cross-section calculation at a single point
71  double xsec = xsecalg->Integral(interaction);
72 
73  LOG("test", pNOTICE)
74  << "xsec[CEvNS; target = " << target << "]"
75  << "(E = " << E << " GeV) = " << xsec/(units::cm2) << " cm^2";
76 */
77 
78  // test integrated cross-section calculation at a range of energies
79  // and compare with published predictions
80  const int n = 100;
81  const double Emin = 0.005;
82  const double Emax = 0.055;
83  const double dE = (Emax-Emin)/(n-1);
84 
85  double e_array[n] = {0};
86  double genie_xsec_array[n] = {0};
87 
88  for(int i=0; i<n; i++) {
89  double E_current = Emin + i*dE;
90  interaction->InitStatePtr()->SetProbeE(E_current);
91  double xsec_current = xsecalg->Integral(interaction)/(units::cm2);
92  e_array[i] = E_current;
93  genie_xsec_array[i] = xsec_current;
94  LOG("test", pNOTICE)
95  << "xsec[CEvNS; target = " << target << "]"
96  << "(E = " << E_current << " GeV) = " << xsec_current << " cm^2";
97  }
98  TGraph * genie_xsec = new TGraph(n, e_array, genie_xsec_array);
99  TGraph * published_xsec = new TGraph("$GENIE/src/contrib/cevns/cevns_arXiv180309183_Ar40_prediction.data");
100  published_xsec->SetMarkerStyle(20);
101  published_xsec->SetMarkerColor(kRed);
102  genie_xsec->SetLineColor(kBlue);
103  TFile f("cevns.root","recreate");
104  genie_xsec->Write("genie_xsec_nuAr40");
105  published_xsec->Write("published_xsec_nuAr40");
106  f.Close();
107  return 0;
108 }
109 //__________________________________________________________________________
Cross Section Calculation Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:99
#define pFATAL
Definition: Messenger.h:56
static Interaction * CEvNS(int tgt, int probe, double E=0)
Algorithm abstract base class.
Definition: Algorithm.h:54
void SetProbeE(double E)
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double cm2
Definition: Units.h:69
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
#define pINFO
Definition: Messenger.h:62
void BuildTune()
build tune and inform XSecSplineList
Definition: RunOpt.cxx:92
Algorithm ID (algorithm name + configuration set name)
Definition: AlgId.h:34
static RunOpt * Instance(void)
Definition: RunOpt.cxx:54
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
#define pNOTICE
Definition: Messenger.h:61
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
virtual double Integral(const Interaction *i) const =0