test_gnumi.C
Go to the documentation of this file.
1 #include "TMath.h"
2 #include "TSystem.h"
3 #include <string>
4 #include <fstream>
5 #include <iostream>
6 #include <iomanip>
7 using namespace std;
8 
9 #include "TStyle.h"
10 #include "TCanvas.h"
11 #include "TH1D.h"
12 #include "TLorentzVector.h"
13 
14 #include "FluxDrivers/GNuMIFlux.h"
15 #include "PDG/PDGCodeList.h"
16 #include "PDG/PDGUtils.h"
17 
18 using namespace genie::flux;
19 
20 #include "FluxDrivers/GNuMIFlux.h"
21 #include "PDG/PDGCodeList.h"
22 #include "PDG/PDGUtils.h"
23 
24 // typedefs and forward declarations
26 
27 double fabserrX(double a, double b)
28 { return TMath::Abs(a-b)/TMath::Max(TMath::Abs(b),1.0e-30); }
29 
30 //////////////////////////////////////////////////////////////////////////////
31 
32 int calc_enuwgt(const fluxentry_t& fluxentry,
33  double x, double y, double z,
34  double& enu, double& wgt_xy, xypartials& partials);
35 
36 //////////////////////////////////////////////////////////////////////////////
37 
38 void test_gnumi(int mxpull = 0)
39 {
40  // mxpull: max number of entries to process 0=all, -N=until N mismatches
41 
42  // currently TFile not TChain, so better resolve to a single file
43  string fluxpatt = "";
44  fluxpatt = "$FLUXPATH/v19/fluka05_le010z185i/root/fluka05_le010z185i_1.root";
45  //fluxpatt = "$FLUXPATH/g4numi/g4fluka_le010z185i_leak_v10_0001.root";
46  //fluxpatt = "$GENIE/src/FluxDrivers/gnumi_test/*g3*.root";
47  //fluxpatt = "*g3*.root";
48  string fluxfilename(gSystem->ExpandPathName(fluxpatt.c_str()));
49 
51  gnumi->LoadBeamSimData(fluxfilename);
52  gnumi->SetGenWeighted(true); // generated entries w/ fWeight != 1.0
53 
54  const genie::PDGCodeList& pdgList = gnumi->FluxParticles();
55  int nfluxtypes = pdgList.size();
56  cout << "GNuMIFlux supplies " << nfluxtypes << " PDG particles of type: ";
57  // don't try iterators .. apparently they aren't usable in CINT
58  for (int indx=0; indx < nfluxtypes ; ++indx)
59  cout << pdgList[indx] << " ";
60  cout << endl;
61 
62  cout << "MaxEnergy " << gnumi->MaxEnergy() << endl;
63 
64  cout << "scan for max weight at NearDet" << endl;
65  //gnumi->UseFluxAtNearDetCenter();
66  gnumi->SetLengthUnits(genie::flux::GNuMIFlux::kcm);
67  gnumi->SetBeamFluxWindow(genie::flux::GNuMIFlux::kMinosNearCenter);
68  gnumi->ScanForMaxWeight();
69 
70  int npull, nxy, file_debug;
71  double xbeam, ybeam, zbeam = 103935.; // distance in cm to near "window"
72 
73  ifstream xypoints("XYPOINTS.TXT");
74 
75  xypoints >> npull >> nxy >> zbeam >> file_debug;
76 
77  TH1D* hist_ferr_enu = new TH1D("ferr_enu","ferr_enu",1000,0.0,0.000005);
78  TH1D* hist_ferr_wgt = new TH1D("ferr_wgt","ferr_wgt",1000,0.0,0.01);
79  hist_ferr_enu->Sumw2();
80  hist_ferr_wgt->Sumw2();
81  TH1D* hist_ferr_enu0 = new TH1D("ferr_enu0","ferr_enu0",1000,0.0,0.000005);
82  TH1D* hist_ferr_wgt0 = new TH1D("ferr_wgt0","ferr_wgt0",1000,0.0,0.01);
83  hist_ferr_enu0->Sumw2();
84  hist_ferr_wgt0->Sumw2();
85  gStyle->SetOptStat(111111);
86 
87  if ( mxpull > 0 && npull > mxpull ) npull = mxpull;
88 
89  for ( int ipull = 1; ipull <= npull ; ++ipull ) {
90 
91  bool genok = gnumi->GenerateNext();
92  const fluxentry_t& fluxentry = gnumi->PassThroughInfo();
93 
94  double wgtprd_gen = gnumi->Weight();
95  double enu_gen = gnumi->Momentum().E();
96  if ( wgtprd_gen != 1.0 ) {
97  // weighted flux ... we can validate "center" values
98  // assume used kMinosNearCenter in SetBeamFluxWindow
99  double wgtprd_ntuple = fluxentry.nimpwt * fluxentry.nwtnear;
100  double enu_ntuple = fluxentry.nenergyn;
101  double ferr_enu0 = fabserrX(enu_gen,enu_ntuple);
102  double ferr_wgtprd0 = fabserrX(wgtprd_gen,wgtprd_ntuple);
103  const double eps0h = 1.0e-6; // 2.5e-5;
104  if ( ferr_enu0 > eps0h || ferr_wgtprd0 > eps0h ) {
105  hist_ferr_enu0->Fill(ferr_enu0,1.0);
106  hist_ferr_wgt0->Fill(ferr_wgtprd0,1.0);
107  }
108  const double eps0 = 2.5e-5; // 2.5e-5;
109  if ( ferr_enu0 > eps0 || ferr_wgtprd0 > eps0 ) {
110  char xenu0 = ( ferr_enu0 > eps0 ) ? '#' : ' ';
111  char xwgt0 = ( ferr_wgtprd0 > eps0 ) ? '#' : ' ';
112  cout << "Checking \"center\" energy/weight "
113  << " enu " << setw(8) << enu_gen << " " << setw(8) << enu_ntuple
114  << " err " << setw(11) << ferr_enu0 << xenu0
115  << " wgtprd " << setw(8) << wgtprd_gen << " " << setw(8) << wgtprd_ntuple
116  << " err " << setw(11) << ferr_wgtprd0 << xwgt0
117  << endl;
118  }
119 
120  }
121 
122 
123  //bool atend = gnumi->End();
124  //cout << "gnumi->GenerateNext() returned " << (genok?"true":"false")
125  // << " End() returns " << (atend?"true":"false")
126  // << endl;
127 
128  //cout << "Current entry: "
129  // << " PdgCode=" << gnumi->PdgCode()
130  // << " Weight=" << gnumi->Weight()
131  // << endl;
132  // gnumi->PrintCurrent();
133 
134  for ( int ixy = 0; ixy <= nxy ; ++ixy ) {
135  xypartials part_file, part_calc;
136  if ( file_debug ) part_file.ReadStream(xypoints);
137 
138  double enu, wgt_xy;
139  double enu_in, wgt_xy_in;
140  int ipull_in, ixy_in;
141  xypoints >> ipull_in >> ixy_in >> xbeam >> ybeam >> enu_in >> wgt_xy_in;
142  if ( ! xypoints.good() || ( ipull_in != ipull ) || ( ixy_in != ixy) ) {
143  cout << "HEY!!! XYPOINTS.TXT "
144  << " file status = " << (int)xypoints.good()
145  << " ipull " << ipull_in << " " << ipull
146  << " ixy " << ixy_in << " " << ixy
147  << endl;
148  return;
149  }
150  fluxentry.CalcEnuWgt(xbeam,ybeam,zbeam,enu,wgt_xy,part_calc);
151 
152  double ferr_enu = fabserrX(enu,enu_in);
153  double ferr_wgt_xy = fabserrX(wgt_xy,wgt_xy_in);
154  const double eps_enu = 2.5e-5; // 0.003; // 6.0e-4; // 2.5e-4; // 1.5e-5; // 1.0e-6;
155  const double eps_wgt = 2.5e-5; // 0.003; // 6.0e-4; // 2.5e-4; // 2.5e-5; //
156  if ( ferr_enu > eps_enu || ferr_wgt_xy > eps_wgt ) {
157  cout << endl;
158  char xenu = ( ferr_enu > eps_enu ) ? '#':' ';
159  char xwgt = ( ferr_wgt_xy > eps_wgt ) ? '#':' ';
160  cout << "pull " << setw(6) << ipull << " ixy " << setw(2) << ixy
161  << " enu " << setw(8) << enu << " " << setw(8) << enu_in
162  << " err " << setw(11) << ferr_enu << xenu
163  << " wgt_xy " << setw(11) << wgt_xy << " " << setw(11) << wgt_xy_in
164  << " err " << setw(11) << ferr_wgt_xy << xwgt
165  << " ptype " << setw(4) << part_file.ptype
166  << " ntype " << setw(3) << part_file.ntype << endl;
167  hist_ferr_enu->Fill(ferr_enu,1.0);
168  hist_ferr_wgt->Fill(ferr_wgt_xy,1.0);
169  }
170  //else {
171  // cout << "pull " << ipull << " ixy " << ixy
172  // << " okay ptype " << part_file.ptype << " ntype " << part_file.ntype << endl;
173  //}
174 
175  if (file_debug) {
176  int np = part_file.Compare(part_calc);
177  if (np>0) {
178  std::cout << fluxentry << endl;
179  part_file.Print();
180  part_calc.Print();
181 
182  if ( mxpull < 0 ) {
183  cout << "COUNTDOWN mismatch error count @ " << mxpull << endl;
184  mxpull += 1;
185  if ( mxpull == 0 ) {
186  cout << "REACHED LIMIT ON MISMATCHES ... exit" << endl;
187  return; // reached limit on errors, asked to stop on error
188  }
189  }
190  }
191  }
192 
193 
194  //if ( part_calc.ptype == 13 || part_calc.ptype == -13 ) {
195  // cout << "STOPPING HERE FOR INSPECTION" << endl;
196  // return;
197  //}
198 
199  } //loop over xy points for given entryƧ
200  } // loop over entries (pulls)
201 
202  TCanvas* c1 = new TCanvas("c1","all positions");
203  c1->Divide(1,2);
204  c1->cd(1);
205  gPad->SetLogy(true);
206  hist_ferr_enu->Draw();
207  c1->cd(2);
208  gPad->SetLogy(true);
209  hist_ferr_wgt->Draw();
210 
211  TCanvas* c2 = new TCanvas("c2","center weight");
212  c2->Divide(1,2);
213  c2->cd(1);
214  gPad->SetLogy(true);
215  hist_ferr_enu0->Draw();
216  c2->cd(2);
217  gPad->SetLogy(true);
218  hist_ferr_wgt0->Draw();
219 
220 }
221 
222 //////////////////////////////////////////////////////////////////////////////
223 int calc_enuwgt(const fluxentry_t& fluxentry,
224  double xpos, double ypos, double zpos,
225  double& enu, double& wgt_xy, xypartials& partials)
226 {
227 
228  // Neutrino Energy and Weigth at arbitrary point
229  // based on:
230  // NuMI-NOTE-BEAM-0109 (MINOS DocDB # 109)
231  // Title: Neutrino Beam Simulation using PAW with Weighted Monte Carlos
232  // Author: Rick Milburn
233  // Date: 1995-10-01
234 
235  // history:
236  // jzh 3/21/96 grab R.H.Milburn's weighing routine
237  // jzh 5/ 9/96 substantially modify the weighting function use dot product
238  // instead of rotation vecs to get theta get all info except
239  // det from ADAMO banks neutrino parent is in Particle.inc
240  // Add weighting factor for polarized muon decay
241  // jzh 4/17/97 convert more code to double precision because of problems
242  // with Enu>30 GeV
243  // rwh 10/ 9/08 transliterate function from f77 to C++
244 
245  // original function description:
246  // Real function for use with PAW Ntuple To transform from destination
247  // detector geometry to the unit sphere moving with decaying hadron with
248  // velocity v, BETA=v/c, etc.. For (pseudo)scalar hadrons the decays will
249  // be isotropic in this sphere so the fractional area (out of 4-pi) is the
250  // fraction of decays that hit the target. For a given target point and
251  // area, and given x-y components of decay transverse location and slope,
252  // and given decay distance from target ans given decay GAMMA and
253  // rest-frame neutrino energy, the lab energy at the target and the
254  // fractional solid angle in the rest-frame are determined.
255  // For muon decays, correction for non-isotropic nature of decay is done.
256 
257  // Arguments:
258  // genie::flux::GNuMIFluxPassThroughInfo fluxentry :: ntuple entry info
259  // double x, y, z :: position to evaluate for (enu,wgt_xy)
260  // in *beam* frame coordinates (?units)
261  // double enu, wgt_xy :: resulting energy and weight
262  // Return:
263  // int :: error code
264  // Assumptions:
265  // Energies given in GeV
266  // Particle codes have been translated from GEANT into PDG codes
267 
268  // for now ... these _should_ come from DB
269  // but use these hard-coded values to "exactly" reproduces old code
270  const double kPIMASS = 0.13957;
271  const double kKMASS = 0.49368;
272  const double kK0MASS = 0.49767;
273  const double kMUMASS = 0.105658389;
274 
275  const int kpdg_nue = 12; // extended Geant 53
276  const int kpdg_nuebar = -12; // extended Geant 52
277  const int kpdg_numu = 14; // extended Geant 56
278  const int kpdg_numubar = -14; // extended Geant 55
279 
280  const int kpdg_muplus = -13; // Geant 5
281  const int kpdg_muminus = 13; // Geant 6
282  const int kpdg_pionplus = 211; // Geant 8
283  const int kpdg_pionminus = -211; // Geant 9
284  const int kpdg_k0long = 130; // Geant 10 ( K0=311, K0S=310 )
285  const int kpdg_k0short = 310; // Geant 16
286  const int kpdg_k0mix = 311;
287  const int kpdg_kaonplus = 321; // Geant 11
288  const int kpdg_kaonminus = -321; // Geant 12
289 
290  const double kRDET = 100.0; // set to flux per 100 cm radius
291 
292  enu = 0.0; // don't know what the final value is
293  wgt_xy = 0.0; // but set these in case we return early due to error
294 
295 
296  // in principle we should get these from the particle DB
297  // but for consistency testing use the hardcoded values
298  double parent_mass = kPIMASS;
299  switch ( fluxentry.ptype ) {
300  case kpdg_pionplus:
301  case kpdg_pionminus:
302  parent_mass = kPIMASS;
303  break;
304  case kpdg_kaonplus:
305  case kpdg_kaonminus:
306  parent_mass = kKMASS;
307  break;
308  case kpdg_k0long:
309  case kpdg_k0short:
310  case kpdg_k0mix:
311  parent_mass = kK0MASS;
312  break;
313  case kpdg_muplus:
314  case kpdg_muminus:
315  parent_mass = kMUMASS;
316  break;
317  default:
318  std::cerr << "NU_REWGT unknown particle type " << fluxentry.ptype
319  << std::endl;
320  return 1;
321  }
322 
323  double parentp2 = ( fluxentry.pdpx*fluxentry.pdpx +
324  fluxentry.pdpy*fluxentry.pdpy +
325  fluxentry.pdpz*fluxentry.pdpz );
326  double parent_energy = TMath::Sqrt( parentp2 +
327  parent_mass*parent_mass);
328  double parentp = TMath::Sqrt( parentp2 );
329 
330  double gamma = parent_energy / parent_mass;
331  double gamma_sqr = gamma * gamma;
332  double beta_mag = TMath::Sqrt( ( gamma_sqr - 1.0 )/gamma_sqr );
333 
334  // Get the neutrino energy in the parent decay CM
335  double enuzr = fluxentry.necm;
336  // Get angle from parent line of flight to chosen point in beam frame
337  double rad = TMath::Sqrt( (xpos-fluxentry.vx)*(xpos-fluxentry.vx) +
338  (ypos-fluxentry.vy)*(ypos-fluxentry.vy) +
339  (zpos-fluxentry.vz)*(zpos-fluxentry.vz) );
340 
341  double costh_pardet = ( fluxentry.pdpx*(xpos-fluxentry.vx) +
342  fluxentry.pdpy*(ypos-fluxentry.vy) +
343  fluxentry.pdpz*(zpos-fluxentry.vz) )
344  / ( parentp * rad);
345  if ( costh_pardet > 1.0 ) costh_pardet = 1.0;
346  if ( costh_pardet < -1.0 ) costh_pardet = -1.0;
347  double theta_pardet = TMath::ACos(costh_pardet);
348 
349  // Weighted neutrino energy in beam, approx, good for small theta
350  //RWH//double emrat = 1.0 / ( gamma * ( 1.0 - beta_mag * TMath::Cos(theta_pardet)));
351  double emrat = 1.0 / ( gamma * ( 1.0 - beta_mag * costh_pardet ));
352  enu = emrat * enuzr; // ! the energy ... normally
353 
354  // RWH-debug
355  bool debug = false;
356  if (debug) {
357  cout << setprecision(15);
358  cout << "ptype " << fluxentry.ptype << " m " << parent_mass
359  << " p " << parentp << " e " << parent_energy << " gamma " << gamma
360  << " beta " << beta_mag << endl;
361 
362  cout << " enuzr " << enuzr << " rad " << rad << " costh " << costh_pardet
363  << " theta " << theta_pardet << " emrat " << emrat
364  << " enu " << enu << endl;
365  }
366  partials.parent_mass = parent_mass;
367  partials.parentp = parentp;
368  partials.parent_energy = parent_energy;
369  partials.gamma = gamma;
370  partials.beta_mag = beta_mag;
371  partials.enuzr = enuzr;
372  partials.rad = rad;
373  partials.costh_pardet = costh_pardet;
374  partials.theta_pardet = theta_pardet;
375  partials.emrat = emrat;
376  partials.eneu = enu;
377 
378  // Get solid angle/4pi for detector element
379  double sangdet = ( kRDET*kRDET / ( (zpos-fluxentry.vz)*(zpos-fluxentry.vz)))/4.0;
380 
381  // Weight for solid angle and lorentz boost
382  wgt_xy = sangdet * ( emrat * emrat ); // ! the weight ... normally
383 
384  partials.sangdet = sangdet;
385  partials.wgt = wgt_xy;
386  partials.ptype = fluxentry.ptype; // assume already PDG
387 
388  // Done for all except polarized muon decay
389  // in which case need to modify weight
390  // (must be done in double precision)
391  if ( fluxentry.ptype == kpdg_muplus || fluxentry.ptype == kpdg_muminus) {
392  double beta[3], p_dcm_nu[4], p_nu[3], p_pcm_mp[3], partial;
393 
394  // Boost neu neutrino to mu decay CM
395  beta[0] = fluxentry.pdpx / parent_energy;
396  beta[1] = fluxentry.pdpy / parent_energy;
397  beta[2] = fluxentry.pdpz / parent_energy;
398  p_nu[0] = (xpos-fluxentry.vx)*enu/rad;
399  p_nu[1] = (ypos-fluxentry.vy)*enu/rad;
400  p_nu[2] = (zpos-fluxentry.vz)*enu/rad;
401  partial = gamma *
402  (beta[0]*p_nu[0] + beta[1]*p_nu[1] + beta[2]*p_nu[2] );
403  partial = enu - partial/(gamma+1.0);
404  p_dcm_nu[0] = p_nu[0] - beta[0]*gamma*partial;
405  p_dcm_nu[1] = p_nu[1] - beta[1]*gamma*partial;
406  p_dcm_nu[2] = p_nu[2] - beta[2]*gamma*partial;
407  p_dcm_nu[3] = TMath::Sqrt( p_dcm_nu[0]*p_dcm_nu[0] +
408  p_dcm_nu[1]*p_dcm_nu[1] +
409  p_dcm_nu[2]*p_dcm_nu[2] );
410 
411  partials.betanu[0] = beta[0];
412  partials.betanu[1] = beta[1];
413  partials.betanu[2] = beta[2];
414  partials.p_nu[0] = p_nu[0];
415  partials.p_nu[1] = p_nu[1];
416  partials.p_nu[2] = p_nu[2];
417  partials.partial1 = partial;
418  partials.p_dcm_nu[0] = p_dcm_nu[0];
419  partials.p_dcm_nu[1] = p_dcm_nu[1];
420  partials.p_dcm_nu[2] = p_dcm_nu[2];
421  partials.p_dcm_nu[3] = p_dcm_nu[3];
422 
423  // Boost parent of mu to mu production CM
424  double particle_energy = fluxentry.ppenergy;
425  gamma = particle_energy/parent_mass;
426  beta[0] = fluxentry.ppdxdz * fluxentry.pppz / particle_energy;
427  beta[1] = fluxentry.ppdydz * fluxentry.pppz / particle_energy;
428  beta[2] = fluxentry.pppz / particle_energy;
429  partial = gamma * ( beta[0]*fluxentry.muparpx +
430  beta[1]*fluxentry.muparpy +
431  beta[2]*fluxentry.muparpz );
432  partial = fluxentry.mupare - partial/(gamma+1.0);
433  p_pcm_mp[0] = fluxentry.muparpx - beta[0]*gamma*partial;
434  p_pcm_mp[1] = fluxentry.muparpy - beta[1]*gamma*partial;
435  p_pcm_mp[2] = fluxentry.muparpz - beta[2]*gamma*partial;
436  double p_pcm = TMath::Sqrt ( p_pcm_mp[0]*p_pcm_mp[0] +
437  p_pcm_mp[1]*p_pcm_mp[1] +
438  p_pcm_mp[2]*p_pcm_mp[2] );
439 
440  //cout << " muparpxyz " << fluxentry.muparpx << " "
441  // << fluxentry.muparpy << " " << fluxentry.muparpz << endl;
442  //cout << " beta " << beta[0] << " " << beta[1] << " " << beta[2] << endl;
443  //cout << " gamma " << gamma << " partial " << partial << endl;
444  //cout << " p_pcm_mp " << p_pcm_mp[0] << " " << p_pcm_mp[1] << " " << p_pcm_mp[2] << " " << p_pcm << endl;
445  partials.muparent_px = fluxentry.muparpx;
446  partials.muparent_py = fluxentry.muparpy;
447  partials.muparent_pz = fluxentry.muparpz;
448  partials.gammamp = gamma;
449  partials.betamp[0] = beta[0];
450  partials.betamp[1] = beta[1];
451  partials.betamp[2] = beta[2];
452  partials.partial2 = partial;
453  partials.p_pcm_mp[0] = p_pcm_mp[0];
454  partials.p_pcm_mp[1] = p_pcm_mp[1];
455  partials.p_pcm_mp[2] = p_pcm_mp[2];
456  partials.p_pcm = p_pcm;
457 
458  const double eps = 1.0e-30; // ? what value to use
459  if ( p_pcm < eps || p_dcm_nu[3] < eps ) {
460  return 3; // mu missing parent info?
461  }
462  // Calc new decay angle w.r.t. (anti)spin direction
463  double costh = ( p_dcm_nu[0]*p_pcm_mp[0] +
464  p_dcm_nu[1]*p_pcm_mp[1] +
465  p_dcm_nu[2]*p_pcm_mp[2] ) /
466  ( p_dcm_nu[3]*p_pcm );
467  if ( costh > 1.0 ) costh = 1.0;
468  if ( costh < -1.0 ) costh = -1.0;
469  // Calc relative weight due to angle difference
470  double wgt_ratio = 0.0;
471  switch ( fluxentry.ntype ) {
472  case kpdg_nue:
473  case kpdg_nuebar:
474  wgt_ratio = 1.0 - costh;
475  break;
476  case kpdg_numu:
477  case kpdg_numubar:
478  double xnu = 2.0 * enuzr / kMUMASS;
479  wgt_ratio = ( (3.0-2.0*xnu ) - (1.0-2.0*xnu)*costh ) / (3.0-2.0*xnu);
480  break;
481  default:
482  return 2; // bad neutrino type
483  }
484  wgt_xy = wgt_xy * wgt_ratio;
485 
486  partials.ntype = fluxentry.ntype; // assume converted to PDG
487  partials.costhmu = costh;
488  partials.wgt_ratio = wgt_ratio;
489 
490  }
491 
492  return 0;
493 }
void SetLengthUnits(double user_units)
Set units assumed by user.
Definition: GNuMIFlux.cxx:1362
double fabserrX(double a, double b)
Definition: test_gnumi.C:27
OStream cerr
Definition: OStream.cxx:7
Double_t beta
Definition: lz4.cxx:387
A list of PDG codes.
Definition: PDGCodeList.h:33
c2
Definition: demo5.py:33
GENIE flux drivers.
A GENIE flux driver encapsulating the NuMI neutrino flux. It reads-in the official GNUMI neutrino flu...
Definition: GNuMIFlux.h:217
double Weight(void)
returns the flux neutrino weight (if any)
Definition: GNuMIFlux.h:234
double enu
Definition: runWimpSim.h:113
double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
Definition: GNuMIFlux.h:231
virtual void LoadBeamSimData(const std::vector< std::string > &filenames, const std::string &det_loc)
Definition: GNuMIFlux.cxx:640
genie::flux::GNuMIFluxPassThroughInfo fluxentry_t
Definition: test_gnumi.C:25
const GNuMIFluxPassThroughInfo & PassThroughInfo(void)
GNuMIFluxPassThroughInfo.
Definition: GNuMIFlux.h:252
const double a
int CalcEnuWgt(const TLorentzVector &xyz, double &enu, double &wgt_xy) const
Definition: GNuMIFlux.cxx:1885
const PDGCodeList & FluxParticles(void)
declare list of flux neutrinos that can be generated (for init. purposes)
Definition: GNuMIFlux.h:230
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
Definition: GNuMIFlux.cxx:295
static constexpr Double_t rad
Definition: Munits.h:162
z
Definition: test.py:28
OStream cout
Definition: OStream.cxx:6
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
const TLorentzVector & Momentum(void)
returns the flux neutrino 4-momentum
Definition: GNuMIFlux.h:235
const hit & b
Definition: hits.cxx:21
c1
Definition: demo5.py:24
void ScanForMaxWeight(void)
scan for max flux weight (before generating unweighted flux neutrinos)
Definition: GNuMIFlux.cxx:836
Float_t e
Definition: plot.C:35
void test_gnumi(int mxpull=0)
Definition: test_gnumi.C:38
int calc_enuwgt(const fluxentry_t &fluxentry, double x, double y, double z, double &enu, double &wgt_xy, xypartials &partials)
Definition: test_gnumi.C:223
void SetGenWeighted(bool genwgt=false)
toggle whether GenerateNext() returns weight=1 flux (initial default false)
Definition: GNuMIFlux.h:292