gtestGiBUUData.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program testGiBUUData
5 
6 \brief Program used for testing and debugging the input GiBUU data tables
7  and interpolation.
8 
9 \syntax gtestGiBUUData
10 
11 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
12  University of Liverpool & STFC Rutherford Appleton Lab
13 
14 \created May 31, 2009
15 
16 \cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
17  For the full text of the license visit http://copyright.genie-mc.org
18  or see $GENIE/LICENSE
19 */
20 //____________________________________________________________________________
21 
22 #include <TFile.h>
23 #include <TTree.h>
24 
25 #include "BaryonResonance/BaryonResUtils.h"
26 #include "GiBUU/GiBUUData.h"
27 #include "Messenger/Messenger.h"
28 #include "Numerical/Spline.h"
29 #include "PDG/PDGCodes.h"
30 
31 using namespace genie;
32 using namespace genie::utils;
33 
34 void SaveToRootFile(void);
35 
36 int main(int /*argc*/, char ** /*argv*/)
37 {
39  return 0;
40 }
41 
42 void SaveToRootFile(void)
43 {
44  //
45  // get GiBUU data
46  //
47 
48  GiBUUData * gd = GiBUUData::Instance();
49 
50  const GiBUUData::FormFactors & ff = gd->FF();
51 
52  //
53  // define inputs
54  //
55 
56  const int kNRes = 13;
57  const int kNNuc = 2;
58  const int kNIntTyp = 3;
59  const int kNQ2 = 200;
60 
61  Resonance_t resonance[kNRes] = {
65 
66  int nucleon[kNNuc] = {
68 
69  InteractionType_t interaction[kNIntTyp] = {
71 
72  double Q2min = 0;
73  double Q2max = 4;
74  double dQ2 = (Q2max-Q2min)/(kNQ2-1);
75 
76  //
77  // define the output tree
78  //
79 
80  int bResId; // resonance id
81  int bNucleon; // nucleon pdg
82  int bIntType; // interaction type (cc, nc, em)
83  int bIsoX2; // resonance isospin x 2
84  double bQ2; // Q2
85  double bC3V; // C3V f/f
86  double bC4V; // C4V f/f
87  double bC5V; // C5V f/f
88  double bC6V; // C6V f/f
89  double bC3A; // C3A f/f
90  double bC4A; // C4A f/f
91  double bC5A; // C5A f/f
92  double bC6A; // C6A f/f
93  double bF1V; // F1V f/f
94  double bF2V; // F2V f/f
95  double bFA; // FA f/f
96  double bFP; // FP f/f
97 
98  TTree * gibuu_ffres = new TTree("gibuu_ffres","GiBUU resonance form factors");
99 
100  gibuu_ffres -> Branch ("res", &bResId, "res/I" );
101  gibuu_ffres -> Branch ("nuc", &bNucleon, "nuc/I" );
102  gibuu_ffres -> Branch ("intyp", &bIntType, "intyp/I" );
103  gibuu_ffres -> Branch ("isoX2", &bIsoX2, "isoX2/I" );
104  gibuu_ffres -> Branch ("Q2", &bQ2, "Q2/D" );
105  gibuu_ffres -> Branch ("C3V", &bC3V, "C3V/D" );
106  gibuu_ffres -> Branch ("C4V", &bC4V, "C4V/D" );
107  gibuu_ffres -> Branch ("C5V", &bC5V, "C5V/D" );
108  gibuu_ffres -> Branch ("C6V", &bC6V, "C6V/D" );
109  gibuu_ffres -> Branch ("C3A", &bC3A, "C3A/D" );
110  gibuu_ffres -> Branch ("C4A", &bC4A, "C4A/D" );
111  gibuu_ffres -> Branch ("C5A", &bC5A, "C5A/D" );
112  gibuu_ffres -> Branch ("C6A", &bC6A, "C6A/D" );
113  gibuu_ffres -> Branch ("F1V", &bF1V, "F1V/D" );
114  gibuu_ffres -> Branch ("F2V", &bF2V, "F2V/D" );
115  gibuu_ffres -> Branch ("FA", &bFA, "FA/D" );
116  gibuu_ffres -> Branch ("FP", &bFP, "FP/D" );
117 
118  //
119  // calculate(interpolate) resonance f/f and fill-in the tree
120  //
121 
122  for(int r=0; r<kNRes; r++) {
123  for(int i=0; i<kNNuc; i++) {
124  for(int j=0; j<kNIntTyp; j++) {
125 
126  bResId = (int) resonance[r];
127  bNucleon = nucleon[i];
128  bIntType = interaction[j];
129  bIsoX2 = res::IsDelta(resonance[r]) ? 3 : 1;
130 
131  for(int iq2=0; iq2<kNQ2; iq2++) {
132  double Q2 = Q2min + iq2 * dQ2;
133 
134  bQ2 = Q2;
135  bC3V = ff.C3V(Q2, resonance[r], nucleon[i], interaction[j]);
136  bC4V = ff.C4V(Q2, resonance[r], nucleon[i], interaction[j]);
137  bC5V = ff.C5V(Q2, resonance[r], nucleon[i], interaction[j]);
138  bC6V = ff.C6V(Q2, resonance[r], nucleon[i], interaction[j]);
139  bC3A = ff.C3A(Q2, resonance[r], nucleon[i], interaction[j]);
140  bC4A = ff.C4A(Q2, resonance[r], nucleon[i], interaction[j]);
141  bC5A = ff.C5A(Q2, resonance[r], nucleon[i], interaction[j]);
142  bC6A = ff.C6A(Q2, resonance[r], nucleon[i], interaction[j]);
143  bF1V = ff.F1V(Q2, resonance[r], nucleon[i], interaction[j]);
144  bF2V = ff.F2V(Q2, resonance[r], nucleon[i], interaction[j]);
145  bFA = ff.FA (Q2, resonance[r], nucleon[i], interaction[j]);
146  bFP = ff.FP (Q2, resonance[r], nucleon[i], interaction[j]);
147 
148  gibuu_ffres->Fill();
149  }//Q2
150  }//interaction
151  }//nucleon
152  }//res
153 
154  //
155  // save
156  //
157  TFile f("./gibuu_data.root", "recreate");
158  gibuu_ffres->Write();
159  f.Close();
160 }
int main(int, char **)
bool IsDelta(Resonance_t res)
is it a Delta resonance?
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
enum genie::EResonance Resonance_t
const double j
Definition: BetheBloch.cxx:29
TRandom3 r(0)
const int kPdgProton
Definition: PDGCodes.h:65
TFile ff[ntarg]
Definition: Style.C:26
const int kNRes
Definition: bwnorm.C:23
const int kPdgNeutron
Definition: PDGCodes.h:67
enum genie::EInteractionType InteractionType_t
Root of GENIE utility namespaces.
void SaveToRootFile(void)