GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GFluxBlender.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Robert Hatcher <rhatcher@fnal.gov>
7  Fermi National Accelerator Laboratory
8 */
9 //____________________________________________________________________________
10 
11 #include <math.h>
12 #include <iostream>
13 #include <iomanip>
14 
15 //GENIE includes
17 #include "Tools/Flux/GNuMIFlux.h"
20 
24 #define LOG_BEGIN(a,b) LOG(a,b)
25 #define LOG_END ""
26 
27 namespace genie {
28 namespace flux {
29 
30 //____________________________________________________________________________
31 
33  GFluxI(),
34  fRealGFluxI(0),
35  fGNuMIFlux(0),
36  fGSimpleFlux(0),
37  fFlavorMixer(0),
38  fPDGListGenerator(),
39  fPDGListMixed(),
40  fNPDGOut(0),
41  fBaselineDist(0),
42  fEnergy(0),
43  fDistance(0),
44  fPdgCGenerated(0),
45  fPdgCMixed(0)
46 { ; }
47 
49 {
50  if ( fRealGFluxI ) { delete fRealGFluxI; fRealGFluxI = 0; }
51  if ( fFlavorMixer ) { delete fFlavorMixer; fFlavorMixer = 0; }
52 }
53 
54 //____________________________________________________________________________
56 {
57  /// Return (reference to) list of possible neutrinos *after* mixing.
58  /// These are PDG code that might be interacted.
59 
60  // this is only ever called by GMCJDriver::GetParticleLists()
61  // during GMCJDriver::Configure() which should happen just once
62  // so don't try to be too clever.
64 
65  // okay, really stupid at this time ...
66  fPDGListMixed.clear();
73 
74  // size the probability arrays to the same number of entries
75  fNPDGOut = fPDGListMixed.size();
76  fProb.resize(fNPDGOut);
77  fSumProb.resize(fNPDGOut);
78 
79  if ( ! fFlavorMixer ) return fRealGFluxI->FluxParticles();
80  else return fPDGListMixed;
81 }
82 
83 //____________________________________________________________________________
85 {
86 
87  bool gen1 = false;
88  while ( ! gen1 ) {
89  if ( ! fRealGFluxI->GenerateNext() ) return false;
90  // have a new entry
92  if ( ! fFlavorMixer ) {
93  // simple case when not configured with a mixing model
95  gen1 = true;
96  } else {
97  // now pick a new flavor
101  fEnergy = fRealGFluxI->Momentum().Energy();
103  // we may have to generate a new neutrino if it oscillates away
104  // don't pass non-particles to GENIE ... it won't like it
105  gen1 = ( fPdgCMixed != 0 );
106  }
107  }
108  return true;
109 }
110 //____________________________________________________________________________
111 void GFluxBlender::Clear(Option_t * opt)
112 {
113 // Clear method needed to conform to GFluxI interface
114 //
115  fRealGFluxI->Clear(opt);
116 }
117 //____________________________________________________________________________
118 long int GFluxBlender::Index(void)
119 {
120 // Index method needed to conform to GFluxI interface
121 //
122  return fRealGFluxI->Index();
123 }
124 //____________________________________________________________________________
125 void GFluxBlender::GenerateWeighted(bool gen_weighted)
126 {
127 // Dummy implementation needed to conform to GFluxI interface
128 //
129  LOG("FluxBlender", pERROR) <<
130  "No GenerateWeighted(bool gen_weighted) method implemented for " <<
131  "gen_weighted: " << gen_weighted;
132 }
133 //____________________________________________________________________________
135 {
136  GFluxI* oldgen = fRealGFluxI;
138  // avoid re-casting
139  fGNuMIFlux = dynamic_cast<GNuMIFlux*>(fRealGFluxI);
140  fGSimpleFlux = dynamic_cast<GSimpleNtpFlux*>(fRealGFluxI);
141  // force evaluation of particle lists
142  this->FluxParticles();
143 
144  return oldgen;
145 }
146 
147 //____________________________________________________________________________
149 {
150  GFlavorMixerI* oldmix = fFlavorMixer;
151  fFlavorMixer = mixer;
152  return oldmix;
153 }
154 
155 //____________________________________________________________________________
156 int GFluxBlender::ChooseFlavor(int pdg_init, double energy, double dist)
157 {
158  // choose a new flavor
159  bool isset = false;
160  int pdg_out = 0;
161  double sumprob = 0;
162 
163  fRndm = RandomGen::Instance()->RndFlux().Rndm();
164  for (size_t indx = 0; indx < fNPDGOut; ++indx ) {
165  int pdg_test = fPDGListMixed[indx];
166  double prob = fFlavorMixer->Probability(pdg_init,pdg_test,energy,dist);
167  fProb[indx] = prob;
168  sumprob += fProb[indx];
169  fSumProb[indx] = sumprob;
170  if ( ! isset && fRndm < sumprob ) {
171  isset = true;
172  pdg_out = pdg_test;
173  }
174  }
175 
176  return pdg_out;
177 }
178 
179 //____________________________________________________________________________
181 {
182  LOG_BEGIN("FluxBlender", pINFO) << "GFluxBlender::PrintConfig()" << LOG_END;
183  if ( fRealGFluxI ) {
184  LOG_BEGIN("FluxBlender", pINFO)
185  << " fRealGFluxI is a \""
186  << typeid(fRealGFluxI).name() << "\""
187  << LOG_END;
188  } else {
189  LOG_BEGIN("FluxBlender", pINFO)
190  << " fRealGFluxI is not initialized" << LOG_END;
191  }
192  if ( fFlavorMixer ) {
193  LOG_BEGIN("FluxBlender", pINFO)
194  << " fFlavorMixer is a \""
195  << typeid(fFlavorMixer).name() << "\""
196  << LOG_END;
197  } else {
198  LOG_BEGIN("FluxBlender", pINFO)
199  << " fFlavorMixer is not initialized" << LOG_END;
200  }
201  LOG_BEGIN("FluxBlender", pINFO)
202  << " BaselineDist " << fBaselineDist << LOG_END;
203  LOG_BEGIN("FluxBlender", pINFO)
204  << "PDG List from Generator" << fPDGListGenerator << LOG_END;
205  LOG_BEGIN("FluxBlender", pINFO)
206  << "PDG List after mixing (n=" << fNPDGOut << ")"
207  << fPDGListMixed << LOG_END;
208 
209 }
210 
211 //____________________________________________________________________________
212 void GFluxBlender::PrintState(bool verbose)
213 {
214  LOG_BEGIN("FluxBlender", pINFO) << "GFluxBlender::PrintState()" << LOG_END;
215  LOG_BEGIN("FluxBlender", pINFO)
216  << " Flavor " << fPdgCGenerated
217  << " ==> " << fPdgCMixed
218  << " (E=" << fEnergy << ", dist=" << fDistance << ")" << LOG_END;
219  if ( verbose ) {
220  LOG_BEGIN("FluxBlender", pINFO) << " Rndm = " << fRndm << LOG_END;
221  for (size_t indx = 0; indx < fNPDGOut; ++indx ) {
222  LOG_BEGIN("FluxBlender", pINFO)
223  << " [" << indx << "] "
224  << std::setw(3) << fPdgCGenerated << " => "
225  << std::setw(3) << fPDGListMixed[indx]
226  << " p = " << std::setw(10) << fProb[indx]
227  << " sum_p = " << std::setw(10) << fSumProb[indx] << LOG_END;
228  }
229  }
230 }
231 
232 //____________________________________________________________________________
233 } // namespace flux
234 } // namespace genie
const int kPdgNuE
Definition: PDGCodes.h:28
#define pERROR
Definition: Messenger.h:59
long int Index(void)
returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current ent...
void Clear(Option_t *opt)
reset state variables based on opt
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
A GENIE flux driver using a simple ntuple format.
void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
GSimpleNtpFlux * fGSimpleFlux
ref to avoid repeat dynamic_cast
Definition: GFluxBlender.h:108
const int kPdgAntiNuE
Definition: PDGCodes.h:29
void PrintState(bool verbose=true)
const int kPdgNuMu
Definition: PDGCodes.h:30
double fEnergy
current neutrino&#39;s energy
Definition: GFluxBlender.h:118
double fDistance
current neutrino&#39;s travel distance
Definition: GFluxBlender.h:119
int fPdgCMixed
current neutrino&#39;s new flavor
Definition: GFluxBlender.h:121
GFluxI * AdoptFluxGenerator(GFluxI *generator)
return previous
GFluxI * fRealGFluxI
actual flux generator
Definition: GFluxBlender.h:106
A list of PDG codes.
Definition: PDGCodeList.h:32
int fPdgCGenerated
current neutrino&#39;s original flavor
Definition: GFluxBlender.h:120
PDGCodeList fPDGListMixed
possible flavors after mixing
Definition: GFluxBlender.h:113
std::vector< double > fProb
individual transition probs
Definition: GFluxBlender.h:123
virtual double Probability(int pdg_initial, int pdg_final, double energy, double dist)=0
virtual long int Index(void)=0
returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current ent...
A GENIE flux driver encapsulating the NuMI neutrino flux. It reads-in the official GNUMI neutrino flu...
Definition: GNuMIFlux.h:217
double fRndm
random # used to make choice
Definition: GFluxBlender.h:125
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
TRandom3 generator(0)
GFlavorMixerI * AdoptFlavorMixer(GFlavorMixerI *mixer)
return previous
GENIE interface for flavor modification.
Definition: GFlavorMixerI.h:42
double fBaselineDist
travel dist for mixing (if flux doesn&#39;t support GetDecayDist())
Definition: GFluxBlender.h:116
size_t fNPDGOut
of possible output flavors
Definition: GFluxBlender.h:114
#define pINFO
Definition: Messenger.h:62
GNuMIFlux * fGNuMIFlux
ref to avoid repeat dynamic_cast
Definition: GFluxBlender.h:107
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
virtual void Clear(Option_t *opt)=0
reset state variables based on opt
virtual bool GenerateNext(void)=0
generate the next flux neutrino (return false in err)
double GetDecayDist() const
dist (user units) from dk to current pos
Definition: GNuMIFlux.cxx:395
virtual const PDGCodeList & FluxParticles(void)=0
declare list of flux neutrinos that can be generated (for init. purposes)
const int kPdgNuTau
Definition: PDGCodes.h:32
int ChooseFlavor(int pdg_init, double energy, double dist)
virtual int PdgCode(void)=0
returns the flux neutrino pdg code
virtual const TLorentzVector & Momentum(void)=0
returns the flux neutrino 4-momentum
GFlavorMixerI * fFlavorMixer
flavor modification schema
Definition: GFluxBlender.h:110
const PDGCodeList & FluxParticles(void)
declare list of flux neutrinos that can be generated (for init. purposes)
double GetDecayDist() const
dist (user units) from dk to current pos
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:71
#define LOG_BEGIN(a, b)
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
std::vector< double > fSumProb
cummulative probability
Definition: GFluxBlender.h:124
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
PDGCodeList fPDGListGenerator
possible flavors from generator
Definition: GFluxBlender.h:112
#define LOG_END
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:29