72 using namespace genie;
89 return TMath::Min(fMaxEv, fMaxEvCut);
96 bool nextok = this->GenerateNext_1try();
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();
108 bool accept = (E<=Emax && E>=Emin && wght>0);
109 if(accept)
return true;
120 this->ResetSelection();
127 double costheta = 0.;
141 double alpha = fSpectralIndex;
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();
149 unsigned int nnu = fPdgCList->size();
150 unsigned int inu = rnd->
RndFlux().Integer(nnu);
151 nu_pdg = (*fPdgCList)[inu];
153 if(Ev < fEnergyBins[0]) {
157 double flux = this->
GetFlux(nu_pdg, Ev, costheta, phi);
162 weight = flux*TMath::Power(Ev,alpha);
170 Axis_t
ax = 0, ay = 0, az = 0;
171 fTotalFluxHisto->GetRandom3(ax, ay, az);
173 costheta = (double)ay;
175 nu_pdg = this->SelectNeutrino(Ev, costheta, phi);
180 double sintheta =
TMath::Sqrt(1-costheta*costheta);
181 double cosphi = TMath::Cos(phi);
182 double sinphi = TMath::Sin(phi);
192 double pz = -1.* Ev * costheta;
193 double py = -1.* Ev * sintheta * sinphi;
194 double px = -1.* Ev * sintheta * cosphi;
206 y += fRl * sintheta * sinphi;
207 x += fRl * sintheta * cosphi;
212 if( !fRotTHz2User.IsIdentity() )
214 TVector3 tx3(x, y, z );
215 TVector3 tp3(px,py,pz);
217 tx3 = fRotTHz2User * tx3;
218 tp3 = fRotTHz2User * tp3;
234 TVector3 dvec1 = vec.Orthogonal();
235 TVector3 dvec2 = dvec1;
236 dvec2.Rotate(-
kPi/2.0,vec);
239 double random = rnd->
RndFlux().Rndm();
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();
249 fgP4.SetPxPyPzE(px, py, pz, Ev);
250 fgX4.SetXYZT (x, y, z, 0.);
257 <<
"Generated neutrino: " 258 <<
"\n pdg-code: " << fgPdgC
272 emin = TMath::Max(0., emin);
278 emax = TMath::Max(0., emax);
286 LOG(
"Flux",
pERROR) <<
"No clear method implemented for option:"<<
opt;
291 fGenWeighted = gen_weighted;
297 fSpectralIndex =
index;
300 LOG(
"Flux",
pWARN) <<
"Warning: cannot use a spectral index of unity";
308 fRotTHz2User = rotation;
313 LOG(
"Flux",
pNOTICE) <<
"Initializing atmospheric flux driver";
318 fNumCosThetaBins = 0;
325 fTotalFluxHistoIntg = 0;
327 bool allow_dup =
false;
332 fFluxHistoMap.clear();
334 fTotalFluxHistoIntg = 0;
343 fSpectralIndex = 2.0;
346 this->GenerateWeighted(
false);
349 this->ForceMinEnergy(0.);
350 this->ForceMaxEnergy(9999999999.);
357 fRotTHz2User.SetToIdentity();
360 this->ResetSelection();
374 fgP4.SetPxPyPzE (0.,0.,0.,0.);
375 fgX4.SetXYZT (0.,0.,0.,0.);
383 map<int,TH3D*>::iterator rawiter = fRawFluxHistoMap.begin();
384 for( ; rawiter != fRawFluxHistoMap.end(); ++rawiter) {
385 TH3D * flux_histogram = rawiter->second;
387 delete flux_histogram;
391 fRawFluxHistoMap.clear();
393 map<int,TH3D*>::iterator iter = fFluxHistoMap.begin();
394 for( ; iter != fFluxHistoMap.end(); ++iter) {
395 TH3D * flux_histogram = iter->second;
397 delete flux_histogram;
401 fFluxHistoMap.clear();
403 if (fTotalFluxHisto)
delete fTotalFluxHisto;
404 if (fPdgCList)
delete fPdgCList;
406 if (fPhiBins ) {
delete[] fPhiBins ; fPhiBins =NULL; }
407 if (fCosThetaBins) {
delete[] fCosThetaBins; fCosThetaBins=NULL; }
408 if (fEnergyBins ) {
delete[] fEnergyBins ; fEnergyBins =NULL; }
414 LOG (
"Flux",
pNOTICE) <<
"Setting R[longitudinal] = " << Rlongitudinal;
415 LOG (
"Flux",
pNOTICE) <<
"Setting R[transverse] = " << Rtransverse;
424 std::ifstream
f(filename.c_str());
430 fFluxFlavour.push_back(nu_pdg);
431 fFluxFile.push_back(filename);
434 <<
"Input particle code: " << nu_pdg <<
" not a neutrino!";
447 std::ifstream
f(filename.c_str());
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);
468 <<
"Loading atmospheric neutrino flux simulation data";
470 fFluxHistoMap.clear();
473 bool loading_status =
true;
475 for(
unsigned int n=0;
n<fFluxFlavour.size();
n++ ){
476 int nu_pdg = fFluxFlavour.at(
n);
484 std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
485 if( myMapEntry == fRawFluxHistoMap.end() ){
488 hist = this->CreateFluxHisto(pname.c_str(), pname.c_str());
489 fRawFluxHistoMap.insert( map<int,TH3D*>::value_type(nu_pdg,hist) );
494 bool loaded = this->FillFluxHisto(nu_pdg, filename);
495 loading_status = loading_status && loaded;
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;
504 TH3D*
hnorm = this->CreateNormalisedFluxHisto( hist );
505 fFluxHistoMap.insert( map<int,TH3D*>::value_type(nu_pdg,hnorm) );
506 fPdgCList->push_back(nu_pdg);
510 <<
"Atmospheric neutrino flux simulation data loaded!";
511 this->AddAllFluxes();
516 <<
"Error loading atmospheric neutrino flux simulation data";
528 TString histname = hist->GetName();
529 histname.Append(
"_IntegratedFlux");
532 TH3D* hist_intg = (TH3D*)(hist->Clone(histname.Data()));
536 Double_t dN_dEdS = 0.0;
541 for(Int_t e_bin = 1; e_bin <= hist->GetXaxis()->GetNbins(); e_bin++)
543 for(Int_t c_bin = 1; c_bin <= hist->GetYaxis()->GetNbins(); c_bin++)
545 for(Int_t p_bin = 1; p_bin <= hist->GetZaxis()->GetNbins(); p_bin++)
547 dN_dEdS = hist->GetBinContent(e_bin,c_bin,p_bin);
549 dE = hist->GetXaxis()->GetBinUpEdge(e_bin)
550 - hist->GetXaxis()->GetBinLowEdge(e_bin);
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) );
559 hist_intg->SetBinContent(e_bin,c_bin,p_bin,dN);
569 LOG(
"Flux",
pDEBUG) <<
"Forcing flux histogram contents to 0";
577 <<
"Computing combined flux & flux normalization factor";
579 if(fTotalFluxHisto)
delete fTotalFluxHisto;
581 fTotalFluxHisto = this->CreateFluxHisto(
"sum",
"combined flux" );
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);
589 fTotalFluxHistoIntg = fTotalFluxHisto->Integral();
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);
609 unsigned int n = fPdgCList->size();
610 double *
flux =
new double[
n];
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);
625 <<
"Sum{Flux(0->" << i <<
")} = " << flux[
i];
629 double R = flux_sum * rnd->
RndFlux().Rndm();
636 int selected_pdg = (*fPdgCList)[
i];
638 <<
"Selected neutrino PDG code = " << selected_pdg;
644 LOG(
"Flux",
pERROR) <<
"Could not select a neutrino species!";
654 std::map<int,TH3D*>::iterator
it = fRawFluxHistoMap.find(flavour);
655 if(it != fRawFluxHistoMap.end())
657 histogram = it->second;
664 TH3D* flux_hist = this->GetFluxHistogram(flavour);
665 if(!flux_hist)
return 0.0;
668 Double_t dN_dEdS = 0.0;
672 for(Int_t e_bin = 1; e_bin <= flux_hist->GetXaxis()->GetNbins(); e_bin++)
674 for(Int_t c_bin = 1; c_bin <= flux_hist->GetYaxis()->GetNbins(); c_bin++)
676 for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
678 dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
680 dE = flux_hist->GetXaxis()->GetBinUpEdge(e_bin)
681 - flux_hist->GetXaxis()->GetBinLowEdge(e_bin);
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) );
688 flux += dN_dEdS*dE*dS;
698 TH3D* flux_hist = this->GetFluxHistogram(flavour);
699 if(!flux_hist)
return 0.0;
702 Double_t dN_dEdS = 0.0;
705 Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
707 for(Int_t c_bin = 1; c_bin <= flux_hist->GetYaxis()->GetNbins(); c_bin++)
709 for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
711 dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
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) );
727 TH3D* flux_hist = this->GetFluxHistogram(flavour);
728 if(!flux_hist)
return 0.0;
731 Double_t dN_dEdS = 0.0;
734 Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
735 Int_t c_bin = flux_hist->GetYaxis()->FindBin(costh);
737 for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
739 dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
741 dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
742 - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) );
752 TH3D* flux_hist = this->GetFluxHistogram(flavour);
753 if(!flux_hist)
return 0.0;
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);
759 Double_t
flux = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
int SelectNeutrino(double Ev, double costheta, double phi)
string P4AsShortString(const TLorentzVector *p)
bool IsNeutrino(int pdgc)
THE MAIN GENIE PROJECT NAMESPACE
static RandomGen * Instance()
Access instance.
TH3D * CreateNormalisedFluxHisto(TH3D *hist)
TH3D * GetFluxHistogram(int flavour)
long int NFluxNeutrinos(void) const
Number of flux nu's generated. Not the same as the number of nu's thrown towards the geometry (if the...
void SetRadii(double Rlongitudinal, double Rtransverse)
A singleton holding random number generator classes. All random number generation in GENIE should tak...
double GetFlux(int flavour)
virtual void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
void ForceMaxEnergy(double emax)
bool IsAntiNeutrino(int pdgc)
void SetSpectralIndex(double index)
virtual double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
void SetUserCoordSystem(TRotation &rotation)
Rotation: Topocentric Horizontal -> User-defined Topocentric Coord System.
static PDGLibrary * Instance(void)
virtual bool GenerateNext(void)
generate the next flux neutrino (return false in err)
Var Sqrt(const Var &v)
Use to take sqrt of a var.
void ZeroFluxHisto(TH3D *hist)
TParticlePDG * Find(int pdgc)
void AddFluxFile(int neutrino_pdg, string filename)
assert(nhit_max >=nhit_nbins)
TH3D * CreateFluxHisto(string name, string title)
void ResetSelection(void)
virtual void Clear(Option_t *opt)
reset state variables based on opt
bool GenerateNext_1try(void)
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
string X4AsString(const TLorentzVector *x)
void ForceMinEnergy(double emin)
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...