GFluxBlender.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file GFluxBlender.h
3 /// \brief GENIE GFluxI adapter to allow flavor modification
4 ///
5 /// \version $Id: GFluxBlender.cxx,v 1.1.1.1 2010/12/22 16:18:52 p-nusoftart Exp $
6 /// \author Robert Hatcher <rhatcher \at fnal.gov>
7 /// Fermi National Accelerator Laboratory
8 ///
9 /// \update 2010-10-31 initial version
10 ///
11 /// \update 2011-02-22 - JD
12 /// Implemented dummy versions of the new GFluxI::Clear, GFluxI::Index
13 /// and GFluxI::GenerateWeighted methods needed for pre-generation of
14 /// flux interaction probabilities in GMCJDriver.
15 ///
16 ////////////////////////////////////////////////////////////////////////
17 #include <math.h>
18 #include <iostream>
19 #include <iomanip>
20 
21 //GENIE includes
23 #include "Tools/Flux/GNuMIFlux.h"
26 
30 #define LOG_BEGIN(a,b) LOG(a,b)
31 #define LOG_END ""
32 
33 namespace genie {
34 namespace flux {
35 
36 //____________________________________________________________________________
37 
39  GFluxI(),
40  fRealGFluxI(0),
41  fGNuMIFlux(0),
42  fGSimpleFlux(0),
43  fFlavorMixer(0),
44  fPDGListGenerator(),
45  fPDGListMixed(),
46  fNPDGOut(0),
47  fBaselineDist(0),
48  fEnergy(0),
49  fDistance(0),
50  fPdgCGenerated(0),
51  fPdgCMixed(0)
52 { ; }
53 
55 {
56  if ( fRealGFluxI ) { delete fRealGFluxI; fRealGFluxI = 0; }
57  if ( fFlavorMixer ) { delete fFlavorMixer; fFlavorMixer = 0; }
58 }
59 
60 //____________________________________________________________________________
62 {
63  /// Return (reference to) list of possible neutrinos *after* mixing.
64  /// These are PDG code that might be interacted.
65 
66  // this is only ever called by GMCJDriver::GetParticleLists()
67  // during GMCJDriver::Configure() which should happen just once
68  // so don't try to be too clever.
70 
71  // okay, really stupid at this time ...
72  fPDGListMixed.clear();
79 
80  // size the probability arrays to the same number of entries
81  fNPDGOut = fPDGListMixed.size();
82  fProb.resize(fNPDGOut);
83  fSumProb.resize(fNPDGOut);
84 
85  if ( ! fFlavorMixer ) return fRealGFluxI->FluxParticles();
86  else return fPDGListMixed;
87 }
88 
89 //____________________________________________________________________________
91 {
92 
93  bool gen1 = false;
94  while ( ! gen1 ) {
95  if ( ! fRealGFluxI->GenerateNext() ) return false;
96  // have a new entry
98  if ( ! fFlavorMixer ) {
99  // simple case when not configured with a mixing model
101  gen1 = true;
102  } else {
103  // now pick a new flavor
107  fEnergy = fRealGFluxI->Momentum().Energy();
109  // we may have to generate a new neutrino if it oscillates away
110  // don't pass non-particles to GENIE ... it won't like it
111  gen1 = ( fPdgCMixed != 0 );
112  }
113  }
114  return true;
115 }
116 //____________________________________________________________________________
117 void GFluxBlender::Clear(Option_t * opt)
118 {
119 // Clear method needed to conform to GFluxI interface
120 //
121  fRealGFluxI->Clear(opt);
122 }
123 //____________________________________________________________________________
124 long int GFluxBlender::Index(void)
125 {
126 // Index method needed to conform to GFluxI interface
127 //
128  return fRealGFluxI->Index();
129 }
130 //____________________________________________________________________________
131 void GFluxBlender::GenerateWeighted(bool gen_weighted)
132 {
133 // Dummy implementation needed to conform to GFluxI interface
134 //
135  LOG("FluxBlender", pERROR) <<
136  "No GenerateWeighted(bool gen_weighted) method implemented for " <<
137  "gen_weighted: " << gen_weighted;
138 }
139 //____________________________________________________________________________
141 {
142  GFluxI* oldgen = fRealGFluxI;
144  // avoid re-casting
145  fGNuMIFlux = dynamic_cast<GNuMIFlux*>(fRealGFluxI);
146  fGSimpleFlux = dynamic_cast<GSimpleNtpFlux*>(fRealGFluxI);
147  // force evaluation of particle lists
148  this->FluxParticles();
149 
150  return oldgen;
151 }
152 
153 //____________________________________________________________________________
155 {
156  GFlavorMixerI* oldmix = fFlavorMixer;
157  fFlavorMixer = mixer;
158  return oldmix;
159 }
160 
161 //____________________________________________________________________________
162 int GFluxBlender::ChooseFlavor(int pdg_init, double energy, double dist)
163 {
164  // choose a new flavor
165  bool isset = false;
166  int pdg_out = 0;
167  double sumprob = 0;
168 
169  fRndm = RandomGen::Instance()->RndFlux().Rndm();
170  for (size_t indx = 0; indx < fNPDGOut; ++indx ) {
171  int pdg_test = fPDGListMixed[indx];
172  double prob = fFlavorMixer->Probability(pdg_init,pdg_test,energy,dist);
173  fProb[indx] = prob;
174  sumprob += fProb[indx];
175  fSumProb[indx] = sumprob;
176  if ( ! isset && fRndm < sumprob ) {
177  isset = true;
178  pdg_out = pdg_test;
179  }
180  }
181 
182  return pdg_out;
183 }
184 
185 //____________________________________________________________________________
187 {
188  LOG_BEGIN("FluxBlender", pINFO) << "GFluxBlender::PrintConfig()" << LOG_END;
189  if ( fRealGFluxI ) {
190  LOG_BEGIN("FluxBlender", pINFO)
191  << " fRealGFluxI is a \""
192  << typeid(fRealGFluxI).name() << "\""
193  << LOG_END;
194  } else {
195  LOG_BEGIN("FluxBlender", pINFO)
196  << " fRealGFluxI is not initialized" << LOG_END;
197  }
198  if ( fFlavorMixer ) {
199  LOG_BEGIN("FluxBlender", pINFO)
200  << " fFlavorMixer is a \""
201  << typeid(fFlavorMixer).name() << "\""
202  << LOG_END;
203  } else {
204  LOG_BEGIN("FluxBlender", pINFO)
205  << " fFlavorMixer is not initialized" << LOG_END;
206  }
207  LOG_BEGIN("FluxBlender", pINFO)
208  << " BaselineDist " << fBaselineDist << LOG_END;
209  LOG_BEGIN("FluxBlender", pINFO)
210  << "PDG List from Generator" << fPDGListGenerator << LOG_END;
211  LOG_BEGIN("FluxBlender", pINFO)
212  << "PDG List after mixing (n=" << fNPDGOut << ")"
213  << fPDGListMixed << LOG_END;
214 
215 }
216 
217 //____________________________________________________________________________
219 {
220  LOG_BEGIN("FluxBlender", pINFO) << "GFluxBlender::PrintState()" << LOG_END;
221  LOG_BEGIN("FluxBlender", pINFO)
222  << " Flavor " << fPdgCGenerated
223  << " ==> " << fPdgCMixed
224  << " (E=" << fEnergy << ", dist=" << fDistance << ")" << LOG_END;
225  if ( verbose ) {
226  LOG_BEGIN("FluxBlender", pINFO) << " Rndm = " << fRndm << LOG_END;
227  for (size_t indx = 0; indx < fNPDGOut; ++indx ) {
228  LOG_BEGIN("FluxBlender", pINFO)
229  << " [" << indx << "] "
230  << std::setw(3) << fPdgCGenerated << " => "
231  << std::setw(3) << fPDGListMixed[indx]
232  << " p = " << std::setw(10) << fProb[indx]
233  << " sum_p = " << std::setw(10) << fSumProb[indx] << LOG_END;
234  }
235  }
236 }
237 
238 //____________________________________________________________________________
239 } // namespace flux
240 } // namespace genie
const XML_Char * name
Definition: expat.h:151
const int kPdgNuE
Definition: PDGCodes.h:28
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
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:79
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:100
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:110
double fDistance
current neutrino&#39;s travel distance
Definition: GFluxBlender.h:111
int fPdgCMixed
current neutrino&#39;s new flavor
Definition: GFluxBlender.h:113
Loaders::FluxType flux
GFluxI * AdoptFluxGenerator(GFluxI *generator)
return previous
GFluxI * fRealGFluxI
actual flux generator
Definition: GFluxBlender.h:98
A list of PDG codes.
Definition: PDGCodeList.h:33
int fPdgCGenerated
current neutrino&#39;s original flavor
Definition: GFluxBlender.h:112
double dist
Definition: runWimpSim.h:113
PDGCodeList fPDGListMixed
possible flavors after mixing
Definition: GFluxBlender.h:105
std::vector< double > fProb
individual transition probs
Definition: GFluxBlender.h:115
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
GENIE GFluxI adapter to allow flavor modification.
double fRndm
random # used to make choice
Definition: GFluxBlender.h:117
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
GFlavorMixerI * AdoptFlavorMixer(GFlavorMixerI *mixer)
return previous
GENIE interface for flavor modification.
Definition: GFlavorMixerI.h:37
#define LOG_BEGIN(a, b)
double energy
Definition: plottest35.C:25
double fBaselineDist
travel dist for mixing (if flux doesn&#39;t support GetDecayDist())
Definition: GFluxBlender.h:108
size_t fNPDGOut
of possible output flavors
Definition: GFluxBlender.h:106
#define pINFO
Definition: Messenger.h:63
GNuMIFlux * fGNuMIFlux
ref to avoid repeat dynamic_cast
Definition: GFluxBlender.h:99
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:551
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:102
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:72
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:116
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
PDGCodeList fPDGListGenerator
possible flavors from generator
Definition: GFluxBlender.h:104
#define LOG_END
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37