testXsec.C
Go to the documentation of this file.
1 // shell% genie
2 // genie[0] .x test_xsec.C
3 //
4 
5 #if !defined(__CINT__) || defined(__MAKECINT__)
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TNtuple.h"
9 
11 #include "EVGCore/InteractionList.h"
12 #include "VLE/StrumiaVissaniIBDPXSec.h"
13 #include "PDG/PDGCodeList.h"
14 #include "PDG/PDGUtils.h"
15 #include "Messenger/Messenger.h"
16 
17 #endif
18 
19 #include "PDG/PDGCodes.h"
20 #include "Conventions/Units.h"
23 
24 const Int_t nknots = 44;
25 const Double_t E[nknots] = {
26  /*1.806,*/ 2.01, 2.25, 2.51, 2.80,
27  3.12, 3.48, 3.89, 4.33, 4.84,
28  5.40, 6.02, 6.72, 7.49, 8.36,
29  8.83, 9.85, 11.0, 12.3, 13.7,
30  15.3, 17.0, 19.0, 21.2, 23.6,
31  26.4, 29.4, 32.8, 36.6, 40.9,
32  43.2, 48.2, 53.7, 59.9, 66.9,
33  74.6, 83.2, 92.9, 104, 116,
34  129, 144, 160, 179, 200
35 };
36 const Double_t s[nknots] = {
37  /*0,*/ 0.00351, 0.00735, 0.0127, 0.0202,
38  0.0304, 0.0440, 0.06190, 0.0854, 0.116,
39  0.155, 0.205, 0.269, 0.349, 0.451,
40  0.511, 0.654, 0.832, 1.05, 1.33,
41  1.67, 2.09, 2.61, 3.24, 4.01,
42  4.95, 6.08, 7.44, 9.08, 11.0,
43  12.1, 14.7, 17.6, 21.0, 25.0,
44  29.6, 34.8, 40.7, 47.3, 54.6,
45  62.7, 71.5, 81.0, 91.3, 102.
46 };
47 Double_t xsec[nknots], xsec_n[nknots];
48 
49 using namespace genie;
50 
51 TFile* outf;
52 TNtuple* nt;
53 
54 void testXsec(const Char_t* outfn="VLExsecNT.root") {
55  gROOT->Macro("loadlibs.C");
56 
57 /*
58  PDGCodeList * targets = new PDGCodeList;
59  targets->push_back(kPdgTgtFreeP);
60 
61  PDGCodeList * neutrinos = new PDGCodeList;
62  neutrinos->push_back(kPdgAntiNuE);
63 */
64 
65  // load VLE xsec parameters
67  const Registry * config =
68  pool->FindRegistry("genie::StrumiaVissaniIBDPXSec/Default");
69 
70  // our xsec calculator
72  alg->Configure(*config);
73 
74  // prepare to save results
75  outf = TFile::Open(outfn,"recreate");
76  nt = new TNtuple("nt","nt","Ev:xsec:xsnun:xspaper");
77  nt->SetDirectory(outf);
78 
79 
80  // nu_e_bar + p --> n + e+
83  genie::Interaction * interaction =
84  new genie::Interaction(init_state, proc_info);
85  Target * target = interaction->InitStatePtr()->TgtPtr();
86  target->SetHitNucPdg(kPdgProton);
87  Int_t i=0;
88  for (i = 0; i < nknots; i++) {
89  TLorentzVector p4(0,0,E[i]*1e-3,E[i]*1e-3);
90  interaction->InitStatePtr()->SetProbeP4(p4);
91  xsec[i] = alg->Integral(interaction);
92  }
93 
94 
95  // nu_e + n --> p + e-
96  InitialState init_state_n(kPdgTgtFreeN, kPdgNuE);
98  genie::Interaction * interaction_n =
99  new genie::Interaction(init_state_n, proc_info_n);
100  Target * target_n = interaction_n->InitStatePtr()->TgtPtr();
101  target_n->SetHitNucPdg(kPdgNeutron);
102  for (i = 0; i < nknots; i++) {
103  TLorentzVector n4(0,0,E[i]*1e-3,E[i]*1e-3);
104  interaction_n->InitStatePtr()->SetProbeP4(n4);
105  xsec_n[i] = alg->Integral(interaction_n);
106  }
107 
108  for (i = 0; i < nknots; i++) {
109  Printf("xsec(E=%g MeV) = %g \t xsec_n = %g", E[i],
110  (1E+41/units::cm2)*xsec[i],
111  (1E+41/units::cm2)*xsec_n[i]);
112  nt->Fill(E[i],
113  (1E+41/units::cm2)*xsec[i],
114  (1E+41/units::cm2)*xsec_n[i],
115  s[i]);
116  }
117 
118  outf->Write();
119 
120 }
Cross Section Calculation Interface.
const int kPdgNuE
Definition: PDGCodes.h:28
const XML_Char * target
Definition: expat.h:268
void SetProbeP4(const TLorentzVector &P4)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
const int kPdgAntiNuE
Definition: PDGCodes.h:29
A singleton class holding all configuration registries built while parsing all loaded XML configurati...
Definition: AlgConfigPool.h:41
Definition: config.py:1
An implementation of the neutrino - (free) nucleon [inverse beta decay] cross section, valid from the threshold energy (1.806MeV) up to hundreds of MeV. Currently cut off at 1/2 nucleon mass. Based on the Strumia/Vissani paper Phys.Lett.B564:42-54,2003.
static const double cm2
Definition: Units.h:77
const XML_Char * s
Definition: expat.h:262
Summary information for an interaction.
Definition: Interaction.h:56
const Double_t E[nknots]
Definition: testXsec.C:25
void testXsec(const Char_t *outfn="VLExsecNT.root")
Definition: testXsec.C:54
const Int_t nknots
Definition: testXsec.C:24
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
virtual double Integral(const Interaction *i) const =0
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
const int kPdgTgtFreeN
Definition: PDGCodes.h:177
const int kPdgTgtFreeP
Definition: PDGCodes.h:176
TFile * outf
Definition: testXsec.C:51
Double_t xsec[nknots]
Definition: testXsec.C:47
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:188
Target * TgtPtr(void) const
Definition: InitialState.h:68
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
Registry * FindRegistry(string key) const
Double_t xsec_n[nknots]
Definition: testXsec.C:47
const int kPdgProton
Definition: PDGCodes.h:65
Float_t e
Definition: plot.C:35
TNtuple * nt
Definition: testXsec.C:52
const int kPdgNeutron
Definition: PDGCodes.h:67
static AlgConfigPool * Instance()
Initial State information.
Definition: InitialState.h:49