gtestFluxAstro.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestFluxAstro
5 
6 \brief Test program for astrophysical neutrino flux drivers
7 
8 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
9  University of Liverpool & STFC Rutherford Appleton Lab
10 
11 \created March 26, 2010
12 
13 \cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
14  For the full text of the license visit http://copyright.genie-mc.org
15  or see $GENIE/LICENSE
16 */
17 //____________________________________________________________________________
18 
19 #include <TFile.h>
20 #include <TNtuple.h>
21 
22 #include "Conventions/Constants.h"
23 #include "Conventions/Units.h"
24 #include "FluxDrivers/GAstroFlux.h"
25 #include "Messenger/Messenger.h"
26 #include "PDG/PDGCodes.h"
27 
28 using namespace genie;
29 using namespace genie::flux;
30 
31 const unsigned int kNEvents = 1000000;
32 
33 //____________________________________________________________________________
34 int main(int /*argc*/, char ** /*argv*/)
35 {
36  const double pi = constants::kPi;
37 
38  const double latitude = pi/5;
39  const double longitude = pi/4;
40  const double depth = 3.0*units::km;
41  const double size = 2.0*units::km;
42 
43  GDiffuseAstroFlux * difflx = new GDiffuseAstroFlux;
44 
45  difflx->ForceMinEnergy(1E+2);
46  difflx->ForceMaxEnergy(1E+8);
47  difflx->GenerateWeighted(true);
48  difflx->SetDetectorPosition(latitude,longitude,depth,size);
49  difflx->SetRelNuPopulations(1,2,0,1,2,0);
50  difflx->SetEnergyPowLawIdx(3.5);
51 
52  TNtuple * fluxntp =
53  new TNtuple("fluxntp", "flux", "x:y:z:t:px:py:pz:E:pdgc:wght");
54 
55  unsigned int ievent = 0;
56  while(ievent++ < kNEvents) {
57  LOG("test", pINFO) << "Event number: " << ievent;
58  difflx->GenerateNext();
59  int pdgc = difflx->PdgCode();
60  double wght = difflx->Weight();
61  const TLorentzVector & x4 = difflx->Position();
62  const TLorentzVector & p4 = difflx->Momentum();
63  fluxntp->Fill(
64  x4.X(), x4.Y(), x4.Z(), x4.T(),
65  p4.Px(), p4.Py(), p4.Pz(), p4.E(), pdgc, wght);
66  }
67 
68  TFile f("./genie-astro-flux.root","recreate");
69  fluxntp->Write();
70  f.Close();
71 
72  LOG("test", pINFO) << "Done!";
73 
74  return 0;
75 }
76 //____________________________________________________________________________
virtual const TLorentzVector & Momentum(void)
returns the flux neutrino 4-momentum
Definition: GAstroFlux.h:142
void ForceMinEnergy(double emin)
Definition: GAstroFlux.cxx:150
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
void SetEnergyPowLawIdx(double n)
Definition: GAstroFlux.cxx:270
void ForceMaxEnergy(double emax)
Definition: GAstroFlux.cxx:156
virtual bool GenerateNext(void)
generate the next flux neutrino (return false in err)
Definition: GAstroFlux.cxx:59
GENIE flux drivers.
static const double km
Definition: Units.h:72
int main(int, char **)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
virtual double Weight(void)
returns the flux neutrino weight (if any)
Definition: GAstroFlux.h:141
Float_t E
Definition: plot.C:20
const unsigned int kNEvents
#define pINFO
Definition: Messenger.h:63
void SetDetectorPosition(double latitude, double longitude, double depth, double size)
Definition: GAstroFlux.cxx:174
virtual const TLorentzVector & Position(void)
returns the flux neutrino 4-position (note: expect SI rather than physical units) ...
Definition: GAstroFlux.h:143
virtual void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
Definition: GAstroFlux.cxx:169
static const double kPi
Definition: Constants.h:38
virtual int PdgCode(void)
returns the flux neutrino pdg code
Definition: GAstroFlux.h:140
void SetRelNuPopulations(double nnue=1, double nnumu=2, double nnutau=0, double nnuebar=1, double nnumubar=2, double nnutaubar=0)
Definition: GAstroFlux.cxx:220