gUpMuFluxGen.cxx
Go to the documentation of this file.
1 //________________________________________________________________________________________
2 /*!
3 
4 \program gevgen_upmu
5 
6 \brief A GENIE upgoing muon event generation application.
7 
8  *** Synopsis :
9 
10  gevgen_upmu [-h]
11  [-r run#]
12  -f flux
13  -n n_of_events
14  -d detector_bounding_box_size
15  -g rock_composition
16  [--seed random_number_seed]
17  --cross-sections xml_file
18  [--message-thresholds xml_file]
19 
20  *** Options :
21 
22  [] Denotes an optional argument.
23 
24  -h
25  Prints out the syntax and exits.
26  -r
27  Specifies the MC run number
28  [default: 1000]
29  -f
30  Specifies the input flux files
31  The general syntax is: `-f simulation:/path/file.data[neutrino_code],...'
32  [Notes]
33  - The `simulation' string can be either `FLUKA' or `BGLRS' (so that
34  input data are binned using the correct FLUKA and BGLRS energy and
35  costheta binning). See comments in
36  - $GENIE/src/Flux/GFLUKAAtmoFlux.h
37  - $GENIE/src/Flux/GBGLRSAtmoFlux.h
38  and follow the links to the FLUKA and BGLRS atmo. flux web pages.
39  - The neutrino codes are the PDG ones, for numu and anumu only
40  - The /path/file.data,neutrino_code part of the option can be
41  repeated multiple times (separated by commas), once for each
42  flux neutrino species you want to consider,
43  eg. '-f FLUKA:~/data/sdave_numu07.dat[14],~/data/sdave_nue07.dat[12]'
44  eg. '-f BGLRS:~/data/flux10_271003_z.kam_nue[12]'
45  -n
46  Specifies how many events to generate.
47  -d
48  Specifies side length (in mm) of the detector bounding box.
49  [default 100m (100000mm)]
50  -g
51  rock composition << NOT IMPLEMENTED YET >>
52  Rock composition should be implemented in terms of materials defined in
53  src/MuELoss/MuELMaterial.h and their weight fraction
54  eg like -g 'silicon[0.30],calcium[0.29],iron[0.02],...'
55  --seed
56  Random number seed.
57  --cross-sections
58  Name (incl. full path) of an XML file with pre-computed
59  cross-section values used for constructing splines.
60  --message-thresholds
61  Allows users to customize the message stream thresholds.
62  The thresholds are specified using an XML file.
63  See $GENIE/config/Messenger.xml for the XML schema.
64 
65  *** Examples:
66 
67  (1) Generate 100k events (run number 999210) for nu_mu only, using the
68  sdave_numu07.dat FLUKA flux file (files in /data/flx/), and a detector of
69  side length 50m.
70 
71  % gevgen_upmu -r 999210 -n 100000 -d 50000
72  -f FLUKA:/data/flx/sdave_numu07.dat[14]
73 
74 
75  You can further control the GENIE behaviour by setting its standard
76  environmental variables.
77  Please read the GENIE User Manual for more information.
78 
79 \created April 15, 2011
80 
81 \author Jen Truby <jen.truby \at sky.com>
82  Oxford University - University of Liverpool & STFC Rutherford Appleton Lab summer student
83 
84  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
85  University of Liverpool & STFC Rutherford Appleton Lab
86 
87 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
88  For the full text of the license visit http://copyright.genie-mc.org
89  or see $GENIE/LICENSE
90 */
91 //_________________________________________________________________________________________
92 
93 #include <cassert>
94 #include <cstdlib>
95 #include <cctype>
96 #include <string>
97 #include <vector>
98 #include <sstream>
99 #include <map>
100 
101 #include <TRotation.h>
102 #include <TH1D.h>
103 #include <TH3D.h>
104 #include <TTree.h>
105 #include <TLorentzVector.h>
106 
125 #include "Framework/Utils/AppInit.h"
126 #include "Framework/Utils/RunOpt.h"
128 
129 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
132 #endif
133 
134 using std::string;
135 using std::vector;
136 using std::map;
137 using std::ostringstream;
138 
139 using namespace genie;
140 using namespace genie::flux;
141 using namespace genie::constants;
142 
143 void GetCommandLineArgs (int argc, char ** argv);
144 void PrintSyntax (void);
145 GFluxI * GetFlux (void);
146 void GenerateUpNu (GFluxI * flux_driver);
147 TH3D * BuildEmuEnuCosThetaPdf (int nu_code);
148 TH1D * GetEmuPdf (double Enu, double costheta, const TH3D* pdf3d);
149 TVector3 GetDetectorVertex (double CosTheta, double Enu);
150 double GetCrossSection (int nu_code, double Enu, double Emu);
151 double ProbabilityEmu (int nu_code, double Enu, double Emu);
152 
153 // User-specified options:
154 //
155 Long_t gOptRunNu; // run number
156 string gOptFluxSim; // flux simulation (FLUKA or BGLRS)
157 map<int,string> gOptFluxFiles; // neutrino pdg code -> flux file map
158 int gOptNev = -1; // exposure - in terms of number of events
159 double gOptDetectorSide; // detector side length, in mm.
160 long int gOptRanSeed; // random number seed
161 string gOptInpXSecFile; // cross-section splines
162 
163 // Constants
164 const double a = 2e+6; // a = 2 MeV / (g cm-2)
165 const double e = 500e+9; // e = 500 GeV
166 
167 // Defaults:
168 //
169 double kDefOptDetectorSide = 1e+5; // side length of detector, 100m (in mm)
170 
171 //________________________________________________________________________________________
172 int main(int argc, char** argv)
173 {
174  // Parse command line arguments
175  GetCommandLineArgs(argc,argv);
176 
177  if ( ! RunOpt::Instance()->Tune() ) {
178  LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
179  exit(-1);
180  }
182 
183  // Init
184  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
187 
188 
189  // Get requested flux driver
190  GFluxI * flux_driver = GetFlux();
191 
192  // Create output tree to store generated up-going muons
193  TTree * ntupmuflux = new TTree("ntupmuflux","GENIE Upgoing Muon Event Tree");
194  // Tree branches
195  int brIev = 0; // Event number
196  int brNuCode = 0; // Neutrino PDG code
197  double brEmu = 0; // Muon energy (GeV)
198  double brEnu = 0; // Neutrino energy (GeV)
199  double brCosTheta = 0; // Neutrino cos(zenith angle), muon assumed to be collinear
200  double brWghtFlxNu = 0; // Weight associated with the current flux neutrino
201  double brWghtEmuPdf = 0; // Weight associated with the Emu pdf for current Enu, costheta bins
202  double brVx = 0; // Muon x (mm) - intersection with box surrounding the detector volume
203  double brVy = 0; // Muon y (mm) - intersection with box surrounding the detector volume
204  double brVz = 0; // Muon z (mm) - intersection with box surrounding the detector volume
205  double brXSec = 0; //
206  ntupmuflux->Branch("iev", &brIev, "iev/I" );
207  ntupmuflux->Branch("nu_code", &brNuCode, "nu_code/I" );
208  ntupmuflux->Branch("Emu", &brEmu, "Emu/D" );
209  ntupmuflux->Branch("Enu", &brEnu, "Enu/D" );
210  ntupmuflux->Branch("costheta", &brCosTheta, "costheta/D" );
211  ntupmuflux->Branch("wght_fluxnu", &brWghtFlxNu, "wght_fluxnu/D" );
212  ntupmuflux->Branch("wght_emupdf", &brWghtEmuPdf, "wght_emupdf/D" );
213  ntupmuflux->Branch("vx", &brVx, "vx/D" );
214  ntupmuflux->Branch("vy", &brVy, "vy/D" );
215  ntupmuflux->Branch("vz", &brVz, "vz/D" );
216  ntupmuflux->Branch("xsec", &brXSec, "xsec/D" );
217 
218  // Build 3-D pdfs describing the the probability of a muon neutrino (or anti-neutrino)
219  // of energy Enu and zenith angle costheta producing a mu- (or mu+) of energy E_mu
220  TH3D * pdf3d_numu = BuildEmuEnuCosThetaPdf (kPdgNuMu );
221  TH3D * pdf3d_numubar = BuildEmuEnuCosThetaPdf (kPdgAntiNuMu);
222 
223  // Up-going muon event loop
224  for(brIev = 0; brIev < gOptNev; brIev++) {
225 
226  // Loop until upgoing nu is generated.
227  GenerateUpNu(flux_driver);
228 
229  // Get neutrino code, Enu, costheta and weight for the generated neutrino
230  brNuCode = flux_driver->PdgCode();
231  brEnu = flux_driver->Momentum().E();
232  brCosTheta = -1. * flux_driver->Momentum().Pz() / flux_driver->Momentum().Vect().Mag();
233  brWghtFlxNu = flux_driver->Weight();
234 
235  LOG("gevgen_upmu", pNOTICE)
236  << "Generated flux neutrino: code = " << brNuCode
237  << ", Ev = " << brEnu << " GeV"
238  << ", cos(theta) = " << brCosTheta
239  << ", weight = " << brWghtFlxNu;
240 
241  // Get Emu pdf my slicing the 3-D Enu,Emu,costheta pdf.
242  TH3D * pdf3d = (brNuCode == kPdgNuMu) ? pdf3d_numu : pdf3d_numubar;
243  TH1D * pdfEmu = GetEmuPdf(brEnu,brCosTheta,pdf3d);
244 
245  // Get a random Emu from the pdf, and get the weight for that Emu.
246  brEmu = pdfEmu->GetRandom();
247  brWghtEmuPdf = pdfEmu->Integral("width");
248  LOG("gevgen_upmu", pNOTICE)
249  << "Selected muon has energy Emu = " << brEmu
250  << " and Emu pdg weight = " << brWghtEmuPdf;
251 
252  // Randomly select the point at which the neutrino crosses the detector.
253  // assumes a simple cube detector of side length gOptDetectorSide.
254  TVector3 Vertex = GetDetectorVertex(brCosTheta,brEnu);
255  brVx = Vertex.X();
256  brVy = Vertex.Y();
257  brVz = Vertex.Z();
258 
259  LOG("gevgen_upmu", pNOTICE)
260  << "Generated muon position: (" << brVx << ", " << brVy << ", " << brVz << ") m";
261 
262  // Save the cross section for the interaction generated.
263  brXSec = GetCrossSection(brNuCode, brEnu, brEmu);
264 
265  // save all relevant values to the ntuple
266  ntupmuflux->Fill();
267 
268  // Clean-up
269  delete pdfEmu;
270  }
271 
272  // Save the muon ntuple and calculate 3-D pdfs
273 
274  ostringstream outfname;
275  outfname << "./genie." << gOptRunNu << ".upmu.root";
276  TFile outf(outfname.str().c_str(), "recreate");
277  ntupmuflux -> Write("ntupmu");
278  pdf3d_numu -> Write("pdf3d_numu");
279  pdf3d_numubar -> Write("pdf3d_numubar");
280  outf.Close();
281 
282  // Clean-up
283  delete flux_driver;
284 
285  return 0;
286 }
287 //________________________________________________________________________________________
288 TVector3 GetDetectorVertex(double costheta, double Enu)
289 {
290 // Find the point at which the muon crosses the detector. Returns 0,0,0 if the muon misses.
291 // Detector is a cube of side length l (=gOptDetectorSide).
292 
293  // Get a RandomGen instance
295 
296  // Generate random phi.
297  double phi = 2.*kPi* rnd->RndFlux().Rndm();
298 
299  // Set distance at which incoming muon is considered
300  double rad = 0.87*gOptDetectorSide;
301 
302  // Set transverse radius of a circle
303  double rad_trans = 0.87*gOptDetectorSide;
304 
305  // Get necessary trig
306  double sintheta = TMath::Sqrt(1-costheta*costheta);
307  double cosphi = TMath::Cos(phi);
308  double sinphi = TMath::Sin(phi);
309 
310  // Compute the muon position at distance Rad.
311  double z = rad * costheta;
312  double y = rad * sintheta * cosphi;
313  double x = rad * sintheta * sinphi;
314 
315  // Displace muon randomly on a circular surface of radius rad_trans,
316  // perpendicular to a sphere radius rad at that position [x,y,z].
317  TVector3 vec(x,y,z); // vector towards selected point
318  TVector3 dvec1 = vec.Orthogonal(); // orthogonal vector
319  TVector3 dvec2 = dvec1; // second orthogonal vector
320  dvec2.Rotate(-kPi/2.0,vec); // rotate second vector by 90deg -> Cartesian coords
321 
322  // Select a random point on the transverse surface, within radius rad_trans
323  double psi = 2 * kPi * rnd->RndFlux().Rndm(); // rndm angle [0,2pi]
324  double random = rnd->RndFlux().Rndm(); // rndm number [0,1]
325  dvec1.SetMag(TMath::Sqrt(random)*rad_trans*TMath::Cos(psi));
326  dvec2.SetMag(TMath::Sqrt(random)*rad_trans*TMath::Sin(psi));
327  x += dvec1.X() + dvec2.X(); // x-coord of a point the muon passes through
328  y += dvec1.Y() + dvec2.Y(); // y-coord of a point the muon passes through
329  z += dvec1.Z() + dvec2.Z(); // z-coord of a point the muon passes through
330 
331  // Find out if the muon passes through any side of the detector.
332 
333  // Get momentum vector
334  double pz = Enu * costheta;
335  double py = Enu * sintheta * sinphi;
336  double px = Enu * sintheta * cosphi;
337  TVector3 p3(-px,-py,-pz);
338 
339  // Set up other vectors needed for line-box intersection.
340  TVector3 x3(x,y,z);
341  TVector3 temp3(x3);
342  TVector3 Hit3(0,0,0);
343 
344  // Find out if the line of the muon intersects the detector.
345  bool HitFound = false;
346  double l = gOptDetectorSide;
347  TVector3 unit(0,0,0);
348  unit = p3.Unit();
349  double unitx = unit.X();
350  double unity = unit.Y();
351  double unitz = unit.Z();
352 
353  // Check to see if muon intersects with z=-l/2 surface.
354  if (x3.Z() < -l/2 && unitz > 0 && !HitFound){
355  while (temp3.Z() < -l/2){
356  temp3 += unit;
357  }
358  if ( -l/2<temp3.X() && l/2>temp3.X() && -l/2<temp3.Y() && l/2>temp3.Y() ){
359  Hit3.SetXYZ(temp3.X(),temp3.Y(),-l/2);
360  HitFound = true;
361  }
362  else temp3 = x3;
363  }
364 
365  // Check to see if muon intersects with x=l/2 surface.
366  if (x3.X() > l/2 && unitx < 0 && !HitFound){
367  while (temp3.X() > l/2){
368  temp3 += p3.Unit();
369  }
370  if ( -l/2<temp3.Y() && l/2>temp3.Y() && -l/2<temp3.Z() && l/2>temp3.Z() ){
371  Hit3.SetXYZ(l/2,temp3.Y(),temp3.Z());
372  HitFound = true;
373  }
374  else temp3 = x3;
375  }
376 
377  // Check to see if muon intersects with y=l/2 surface.
378  if (x3.Y() > l/2 && unity < 0 && !HitFound){
379  while (temp3.Y() > l/2){
380  temp3 += p3.Unit();
381  }
382  if ( -l/2<temp3.X() && l/2>temp3.X() && -l/2<temp3.Z() && l/2>temp3.Z() ){
383  Hit3.SetXYZ(temp3.X(),l/2,temp3.Z());
384  HitFound = true;
385  }
386  else temp3 = x3;
387  }
388 
389  // Check to see if muon intersects with x=-l/2 surface.
390  if (x3.X() < -l/2 && unitx > 0 && !HitFound){
391  while (temp3.X() < -l/2){
392  temp3 += p3.Unit();
393  }
394  if ( -l/2<temp3.Y() && l/2>temp3.Y() && -l/2<temp3.Z() && l/2>temp3.Z() ){
395  Hit3.SetXYZ(-l/2,temp3.Y(),temp3.Z());
396  HitFound = true;
397  }
398  else temp3 = x3;
399  }
400 
401  // Check to see if muon intersects with y=-l/2 surface.
402  if (x3.Y() < -l/2 && unity > 0 && !HitFound){
403  while (temp3.Y() < -l/2){
404  temp3 += p3.Unit();
405  }
406  if ( -l/2<temp3.X() && l/2>temp3.X() && -l/2<temp3.Z() && l/2>temp3.Z() ){
407  Hit3.SetXYZ(temp3.X(),-l/2,temp3.Z());
408  HitFound = true;
409  }
410  else temp3 = x3;
411  }
412 
413  return Hit3;
414 }
415 //________________________________________________________________________________________
416 void GenerateUpNu(GFluxI * flux_driver)
417 {
418 // Generate a Neutrino. Keep generating neutrinos until an upgoing one is generated.
419 
420  while (1)
421  {
422  LOG("gevgen_upmu", pINFO) << "Pulling next neurtino from the flux driver...";
423  flux_driver->GenerateNext();
424 
425  int nu_code = flux_driver->PdgCode();
426 
427  const TLorentzVector & p4 = flux_driver->Momentum();
428  double costheta = -p4.Pz() / p4.Vect().Mag();
429 
430  bool keep = (costheta<0) && (nu_code==kPdgNuMu || nu_code==kPdgAntiNuMu);
431  if(keep) return;
432  }
433 }
434 //________________________________________________________________________________________
435 double ProbabilityEmu(int nu_code, double Enu, double Emu)
436 {
437 // Calculate the probability of an incoming neutrino of energy Enu
438 // generating a muon of energy Emu.
439  double dxsec_dxdy = GetCrossSection(nu_code,Enu,Emu);
440  double Int = e * constants::kNA * dxsec_dxdy / (a * (1 + Emu/e));
441  return Int;
442 }
443 //________________________________________________________________________________________
444 TH3D* BuildEmuEnuCosThetaPdf(int nu_code)
445 {
446 // Set up a 3D histogram, with axes Emu, Enu, CosTheta.
447 // Bin convention is defined at the start.
448 // Content of each bin is given by the probability of getting a muon
449 // of energy Emu from a neutrino of energy Enu and zenith ange costheta
450 
451  // Bin convention for Enu, consistent with BGLRS
452  const double Enumin = 0.1;
453  const int nEnubinsPerDecade = 10;
454  const int nEnubins = 31;
455 
456  // Bin convention for Emu, consistent with BGLRS
457  const double Emumin = 0.1;
458  const int nEmubinsPerDecade = 10;
459  const int nEmubins = 31;
460 
461  // Bin convention for CosTheta, consistent with BGLRS
462  const double costheta_min = -1;
463  const double costheta_max = +1;
464  const int ncostheta = 20;
465 
466  // Set up an array of CosTheta bins.
467  double costhetabinwidth = ((costheta_max-costheta_min)/ncostheta);
468  double CosThetaBinEdges[ncostheta+1];
469  for (int i=0; i<=ncostheta; i++)
470  {
471  CosThetaBinEdges[i] = costheta_min + i*costhetabinwidth;
472  }
473 
474  // Set up an array of Enu bins.
475  Double_t MinLogEnu = log(Enumin);
476  Double_t MaxLogEnu = log(10*Enumin);
477  Double_t LogBinWidthEnu = ((MaxLogEnu-MinLogEnu)/nEnubinsPerDecade);
478  Double_t EnuBinEdges[nEnubins+1];
479  for (int i=0; i<=nEnubins; i++)
480  {
481  EnuBinEdges[i] = exp(MinLogEnu + i*LogBinWidthEnu);
482  }
483 
484  // Set up an array of Emu bins.
485  Double_t MinLogEmu = log(Emumin);
486  Double_t MaxLogEmu = log(10*Emumin);
487  Double_t LogBinWidthEmu = ((MaxLogEmu-MinLogEmu)/nEmubinsPerDecade);
488  Double_t EmuBinEdges[nEmubins+2];
489  for (int i=0; i<=nEmubins+1; i++)
490  {
491  EmuBinEdges[i] = exp(MinLogEmu + i*LogBinWidthEmu);
492  }
493 
494  // Create 3D histogram. X-axis: Emu; Y-axis: Enu; Z-axis: CosTheta.
495  TH3D *h3 = new TH3D("h3","",nEmubins+1,EmuBinEdges,nEnubins,EnuBinEdges,ncostheta,CosThetaBinEdges);
496 
497  // Draw histogram.
498  //h3->Draw();
499 
500  // Calculate Probability for each triplet.
501  for (int i=1; i< h3->GetNbinsX(); i++)
502  {
503  double Emu = h3->GetXaxis()->GetBinCenter(i);
504  for (int j=1; j<= h3->GetNbinsY(); j++)
505  {
506  double Enu = h3->GetBinCenter(j);
507  double Int = ProbabilityEmu(nu_code,Enu,Emu);
508  for (int k=1; k<= h3->GetNbinsZ(); k++)
509  {
510  h3->SetBinContent(i,j,k,Int); // Bin Content will is Int.
511  //h3->SetBinContent(i,j,k,1); // Bin content constant (to test)
512  }
513  }
514  }
515  return h3;
516 }
517 //________________________________________________________________________________________
518 TH1D * GetEmuPdf(double Enu, double costheta, const TH3D* pdf3d)
519 {
520  // Build 1-D Emu pdf
521  int Emu_nbins = pdf3d->GetXaxis()->GetNbins();
522  const double * Emu_bins = pdf3d->GetXaxis()->GetXbins()->GetArray();
523  TH1D * pdf_Emu = new TH1D("pdf_Emu","",Emu_nbins,Emu_bins);
524  pdf_Emu->SetDirectory(0);
525 
526  // Find appropriate bins for Enu and CosTheta.
527  int Enu_bin = pdf3d->GetYaxis()->FindBin(Enu);
528  int costheta_bin = pdf3d->GetZaxis()->FindBin(costheta);
529 
530  // For those Enu and costheta bins, build Emu pdf
531  for (int Emu_bin = 1; Emu_bin < pdf_Emu->GetNbinsX(); Emu_bin++) {
532  pdf_Emu->SetBinContent(Emu_bin, pdf3d->GetBinContent(Emu_bin,Enu_bin,costheta_bin));
533  }
534 
535  return pdf_Emu;
536 }
537 //________________________________________________________________________________________
538 double GetCrossSection(int nu_code, double Enu, double Emu)
539 {
540 // Get the interaction cross section for a neutrino of energy Enu
541 // to produce a muon of energy Emu.
542 
543  double dxsec_dxdy = 0;
544 
545  if ( Emu >= Enu ) {
546  dxsec_dxdy = 0;
547  }
548  else {
549  // get the algorithm factory & config pool
550  AlgFactory * algf = AlgFactory::Instance();
551 
552  // instantiate some algorithms
553  AlgId id("genie::QPMDISPXSec","Default");
554  Algorithm * alg = algf->AdoptAlgorithm(id);
555  XSecAlgorithmI * xsecalg = dynamic_cast<XSecAlgorithmI*>(alg);
556 
557  Interaction * vp = Interaction::DISCC(kPdgTgtFreeP, kPdgProton, nu_code, Enu);
559 
560  // Integrate over x [0,1] and y [0,1-Emu/Enu], 100 steps for each.
561  double ymax = 1 - Emu/Enu;
562 
563  double dy = ymax/10;
564  double dx = 0.1;
565  double x = 0;
566  double y = 0;
567 
568  for (int i = 0; i<=10; i++){
569  x = i*dx;
570  for (int j = 0; j<=10; j++){
571  y = j * dy;
572  double W = 0;
573  double Q2 = 0;
575  vp->KinePtr()->Setx(x);
576  vp->KinePtr()->Sety(y);
577  vp->KinePtr()->SetQ2(Q2);
578  vp->KinePtr()->SetW(W);
579  vn->KinePtr()->Setx(x);
580  vn->KinePtr()->Sety(y);
581  vn->KinePtr()->SetQ2(Q2);
582  vn->KinePtr()->SetW(W);
583 
584  dxsec_dxdy += 0.5*(xsecalg->XSec(vp,kPSxyfE) + xsecalg->XSec(vn,kPSxyfE)) / units::cm2;
585  }
586  }
587 
588  delete vp;
589  delete vn;
590  }
591 
592  return dxsec_dxdy;
593 }
594 //________________________________________________________________________________________
596 {
597  GFluxI * flux_driver = 0;
598 
599 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
600 
601  // Instantiate appropriate concrete flux driver
602  GAtmoFlux * atmo_flux_driver = 0;
603  if(gOptFluxSim == "FLUKA") {
604  GFLUKAAtmoFlux * fluka_flux = new GFLUKAAtmoFlux;
605  atmo_flux_driver = dynamic_cast<GAtmoFlux *>(fluka_flux);
606  } else
607  if(gOptFluxSim == "BGLRS") {
608  GBGLRSAtmoFlux * bartol_flux = new GBGLRSAtmoFlux;
609  atmo_flux_driver = dynamic_cast<GAtmoFlux *>(bartol_flux);
610  } else {
611  LOG("gevgen_upmu", pFATAL) << "Uknonwn flux simulation: " << gOptFluxSim;
612  gAbortingInErr = true;
613  exit(1);
614  }
615 
616 
617  atmo_flux_driver->GenerateWeighted(true);
618 
619  // Configure GAtmoFlux options (common to all concrete atmospheric flux drivers)
620  // set flux files:
621  map<int,string>::const_iterator file_iter = gOptFluxFiles.begin();
622  for( ; file_iter != gOptFluxFiles.end(); ++file_iter) {
623  int neutrino_code = file_iter->first;
624  string filename = file_iter->second;
625  atmo_flux_driver->AddFluxFile(neutrino_code, filename);
626  }
627  atmo_flux_driver->LoadFluxData();
628  // configure flux generation surface:
629  atmo_flux_driver->SetRadii(1, 1);
630  // Cast to GFluxI, the generic flux driver interface
631  flux_driver = dynamic_cast<GFluxI *>(atmo_flux_driver);
632 
633 #else
634  LOG("gevgen_upmu", pFATAL) << "You need to enable the GENIE flux drivers first!";
635  LOG("gevgen_upmu", pFATAL) << "Use --enable-flux-drivers at the configuration step.";
636  gAbortingInErr = true;
637  exit(1);
638 #endif
639 
640  return flux_driver;
641 }
642 //________________________________________________________________________________________
643 void GetCommandLineArgs(int argc, char ** argv)
644 {
645 // Get the command line arguments
646 
647  LOG("gevgen_upmu", pINFO) << "Parsing command line arguments";
648 
649  // Common run options. Set defaults and read.
651 
652  // Parse run options for this app
653 
654  CmdLnArgParser parser(argc,argv);
655 
656  // help?
657  bool help = parser.OptionExists('h');
658  if(help) {
659  PrintSyntax();
660  exit(0);
661  }
662 
663  //
664  // *** run number
665  //
666  if( parser.OptionExists('r') ) {
667  LOG("gevgen_upmu", pDEBUG) << "Reading MC run number";
668  gOptRunNu = parser.ArgAsLong('r');
669  } else {
670  LOG("gevgen_upmu", pDEBUG) << "Unspecified run number - Using default";
671  gOptRunNu = 1000;
672  } //-r
673 
674  //
675  // *** exposure
676  //
677 
678  // in number of events
679  bool have_required_statistics = false;
680  if( parser.OptionExists('n') ) {
681  LOG("gevgen_upmu", pDEBUG)
682  << "Reading number of events to generate";
683  gOptNev = parser.ArgAsInt('n');
684  have_required_statistics = true;
685  }//-n
686  if(!have_required_statistics) {
687  LOG("gevgen_upmu", pFATAL)
688  << "You must request exposure either in terms of number of events"
689  << "\nUse the -n option";
690  PrintSyntax();
691  gAbortingInErr = true;
692  exit(1);
693  }
694 
695  //
696  // *** detector side length
697  //
698 
699  if( parser.OptionExists('d') ) {
700  LOG("gevgen_upmu", pDEBUG)
701  << "Reading detector side length";
702  gOptDetectorSide = parser.ArgAsDouble('d');
703  } else {
704  LOG("gevgen_upmu", pDEBUG)
705  << "Unspecified detector side length - Using default";
707  }//-d
708 
709  //
710  // *** flux files
711  //
712 
713  // syntax:
714  // simulation:/path/file.data[neutrino_code],/path/file.data[neutrino_code],...
715  //
716  if( parser.OptionExists('f') ) {
717  LOG("gevgen_upmu", pDEBUG) << "Getting input flux files";
718  string flux = parser.ArgAsString('f');
719 
720  // get flux simulation info (FLUKA,BGLRS) so as to instantiate the
721  // appropriate flux driver
722  string::size_type jsimend = flux.find_first_of(":",0);
723  if(jsimend==string::npos) {
724  LOG("gevgen_upmu", pFATAL)
725  << "You need to specify the flux file source";
726  PrintSyntax();
727  gAbortingInErr = true;
728  exit(1);
729  }
730  gOptFluxSim = flux.substr(0,jsimend);
731  for(string::size_type i=0; i<gOptFluxSim.size(); i++) {
732  gOptFluxSim[i] = toupper(gOptFluxSim[i]);
733  }
734  if((gOptFluxSim != "FLUKA") && (gOptFluxSim != "BGLRS")) {
735  LOG("gevgen_upmu", pFATAL)
736  << "The flux file source needs to be one of <FLUKA,BGLRS>";
737  PrintSyntax();
738  gAbortingInErr = true;
739  exit(1);
740  }
741  // now get the list of input files and the corresponding neutrino codes.
742  flux.erase(0,jsimend+1);
743  vector<string> fluxv = utils::str::Split(flux,",");
744  vector<string>::const_iterator fluxiter = fluxv.begin();
745  for( ; fluxiter != fluxv.end(); ++fluxiter) {
746  string filename_and_pdg = *fluxiter;
747  string::size_type open_bracket = filename_and_pdg.find("[");
748  string::size_type close_bracket = filename_and_pdg.find("]");
749  if (open_bracket ==string::npos ||
750  close_bracket==string::npos)
751  {
752  LOG("gevgen_upmu", pFATAL)
753  << "You made an error in specifying the flux info";
754  PrintSyntax();
755  gAbortingInErr = true;
756  exit(1);
757  }
758  string::size_type ibeg = 0;
759  string::size_type iend = open_bracket;
760  string::size_type jbeg = open_bracket+1;
761  string::size_type jend = close_bracket;
762  string flux_filename = filename_and_pdg.substr(ibeg,iend-ibeg);
763  string neutrino_pdg = filename_and_pdg.substr(jbeg,jend-jbeg);
764  gOptFluxFiles.insert(
765  map<int,string>::value_type(atoi(neutrino_pdg.c_str()), flux_filename));
766  }
767  if(gOptFluxFiles.size() == 0) {
768  LOG("gevgen_upmu", pFATAL)
769  << "You must specify at least one flux file!";
770  PrintSyntax();
771  gAbortingInErr = true;
772  exit(1);
773  }
774 
775  } else {
776  LOG("gevgen_upmu", pFATAL) << "No flux info was specified! Use the -f option.";
777  PrintSyntax();
778  gAbortingInErr = true;
779  exit(1);
780  }
781 
782  //
783  // *** random number seed
784  //
785  if( parser.OptionExists("seed") ) {
786  LOG("gevgen_upmu", pINFO) << "Reading random number seed";
787  gOptRanSeed = parser.ArgAsLong("seed");
788  } else {
789  LOG("gevgen_upmu", pINFO) << "Unspecified random number seed - Using default";
790  gOptRanSeed = -1;
791  }
792 
793  //
794  // *** input cross-section file
795  //
796  if( parser.OptionExists("cross-sections") ) {
797  LOG("gevgen_upmu", pINFO) << "Reading cross-section file";
798  gOptInpXSecFile = parser.ArgAsString("cross-sections");
799  } else {
800  LOG("gevgen_upmu", pINFO) << "Unspecified cross-section file";
801  gOptInpXSecFile = "";
802  }
803 
804 
805  //
806  // print-out summary
807  //
808 
809  PDGLibrary * pdglib = PDGLibrary::Instance();
810 
811  ostringstream fluxinfo;
812  fluxinfo << "Using " << gOptFluxSim << " flux files: ";
813  map<int,string>::const_iterator file_iter = gOptFluxFiles.begin();
814  for( ; file_iter != gOptFluxFiles.end(); ++file_iter) {
815  int neutrino_code = file_iter->first;
816  string filename = file_iter->second;
817  TParticlePDG * p = pdglib->Find(neutrino_code);
818  if(p) {
819  string name = p->GetName();
820  fluxinfo << "(" << name << ") -> " << filename << " / ";
821  }
822  }
823 
824  ostringstream expinfo;
825  if(gOptNev > 0) { expinfo << gOptNev << " events"; }
826 
827  LOG("gevgen_atmo", pNOTICE)
828  << "\n\n"
829  << utils::print::PrintFramedMesg("gevgen_upmu job configuration");
830 
831  LOG("gevgen_upmu", pNOTICE)
832  << "\n"
833  << "\n @@ Run number: " << gOptRunNu
834  << "\n @@ Random number seed: " << gOptRanSeed
835  << "\n @@ Using cross-section file: " << gOptInpXSecFile
836  << "\n @@ Flux"
837  << "\n\t" << fluxinfo.str()
838  << "\n @@ Exposure"
839  << "\n\t" << expinfo.str()
840  << "\n\n";
841 }
842 //________________________________________________________________________________________
843 void PrintSyntax(void)
844 {
845  LOG("gevgen_upmu", pFATAL)
846  << "\n **Syntax**"
847  << "\n gevgen_upmu [-h]"
848  << "\n [-r run#]"
849  << "\n -f simulation:flux_file[neutrino_code],..."
850  << "\n -n n_of_events,"
851  << "\n [-d detector side length (mm)]"
852  << "\n [--seed random_number_seed]"
853  << "\n --cross-sections xml_file"
854  << "\n [--message-thresholds xml_file]"
855  << "\n"
856  << " Please also read the detailed documentation at http://www.genie-mc.org"
857  << "\n";
858 }
859 //________________________________________________________________________________________
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
Cross Section Calculation Interface.
void RandGen(long int seed)
Definition: AppInit.cxx:31
const double kPi
double ProbabilityEmu(int nu_code, double Enu, double Emu)
const XML_Char * name
Definition: expat.h:151
void XSecTable(string inpfile, bool require_table)
Definition: AppInit.cxx:39
long ArgAsLong(char opt)
Basic constants.
Long_t gOptRunNu
double ArgAsDouble(char opt)
string gOptInpXSecFile
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
static const double kNucleonMass
Definition: Constants.h:78
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
void GenerateUpNu(GFluxI *flux_driver)
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:265
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
int gOptNev
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
TH1F * h3
Definition: berger.C:36
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:107
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
double GetCrossSection(int nu_code, double Enu, double Emu)
const int kPdgNuMu
Definition: PDGCodes.h:30
double kDefOptDetectorSide
Algorithm abstract base class.
Definition: Algorithm.h:54
string filename
Definition: shutoffs.py:106
void SetRadii(double Rlongitudinal, double Rtransverse)
Definition: GAtmoFlux.cxx:412
void GetCommandLineArgs(int argc, char **argv)
static const double cm2
Definition: Units.h:77
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
virtual void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
Definition: GAtmoFlux.cxx:289
Loaders::FluxType flux
int main(int argc, char **argv)
TH1D * GetEmuPdf(double Enu, double costheta, const TH3D *pdf3d)
long int gOptRanSeed
Double_t ymax
Definition: plot.C:25
GENIE flux drivers.
Summary information for an interaction.
Definition: Interaction.h:56
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Definition: typedefs.hpp:11
double dy[NP][NC]
const double e
double dx[NP][NC]
A flux driver for the FLUKA 3-D Atmospheric Neutrino Flux.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
GFluxI * GetFlux(void)
const double a
void XYtoWQ2(double Ev, double M, double &W, double &Q2, double x, double y)
Definition: KineUtils.cxx:1069
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
const int kPdgTgtFreeN
Definition: PDGCodes.h:177
const int kPdgTgtFreeP
Definition: PDGCodes.h:176
static const double kNA
Definition: Constants.h:50
TFile * outf
Definition: testXsec.C:51
const double j
Definition: BetheBloch.cxx:29
#define pINFO
Definition: Messenger.h:63
Eigen::VectorXd vec
string gOptFluxSim
Algorithm * AdoptAlgorithm(const AlgId &algid) const
Definition: AlgFactory.cxx:127
static constexpr Double_t rad
Definition: Munits.h:162
z
Definition: test.py:28
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
void BuildTune()
build tune and inform XSecSplineList
Definition: RunOpt.cxx:100
virtual bool GenerateNext(void)=0
generate the next flux neutrino (return false in err)
void PrintSyntax(void)
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:241
Algorithm ID (algorithm name + configuration set name)
Definition: AlgId.h:35
TH3D * BuildEmuEnuCosThetaPdf(int nu_code)
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
static RunOpt * Instance(void)
Definition: RunOpt.cxx:62
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:253
static Interaction * DISCC(int tgt, int nuc, int probe, double E=0)
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
void AddFluxFile(int neutrino_pdg, string filename)
Definition: GAtmoFlux.cxx:421
virtual int PdgCode(void)=0
returns the flux neutrino pdg code
exit(0)
virtual const TLorentzVector & Momentum(void)=0
returns the flux neutrino 4-momentum
virtual double Weight(void)=0
returns the flux neutrino weight (if any)
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
Definition: PrintUtils.cxx:171
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:100
const int kPdgProton
Definition: PDGCodes.h:65
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:62
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:72
A base class for the FLUKA, BGLRS and ATMNC atmo. nu. flux drivers. The driver depends on data files ...
Definition: GAtmoFlux.h:61
bool gAbortingInErr
Definition: Messenger.cxx:56
#define W(x)
A flux driver for the Bartol Atmospheric Neutrino Flux.
const int kPdgNeutron
Definition: PDGCodes.h:67
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
bool OptionExists(char opt)
was option set?
map< int, string > gOptFluxFiles
TVector3 GetDetectorVertex(double CosTheta, double Enu)
#define pDEBUG
Definition: Messenger.h:64
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
double gOptDetectorSide
gm Write()
enum BeamMode string