GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gtestFluxSimple.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestFluxSimple
5 
6 \brief Test program for simple flux drivers
7 
8 \author Costas Andreopoulos <c.andreopoulos \at cern.ch>
9  University of Liverpool
10 
11 \created August 22, 2005
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 
19 #include <TFile.h>
20 #include <TNtuple.h>
21 #include <TSystem.h>
22 #include <TH1D.h>
23 #include <TF1.h>
24 
28 
29 using namespace genie;
30 using namespace genie::flux;
31 
32 const unsigned int kNEvents = 10000;
33 
34 TNtuple * runGCylindTH1FluxDriver (void);
35 TNtuple * createFluxNtuple (GFluxI * flux);
36 
37 //___________________________________________________________________
38 int main(int /*argc*/, char ** /*argv*/)
39 {
40  LOG("test", pINFO) << "Running GCylindTH1Flux driver test";
41  TNtuple * ntcylh1f = runGCylindTH1FluxDriver();
42  ntcylh1f->SetTitle("GCylindTH1Flux driver data");
43 
44  LOG("test", pINFO) << "Saving flux ntuples";
45 
46  TFile f("./genie-flux-drivers.root","recreate");
47  ntcylh1f -> Write("ntcylh1f");
48  f.Close();
49 
50  delete ntcylh1f;
51 
52  LOG("test", pINFO) << "Done!";
53 
54  return 0;
55 }
56 //___________________________________________________________________
57 TNtuple * runGCylindTH1FluxDriver(void)
58 {
59  LOG("test", pINFO) << "Creating GCylindTH1Flux flux driver";
60 
61  GCylindTH1Flux * flux = new GCylindTH1Flux;
62 
63  LOG("test", pINFO) << "Setting configuration data";
64 
65  TF1 * f1 = new TF1("f1","1./x",0.5,5.0);
66  TH1D * spectrum1 = new TH1D("spectrum1","numu E", 20,0.5,5);
67  spectrum1->FillRandom("f1",100000);
68 
69  TF1 * f2 = new TF1("f2","x",0.5,5.0);
70  TH1D * spectrum2 = new TH1D("spectrum2","numubar E", 20,0.5,5);
71  spectrum2->FillRandom("f2",10000);
72 
73  TVector3 direction(0,0,1);
74  TVector3 beam_spot(0,0,-10);
75 
76  double Rtransverse = 0.5;
77 
78  LOG("test", pINFO) << "Configuring GCylindTH1Flux flux driver";
79 
80  flux -> SetNuDirection (direction);
81  flux -> SetBeamSpot (beam_spot);
82  flux -> SetTransverseRadius (Rtransverse);
83  flux -> AddEnergySpectrum (kPdgNuMu, spectrum1);
84  flux -> AddEnergySpectrum (kPdgAntiNuMu, spectrum2);
85 
86  LOG("test", pINFO) << "Creating flux ntuple";
87  GFluxI * fluxi = dynamic_cast<GFluxI*>(flux);
88 
89  TNtuple * fluxntp = createFluxNtuple(fluxi);
90 
91  delete f1;
92  delete f2;
93  delete flux;
94 
95  return fluxntp;
96 }
97 //___________________________________________________________________
98 TNtuple * createFluxNtuple(GFluxI * flux)
99 {
100  TNtuple * fluxntp =
101  new TNtuple("fluxntp",
102  "flux data", "x:y:z:t:px:py:pz:E:pdgc");
103 
104  LOG("test", pINFO) << "Generating flux neutrinos";
105 
106  unsigned int ievent = 0;
107  while(ievent++ < kNEvents) {
108 
109  flux->GenerateNext();
110 
111  int pdgc = flux->PdgCode();
112  const TLorentzVector & x4 = flux->Position();
113  const TLorentzVector & p4 = flux->Momentum();
114 
115  fluxntp->Fill( x4.X(), x4.Y(), x4.Z(), x4.T(),
116  p4.Px(), p4.Py(), p4.Pz(), p4.E(), pdgc);
117  }
118  return fluxntp;
119 }
120 //___________________________________________________________________
121 
TNtuple * runGCylindTH1FluxDriver(void)
virtual const TLorentzVector & Position(void)=0
returns the flux neutrino 4-position (note: expect SI rather than physical units) ...
const int kPdgNuMu
Definition: PDGCodes.h:30
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
A generic GENIE flux driver. Generates a &#39;cylindrical&#39; neutrino beam along the input direction...
const unsigned int kNEvents
#define pINFO
Definition: Messenger.h:62
TNtuple * createFluxNtuple(GFluxI *flux)
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
virtual bool GenerateNext(void)=0
generate the next flux neutrino (return false in err)
virtual int PdgCode(void)=0
returns the flux neutrino pdg code
virtual const TLorentzVector & Momentum(void)=0
returns the flux neutrino 4-momentum
hadnt Write("hadnt")
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:29