GAtmoFlux.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 05, 2008 - CA
14  This class was added in 2.3.1 by code factored out from the concrete
15  GFLUKAAtmoFlux driver & the newer, largely similar, GBGLRSAtmoFlux
16  driver.
17  @ Feb 23, 2010 - CA
18  Re-structuring and clean-up. Added option to generate weighted flux.
19  Added option to specify a maximum energy cut.
20  @ Feb 24, 2010 - CA
21  Added option to specify a minimum energy cut.
22  @ Sep 22, 2010 - TF, CA
23  Added SetUserCoordSystem(TRotation &) to specify a rotation from the
24  Topocentric Horizontal (THZ) coordinate system to a user-defined
25  topocentric coordinate system. Added NFluxNeutrinos() to get number of
26  flux neutrinos generated for sample normalization purposes (note that, in
27  the presence of cuts, this is not the same as the number of flux neutrinos
28  thrown towards the geometry).
29  @ Feb 22, 2011 - JD
30  Implemented dummy versions of the new GFluxI::Clear and GFluxI::Index as
31  these methods needed for pre-generation of flux interaction probabilities
32  in GMCJDriver.
33  @ Feb 03, 2011 - TF
34  Bug fixed: events are now generated randomly and uniformly on a disc with
35  R = R_{transverse}
36  @ Feb 23, 2012 - AB
37  Bug fixed: events were being generated according to the differential flux
38  in each energy bin, dPhi/dE, rather than the total flux, Phi, in each bin.
39  This has now been fixed.
40  @ Jul 13, 2015 - CA
41  Changed all internal flux histograms from 2-D to 3-D to accomodate fluxes
42  (HAKKM) with phi dependence. Changed several names accordingly. For FLUKA
43  and BGLRS fluxes, just add a single phi bin from 0 to 2*pi.
44  @ Apr 13, 2016 - CA
45  Added AddFluxFile() without neutrino PDG code argument for the case that
46  the input data file (as it is the case for HAKKM) includes all species.
47  The pure virtual method FillFluxHisto() now takes a neutrino PDG code
48  rather than a TH3D ptr input and it is expected to retrieve the TH3D flux
49  itself. This change was made to easily fit HAKKM in the code already used
50  by FLUKA and BGLRS.
51 
52 */
53 //____________________________________________________________________________
54 
55 #include <cassert>
56 #include <iostream>
57 #include <fstream>
58 
59 #include <TH3D.h>
60 #include <TMath.h>
61 
63 #include "Tools/Flux/GAtmoFlux.h"
71 
72 using namespace genie;
73 using namespace genie::flux;
74 using namespace genie::constants;
75 
76 //____________________________________________________________________________
78 {
79  fInitialized = 0;
80 }
81 //___________________________________________________________________________
83 {
84  this->CleanUp();
85 }
86 //___________________________________________________________________________
88 {
89  return TMath::Min(fMaxEv, fMaxEvCut);
90 }
91 //___________________________________________________________________________
93 {
94  while(1) {
95  // Attempt to generate next flux neutrino
96  bool nextok = this->GenerateNext_1try();
97  if(!nextok) continue;
98 
99  // Check generated neutrino energy against max energy.
100  // We may have to reject the current neutrino if a user-defined max
101  // energy cut restricts the available range of energies.
102  const TLorentzVector & p4 = this->Momentum();
103  double E = p4.Energy();
104  double Emin = this->MinEnergy();
105  double Emax = this->MaxEnergy();
106  double wght = this->Weight();
107 
108  bool accept = (E<=Emax && E>=Emin && wght>0);
109  if(accept) return true;
110  }
111  return false;
112 }
113 //___________________________________________________________________________
115 {
116  // Must have run intitialization
118 
119  // Reset previously generated neutrino code / 4-p / 4-x
120  this->ResetSelection();
121 
122  // Get a RandomGen instance
124 
125  // Generate (Ev, costheta, phi)
126  double Ev = 0.;
127  double costheta = 0.;
128  double phi = 0;
129  double weight = 0;
130  int nu_pdg = 0;
131 
132  if(fGenWeighted) {
133 
134  //
135  // generate weighted flux
136  //
137 
138  // generate events according to a power law spectrum,
139  // then weight events by flux and inverse power law
140  // (note: cannot use index alpha=1)
141  double alpha = fSpectralIndex;
142 
143  double emin = TMath::Power(fEnergyBins[0],1.0-alpha);
144  double emax = TMath::Power(fEnergyBins[fNumEnergyBins],1.0-alpha);
145  Ev = TMath::Power(emin+(emax-emin)*rnd->RndFlux().Rndm(),1.0/(1.0-alpha));
146  costheta = -1+2*rnd->RndFlux().Rndm();
147  phi = 2.*kPi* rnd->RndFlux().Rndm();
148 
149  unsigned int nnu = fPdgCList->size();
150  unsigned int inu = rnd->RndFlux().Integer(nnu);
151  nu_pdg = (*fPdgCList)[inu];
152 
153  if(Ev < fEnergyBins[0]) {
154  LOG("Flux", pINFO) << "E < Emin";
155  return false;
156  }
157  double flux = this->GetFlux(nu_pdg, Ev, costheta, phi);
158  if(flux<=0) {
159  LOG("Flux", pINFO) << "Flux <= 0";
160  return false;
161  }
162  weight = flux*TMath::Power(Ev,alpha);
163  }
164  else {
165 
166  //
167  // generate nominal flux
168  //
169 
170  Axis_t ax = 0, ay = 0, az = 0;
171  fTotalFluxHisto->GetRandom3(ax, ay, az);
172  Ev = (double)ax;
173  costheta = (double)ay;
174  phi = (double)az;
175  nu_pdg = this->SelectNeutrino(Ev, costheta, phi);
176  weight = 1.0;
177  }
178 
179  // Compute etc trigonometric numbers
180  double sintheta = TMath::Sqrt(1-costheta*costheta);
181  double cosphi = TMath::Cos(phi);
182  double sinphi = TMath::Sin(phi);
183 
184  // Set the neutrino pdg code
185  fgPdgC = nu_pdg;
186 
187  // Set the neutrino weight
188  fWeight = weight;
189 
190  // Compute the neutrino momentum
191  // The `-1' means it is directed towards the detector.
192  double pz = -1.* Ev * costheta;
193  double py = -1.* Ev * sintheta * sinphi;
194  double px = -1.* Ev * sintheta * cosphi;
195 
196  // Default vertex is at the origin
197  double z = 0.0;
198  double y = 0.0;
199  double x = 0.0;
200 
201  // Shift the neutrino position onto the flux generation surface.
202  // The position is computed at the surface of a sphere with R=fRl
203  // at the topocentric horizontal (THZ) coordinate system.
204  if( fRl>0.0 ){
205  z += fRl * costheta;
206  y += fRl * sintheta * sinphi;
207  x += fRl * sintheta * cosphi;
208  }
209 
210  // Apply user-defined rotation from THZ -> user-defined topocentric
211  // coordinate system.
212  if( !fRotTHz2User.IsIdentity() )
213  {
214  TVector3 tx3(x, y, z );
215  TVector3 tp3(px,py,pz);
216 
217  tx3 = fRotTHz2User * tx3;
218  tp3 = fRotTHz2User * tp3;
219 
220  x = tx3.X();
221  y = tx3.Y();
222  z = tx3.Z();
223  px = tp3.X();
224  py = tp3.Y();
225  pz = tp3.Z();
226  }
227 
228  // If the position is left as is, then all generated neutrinos
229  // would point towards the origin.
230  // Displace the position randomly on the surface that is
231  // perpendicular to the selected point P(xo,yo,zo) on the sphere
232  if( fRt>0.0 ){
233  TVector3 vec(x,y,z); // vector towards selected point
234  TVector3 dvec1 = vec.Orthogonal(); // orthogonal vector
235  TVector3 dvec2 = dvec1; // second orthogonal vector
236  dvec2.Rotate(-kPi/2.0,vec); // rotate second vector by 90deg,
237  // now forming a new orthogonal cartesian coordinate system
238  double psi = 2.*kPi* rnd->RndFlux().Rndm(); // rndm angle [0,2pi]
239  double random = rnd->RndFlux().Rndm(); // rndm number [0,1]
240  dvec1.SetMag(TMath::Sqrt(random)*fRt*TMath::Cos(psi));
241  dvec2.SetMag(TMath::Sqrt(random)*fRt*TMath::Sin(psi));
242  x += dvec1.X() + dvec2.X();
243  y += dvec1.Y() + dvec2.Y();
244  z += dvec1.Z() + dvec2.Z();
245  }
246 
247  // Set the neutrino momentum and position 4-vectors with values
248  // calculated at previous steps.
249  fgP4.SetPxPyPzE(px, py, pz, Ev);
250  fgX4.SetXYZT (x, y, z, 0.);
251 
252  // Increment flux neutrino counter used for sample normalization purposes.
253  fNNeutrinos++;
254 
255  // Report and exit
256  LOG("Flux", pINFO)
257  << "Generated neutrino: "
258  << "\n pdg-code: " << fgPdgC
259  << "\n p4: " << utils::print::P4AsShortString(&fgP4)
260  << "\n x4: " << utils::print::X4AsString(&fgX4);
261 
262  return true;
263 }
264 //___________________________________________________________________________
265 long int GAtmoFlux::NFluxNeutrinos(void) const
266 {
267  return fNNeutrinos;
268 }
269 //___________________________________________________________________________
271 {
272  emin = TMath::Max(0., emin);
273  fMinEvCut = emin;
274 }
275 //___________________________________________________________________________
277 {
278  emax = TMath::Max(0., emax);
279  fMaxEvCut = emax;
280 }
281 //___________________________________________________________________________
282 void GAtmoFlux::Clear(Option_t * opt)
283 {
284 // Dummy clear method needed to conform to GFluxI interface
285 //
286  LOG("Flux", pERROR) << "No clear method implemented for option:"<< opt;
287 }
288 //___________________________________________________________________________
289 void GAtmoFlux::GenerateWeighted(bool gen_weighted)
290 {
291  fGenWeighted = gen_weighted;
292 }
293 //___________________________________________________________________________
295 {
296  if( index != 1.0 ){
297  fSpectralIndex = index;
298  }
299  else {
300  LOG("Flux", pWARN) << "Warning: cannot use a spectral index of unity";
301  }
302 
303  LOG("Flux", pNOTICE) << "Using Spectral Index = " << index;
304 }
305 //___________________________________________________________________________
306 void GAtmoFlux::SetUserCoordSystem(TRotation & rotation)
307 {
308  fRotTHz2User = rotation;
309 }
310 //___________________________________________________________________________
312 {
313  LOG("Flux", pNOTICE) << "Initializing atmospheric flux driver";
314 
315  fMaxEv = -1;
316 
317  fNumPhiBins = 0;
318  fNumCosThetaBins = 0;
319  fNumEnergyBins = 0;
320  fPhiBins = 0;
321  fCosThetaBins = 0;
322  fEnergyBins = 0;
323 
324  fTotalFluxHisto = 0;
325  fTotalFluxHistoIntg = 0;
326 
327  bool allow_dup = false;
328  fPdgCList = new PDGCodeList(allow_dup);
329 
330  // initializing flux TH3D histos [ flux = f(Ev,costheta,phi) ] & files
331  fFluxFile.clear();
332  fFluxHistoMap.clear();
333  fTotalFluxHisto = 0;
334  fTotalFluxHistoIntg = 0;
335 
336  // Default option is to generate unweighted flux neutrinos
337  // (flux = f(E,costheta) will be used as PDFs)
338  // User can enable option to generate weighted neutrinos
339  // (neutrinos will be generated uniformly over costheta,
340  // and using a power law function in neutrino energy.
341  // The input flux = f(E,costheta) will be used for calculating a weight).
342  // Using a weighted flux avoids statistical fluctuations at high energies.
343  fSpectralIndex = 2.0;
344 
345  // weighting switched off by default
346  this->GenerateWeighted(false);
347 
348  // Default: No min/max energy cut
349  this->ForceMinEnergy(0.);
350  this->ForceMaxEnergy(9999999999.);
351 
352  // Default radii
353  fRl = 0.0;
354  fRt = 0.0;
355 
356  // Default detector coord system: Topocentric Horizontal Coordinate system
357  fRotTHz2User.SetToIdentity();
358 
359  // Reset `current' selected flux neutrino
360  this->ResetSelection();
361 
362  // Reset number of neutrinos thrown so far
363  fNNeutrinos = 0;
364 
365  // Done!
366  fInitialized = 1;
367 }
368 //___________________________________________________________________________
370 {
371 // initializing running neutrino pdg-code, 4-position, 4-momentum
372 
373  fgPdgC = 0;
374  fgP4.SetPxPyPzE (0.,0.,0.,0.);
375  fgX4.SetXYZT (0.,0.,0.,0.);
376  fWeight = 0;
377 }
378 //___________________________________________________________________________
380 {
381  LOG("Flux", pNOTICE) << "Cleaning up...";
382 
383  map<int,TH3D*>::iterator rawiter = fRawFluxHistoMap.begin();
384  for( ; rawiter != fRawFluxHistoMap.end(); ++rawiter) {
385  TH3D * flux_histogram = rawiter->second;
386  if(flux_histogram) {
387  delete flux_histogram;
388  flux_histogram = 0;
389  }
390  }
391  fRawFluxHistoMap.clear();
392 
393  map<int,TH3D*>::iterator iter = fFluxHistoMap.begin();
394  for( ; iter != fFluxHistoMap.end(); ++iter) {
395  TH3D * flux_histogram = iter->second;
396  if(flux_histogram) {
397  delete flux_histogram;
398  flux_histogram = 0;
399  }
400  }
401  fFluxHistoMap.clear();
402 
403  if (fTotalFluxHisto) delete fTotalFluxHisto;
404  if (fPdgCList) delete fPdgCList;
405 
406  if (fPhiBins ) { delete[] fPhiBins ; fPhiBins =NULL; }
407  if (fCosThetaBins) { delete[] fCosThetaBins; fCosThetaBins=NULL; }
408  if (fEnergyBins ) { delete[] fEnergyBins ; fEnergyBins =NULL; }
409 
410 }
411 //___________________________________________________________________________
412 void GAtmoFlux::SetRadii(double Rlongitudinal, double Rtransverse)
413 {
414  LOG ("Flux", pNOTICE) << "Setting R[longitudinal] = " << Rlongitudinal;
415  LOG ("Flux", pNOTICE) << "Setting R[transverse] = " << Rtransverse;
416 
417  fRl = Rlongitudinal;
418  fRt = Rtransverse;
419 }
420 //___________________________________________________________________________
421 void GAtmoFlux::AddFluxFile(int nu_pdg, string filename)
422 {
423  // Check file exists
424  std::ifstream f(filename.c_str());
425  if (!f.good()) {
426  LOG("Flux", pFATAL) << "Flux file does not exist "<<filename;
427  exit(-1);
428  }
429  if ( pdg::IsNeutrino(nu_pdg) || pdg::IsAntiNeutrino(nu_pdg) ) {
430  fFluxFlavour.push_back(nu_pdg);
431  fFluxFile.push_back(filename);
432  } else {
433  LOG ("Flux", pWARN)
434  << "Input particle code: " << nu_pdg << " not a neutrino!";
435  }
436 }
437 //___________________________________________________________________________
439 {
440 // FLUKA and BGLRS provide one file per neutrino species.
441 // HAKKKM provides a single file for all nue,nuebar,numu,numubar.
442 // If no neutrino species is provided, assume that the file contains all 4
443 // but fit it into the franework developed for FLUKA and BGLRS,
444 // i.e. add the file 4 times
445 
446 // Check file exists
447  std::ifstream f(filename.c_str());
448  if (!f.good()) {
449  LOG("Flux", pFATAL) << "Flux file does not exist "<<filename;
450  exit(-1);
451  }
452 
453  fFluxFlavour.push_back(kPdgNuE); fFluxFile.push_back(filename);
454  fFluxFlavour.push_back(kPdgAntiNuE); fFluxFile.push_back(filename);
455  fFluxFlavour.push_back(kPdgNuMu); fFluxFile.push_back(filename);
456  fFluxFlavour.push_back(kPdgAntiNuMu); fFluxFile.push_back(filename);
457 
458 }
459 //___________________________________________________________________________
460 //void GAtmoFlux::SetFluxFile(int nu_pdg, string filename)
461 //{
462 // return AddFluxFile( nu_pdg, filename );
463 //}
464 //___________________________________________________________________________
466 {
467  LOG("Flux", pNOTICE)
468  << "Loading atmospheric neutrino flux simulation data";
469 
470  fFluxHistoMap.clear();
471  fPdgCList->clear();
472 
473  bool loading_status = true;
474 
475  for( unsigned int n=0; n<fFluxFlavour.size(); n++ ){
476  int nu_pdg = fFluxFlavour.at(n);
477  string filename = fFluxFile.at(n);
478  string pname = PDGLibrary::Instance()->Find(nu_pdg)->GetName();
479 
480  LOG("Flux", pNOTICE) << "Loading data for: " << pname;
481 
482  // create histogram
483  TH3D* hist = 0;
484  std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
485  if( myMapEntry == fRawFluxHistoMap.end() ){
486 // hist = myMapEntry->second;
487 // if(hist==0) {
488  hist = this->CreateFluxHisto(pname.c_str(), pname.c_str());
489  fRawFluxHistoMap.insert( map<int,TH3D*>::value_type(nu_pdg,hist) );
490 // }
491  }
492  // now let concrete instances to read the flux-specific data files
493  // and fill the histogram
494  bool loaded = this->FillFluxHisto(nu_pdg, filename);
495  loading_status = loading_status && loaded;
496  }
497 
498  if(loading_status) {
499  map<int,TH3D*>::iterator hist_iter = fRawFluxHistoMap.begin();
500  for ( ; hist_iter != fRawFluxHistoMap.end(); ++hist_iter) {
501  int nu_pdg = hist_iter->first;
502  TH3D* hist = hist_iter->second;
503 
504  TH3D* hnorm = this->CreateNormalisedFluxHisto( hist );
505  fFluxHistoMap.insert( map<int,TH3D*>::value_type(nu_pdg,hnorm) );
506  fPdgCList->push_back(nu_pdg);
507  }
508 
509  LOG("Flux", pNOTICE)
510  << "Atmospheric neutrino flux simulation data loaded!";
511  this->AddAllFluxes();
512  return true;
513  }
514 
515  LOG("Flux", pERROR)
516  << "Error loading atmospheric neutrino flux simulation data";
517  return false;
518 }
519 //___________________________________________________________________________
521 {
522 // return integrated flux
523 
524  // sanity check
525  if(!hist) return 0;
526 
527  // make new histogram name
528  TString histname = hist->GetName();
529  histname.Append("_IntegratedFlux");
530 
531  // make new histogram
532  TH3D* hist_intg = (TH3D*)(hist->Clone(histname.Data()));
533  hist_intg->Reset();
534 
535  // integrate flux in each bin
536  Double_t dN_dEdS = 0.0;
537  Double_t dS = 0.0;
538  Double_t dE = 0.0;
539  Double_t dN = 0.0;
540 
541  for(Int_t e_bin = 1; e_bin <= hist->GetXaxis()->GetNbins(); e_bin++)
542  {
543  for(Int_t c_bin = 1; c_bin <= hist->GetYaxis()->GetNbins(); c_bin++)
544  {
545  for(Int_t p_bin = 1; p_bin <= hist->GetZaxis()->GetNbins(); p_bin++)
546  {
547  dN_dEdS = hist->GetBinContent(e_bin,c_bin,p_bin);
548 
549  dE = hist->GetXaxis()->GetBinUpEdge(e_bin)
550  - hist->GetXaxis()->GetBinLowEdge(e_bin);
551 
552  dS = ( hist->GetZaxis()->GetBinUpEdge(p_bin)
553  - hist->GetZaxis()->GetBinLowEdge(p_bin) )
554  * ( hist->GetYaxis()->GetBinUpEdge(c_bin)
555  - hist->GetYaxis()->GetBinLowEdge(c_bin) );
556 
557  dN = dN_dEdS*dE*dS;
558 
559  hist_intg->SetBinContent(e_bin,c_bin,p_bin,dN);
560  }
561  }
562  }
563 
564  return hist_intg;
565 }
566 //___________________________________________________________________________
568 {
569  LOG("Flux", pDEBUG) << "Forcing flux histogram contents to 0";
570  if(!histo) return;
571  histo->Reset();
572 }
573 //___________________________________________________________________________
575 {
576  LOG("Flux", pNOTICE)
577  << "Computing combined flux & flux normalization factor";
578 
579  if(fTotalFluxHisto) delete fTotalFluxHisto;
580 
581  fTotalFluxHisto = this->CreateFluxHisto("sum", "combined flux" );
582 
583  map<int,TH3D*>::iterator it = fFluxHistoMap.begin();
584  for( ; it != fFluxHistoMap.end(); ++it) {
585  TH3D * flux_histogram = it->second;
586  fTotalFluxHisto->Add(flux_histogram);
587  }
588 
589  fTotalFluxHistoIntg = fTotalFluxHisto->Integral();
590 }
591 //___________________________________________________________________________
592 TH3D * GAtmoFlux::CreateFluxHisto(string name, string title)
593 {
594  LOG("Flux", pNOTICE) << "Instantiating histogram: [" << name << "]";
595  TH3D * hist = new TH3D(
596  name.c_str(), title.c_str(),
597  fNumEnergyBins, fEnergyBins,
598  fNumCosThetaBins, fCosThetaBins,
599  fNumPhiBins, fPhiBins);
600  return hist;
601 }
602 //___________________________________________________________________________
603 int GAtmoFlux::SelectNeutrino(double Ev, double costheta, double phi)
604 {
605 // Select a neutrino species at the input Ev and costheta given their
606 // relatve flux at this bin.
607 // Returns a neutrino PDG code
608 
609  unsigned int n = fPdgCList->size();
610  double * flux = new double[n];
611 
612  unsigned int i=0;
613  map<int,TH3D*>::iterator it = fFluxHistoMap.begin();
614  for( ; it != fFluxHistoMap.end(); ++it) {
615  TH3D * flux_histogram = it->second;
616  int ibin = flux_histogram->FindBin(Ev,costheta,phi);
617  flux[i] = flux_histogram->GetBinContent(ibin);
618  i++;
619  }
620  double flux_sum = 0;
621  for(i=0; i<n; i++) {
622  flux_sum += flux[i];
623  flux[i] = flux_sum;
624  LOG("Flux", pDEBUG)
625  << "Sum{Flux(0->" << i <<")} = " << flux[i];
626  }
627 
629  double R = flux_sum * rnd->RndFlux().Rndm();
630 
631  LOG("Flux", pDEBUG) << "R = " << R;
632 
633  for(i=0; i<n; i++) {
634  if( R < flux[i] ) {
635  delete [] flux;
636  int selected_pdg = (*fPdgCList)[i];
637  LOG("Flux", pINFO)
638  << "Selected neutrino PDG code = " << selected_pdg;
639  return selected_pdg;
640  }
641  }
642 
643  // shouldn't be here
644  LOG("Flux", pERROR) << "Could not select a neutrino species!";
645  assert(false);
646 
647  return -1;
648 }
649 
650 //___________________________________________________________________________
651 TH3D* GAtmoFlux::GetFluxHistogram(int flavour)
652 {
653  TH3D* histogram = 0;
654  std::map<int,TH3D*>::iterator it = fRawFluxHistoMap.find(flavour);
655  if(it != fRawFluxHistoMap.end())
656  {
657  histogram = it->second;
658  }
659  return histogram;
660 }
661 //___________________________________________________________________________
662 double GAtmoFlux::GetFlux(int flavour)
663 {
664  TH3D* flux_hist = this->GetFluxHistogram(flavour);
665  if(!flux_hist) return 0.0;
666 
667  Double_t flux = 0.0;
668  Double_t dN_dEdS = 0.0;
669  Double_t dS = 0.0;
670  Double_t dE = 0.0;
671 
672  for(Int_t e_bin = 1; e_bin <= flux_hist->GetXaxis()->GetNbins(); e_bin++)
673  {
674  for(Int_t c_bin = 1; c_bin <= flux_hist->GetYaxis()->GetNbins(); c_bin++)
675  {
676  for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
677  {
678  dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
679 
680  dE = flux_hist->GetXaxis()->GetBinUpEdge(e_bin)
681  - flux_hist->GetXaxis()->GetBinLowEdge(e_bin);
682 
683  dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
684  - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) )
685  * ( flux_hist->GetYaxis()->GetBinUpEdge(c_bin)
686  - flux_hist->GetYaxis()->GetBinLowEdge(c_bin) );
687 
688  flux += dN_dEdS*dE*dS;
689  }
690  }
691  }
692 
693  return flux;
694 }
695 //___________________________________________________________________________
696 double GAtmoFlux::GetFlux(int flavour, double energy)
697 {
698  TH3D* flux_hist = this->GetFluxHistogram(flavour);
699  if(!flux_hist) return 0.0;
700 
701  Double_t flux = 0.0;
702  Double_t dN_dEdS = 0.0;
703  Double_t dS = 0.0;
704 
705  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
706 
707  for(Int_t c_bin = 1; c_bin <= flux_hist->GetYaxis()->GetNbins(); c_bin++)
708  {
709  for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
710  {
711  dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
712 
713  dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
714  - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) )
715  * ( flux_hist->GetYaxis()->GetBinUpEdge(c_bin)
716  - flux_hist->GetYaxis()->GetBinLowEdge(c_bin) );
717 
718  flux += dN_dEdS*dS;
719  }
720  }
721 
722  return flux;
723 }
724 //___________________________________________________________________________
725 double GAtmoFlux::GetFlux(int flavour, double energy, double costh)
726 {
727  TH3D* flux_hist = this->GetFluxHistogram(flavour);
728  if(!flux_hist) return 0.0;
729 
730  Double_t flux = 0.0;
731  Double_t dN_dEdS = 0.0;
732  Double_t dS = 0.0;
733 
734  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
735  Int_t c_bin = flux_hist->GetYaxis()->FindBin(costh);
736 
737  for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
738  {
739  dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
740 
741  dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
742  - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) );
743 
744  flux += dN_dEdS*dS;
745  }
746 
747  return flux;
748 }
749 //___________________________________________________________________________
750 double GAtmoFlux::GetFlux(int flavour, double energy, double costh, double phi)
751 {
752  TH3D* flux_hist = this->GetFluxHistogram(flavour);
753  if(!flux_hist) return 0.0;
754 
755  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
756  Int_t c_bin = flux_hist->GetYaxis()->FindBin(costh);
757  Int_t p_bin = flux_hist->GetZaxis()->FindBin(phi);
758 
759  Double_t flux = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
760  return flux;
761 }
762 //___________________________________________________________________________
int SelectNeutrino(double Ev, double costheta, double phi)
Definition: GAtmoFlux.cxx:603
const double kPi
const XML_Char * name
Definition: expat.h:151
Basic constants.
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:52
const int kPdgNuE
Definition: PDGCodes.h:28
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
set< int >::iterator it
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const Var weight
TH3D * CreateNormalisedFluxHisto(TH3D *hist)
Definition: GAtmoFlux.cxx:520
const int kPdgAntiNuE
Definition: PDGCodes.h:29
#define pFATAL
Definition: Messenger.h:57
TH3D * GetFluxHistogram(int flavour)
Definition: GAtmoFlux.cxx:651
const int kPdgNuMu
Definition: PDGCodes.h:30
long int NFluxNeutrinos(void) const
Number of flux nu&#39;s generated. Not the same as the number of nu&#39;s thrown towards the geometry (if the...
Definition: GAtmoFlux.cxx:265
string filename
Definition: shutoffs.py:106
void SetRadii(double Rlongitudinal, double Rtransverse)
Definition: GAtmoFlux.cxx:412
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
double GetFlux(int flavour)
Definition: GAtmoFlux.cxx:662
virtual void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
Definition: GAtmoFlux.cxx:289
double dE
Loaders::FluxType flux
A list of PDG codes.
Definition: PDGCodeList.h:33
GFluxI * GetFlux(void)
Definition: gAtmoEvGen.cxx:441
const double emin
GENIE flux drivers.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void ForceMaxEnergy(double emax)
Definition: GAtmoFlux.cxx:276
bool fInitialized
Float_t E
Definition: plot.C:20
const double emax
TH1F * hnorm
Definition: plots_total.C:4
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
double energy
Definition: plottest35.C:25
void SetSpectralIndex(double index)
Definition: GAtmoFlux.cxx:294
#define R(x)
TH2D * histo
#define pINFO
Definition: Messenger.h:63
Eigen::VectorXd vec
string pname
Definition: eplot.py:33
virtual double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
Definition: GAtmoFlux.cxx:87
z
Definition: test.py:28
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
#define pWARN
Definition: Messenger.h:61
void SetUserCoordSystem(TRotation &rotation)
Rotation: Topocentric Horizontal -> User-defined Topocentric Coord System.
Definition: GAtmoFlux.cxx:306
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
virtual bool GenerateNext(void)
generate the next flux neutrino (return false in err)
Definition: GAtmoFlux.cxx:92
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
void ZeroFluxHisto(TH3D *hist)
Definition: GAtmoFlux.cxx:567
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
void AddFluxFile(int neutrino_pdg, string filename)
Definition: GAtmoFlux.cxx:421
exit(0)
assert(nhit_max >=nhit_nbins)
double Emax
TH3D * CreateFluxHisto(string name, string title)
Definition: GAtmoFlux.cxx:592
void ResetSelection(void)
Definition: GAtmoFlux.cxx:369
virtual void Clear(Option_t *opt)
reset state variables based on opt
Definition: GAtmoFlux.cxx:282
#define pNOTICE
Definition: Messenger.h:62
bool GenerateNext_1try(void)
Definition: GAtmoFlux.cxx:114
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:72
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:64
void ForceMinEnergy(double emin)
Definition: GAtmoFlux.cxx:270
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define pDEBUG
Definition: Messenger.h:64