GDk2NuFlux.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2012, GENIE Neutrino MC Generator Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Robert Hatcher <rhatcher@fnal.gov>
8  Fermi National Accelerator Laboratory
9 
10  For the class documentation see the corresponding header file.
11 
12 
13 */
14 //____________________________________________________________________________
15 
16 //#define __GENIE_LOW_LEVEL_MESG_ENABLED__
17 
18 #include <cstdlib>
19 #include <fstream>
20 #include <vector>
21 #include <sstream>
22 #include <cassert>
23 #include <climits>
24 
25 #include "libxml/xmlmemory.h"
26 #include "libxml/parser.h"
27 
28 // ROOT headers
29 #include <TFile.h>
30 #include <TChain.h>
31 #include <TChainElement.h>
32 #include <TSystem.h>
33 #include <TStopwatch.h>
34 
35 #ifdef GENIE_PRE_R3
36  // GENIE headers
37  #include "Conventions/GBuild.h"
38  #include "Conventions/Units.h"
39  #include "Messenger/Messenger.h"
40  #include "Numerical/RandomGen.h"
41  #include "PDG/PDGCodes.h"
42  #include "PDG/PDGCodeList.h"
43  #include "Utils/MathUtils.h"
44  #include "Utils/PrintUtils.h"
45  #include "Utils/UnitUtils.h"
46  #include "Utils/XmlParserUtils.h"
47  #include "Utils/StringUtils.h"
48 #else
49  // GENIE R-3 headers
61 #endif
62 
63 // Dk2Nu headers
64 #include "genie/GDk2NuFlux.h"
65 #include "tree/dk2nu.h"
66 #include "tree/dkmeta.h"
67 #include "tree/NuChoice.h"
68 #include "tree/calcLocationWeights.h"
69 
70 // need GDk2NuFlux header to register w/ factory
71 #ifdef GENIE_PRE_R3
72  #include "Conventions/GVersion.h"
73  #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,9,0)
74  #include "FluxDrivers/GFluxDriverFactory.h"
76  #endif
77 #else
80 #endif
81 
82 #include <vector>
83 #include <algorithm>
84 #include <iomanip>
85 #include "TRegexp.h"
86 #include "TString.h"
87 
88 using namespace genie;
89 using namespace genie::flux;
90 
91 // declaration of helper class
92 namespace genie {
93  namespace flux {
94  class GDk2NuFluxXMLHelper {
95  public:
97  : fVerbose(0), fGDk2NuFlux(dk2nuFlux) { ; }
99  bool LoadConfig(std::string cfg);
100 
101  // these should go in a more general package
102  std::vector<double> GetDoubleVector(std::string str);
103  std::vector<long int> GetIntVector(std::string str);
104 
105  private:
106  bool LoadParamSet(xmlDocPtr&, std::string cfg);
107  void ParseParamSet(xmlDocPtr&, xmlNodePtr&);
108  void ParseBeamDir(xmlDocPtr&, xmlNodePtr&);
109  void ParseBeamPos(std::string);
110  void ParseRotSeries(xmlDocPtr&, xmlNodePtr&);
111  void ParseWindowSeries(xmlDocPtr&, xmlNodePtr&);
112  void ParseEnuMax(std::string);
113  void ParseMaxWgtFail(std::string);
114  TVector3 AnglesToAxis(double theta, double phi, std::string units = "deg");
115  TVector3 ParseTV3(const std::string& );
116 
117  int fVerbose; ///< how noisy to be when parsing XML
118  // who to apply these changes to
119  GDk2NuFlux* fGDk2NuFlux;
120 
121  // derived offsets/rotations
122  TVector3 fBeamPosXML;
123  TRotation fBeamRotXML;
124  TVector3 fFluxWindowPtXML[3];
125  };
126  }
127 }
128 
129 //____________________________________________________________________________
131 // #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,9,0)
133 // #endif
134 {
135  this->Initialize();
136 }
137 //___________________________________________________________________________
139 {
140  this->CleanUp();
141 }
142 //___________________________________________________________________________
143 // #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,9,0)
144 double GDk2NuFlux::GetTotalExposure() const
145 {
146  // complete the GFluxExposureI interface
147  return UsedPOTs();
148 }
149 // #endif
150 //___________________________________________________________________________
151 int GDk2NuFlux::PdgCode(void)
152 {
153  return fCurNuChoice->pdgNu;
154 }
155 //___________________________________________________________________________
156 const TLorentzVector& GDk2NuFlux::Momentum(void)
157 {
158  return fCurNuChoice->p4NuUser;
159 }
160 //___________________________________________________________________________
161 const TLorentzVector& GDk2NuFlux::Position(void)
162 {
163  return fCurNuChoice->x4NuUser;
164 }
165 //___________________________________________________________________________
166 bool GDk2NuFlux::GenerateNext(void)
167 {
168 // Get next (unweighted) flux ntuple entry on the specified detector location
169 //
171  while ( true ) {
172  // Check for end of flux ntuple
173  bool end = this->End();
174  if ( end ) {
175  LOG("Flux", pNOTICE) << "GenerateNext signaled End() ";
176  return false;
177  }
178 
179  // Get next weighted flux ntuple entry
180  bool nextok = this->GenerateNext_weighted();
181  if ( fGenWeighted ) return nextok;
182  if ( ! nextok ) continue;
183 
184  /* RWH - debug purposes
185  if ( fNCycles == 0 ) {
186  LOG("Flux", pNOTICE)
187  << "Got flux entry: " << fIEntry
188  << " - Cycle: "<< fICycle << "/ infinite";
189  } else {
190  LOG("Flux", pNOTICE)
191  << "Got flux entry: "<< fIEntry
192  << " - Cycle: "<< fICycle << "/"<< fNCycles;
193  }
194  */
195 
196  // Get fractional weight & decide whether to accept curr flux neutrino
197  double f = this->Weight() / fMaxWeight;
198  //LOG("Flux", pNOTICE)
199  // << "Curr flux neutrino fractional weight = " << f;
200  if (f > 1.) {
201  // Oh,dear me ... what the weight max _should_ have been
202  double bumped = this->Weight() * fMaxWgtFudge;
203  if ( bumped > fMaxWeightMax ) fMaxWeightMax = bumped;
204  ++fMaxWgtExceeded;
205 
206  // what to do ...
207  if ( fMaxWgtFailModel <= 0 ) {
208  fMaxWeight = bumped;
209  LOG("Flux", pERROR)
210  << "** Fractional weight = " << f
211  << " > 1 !! Bump fMaxWeight estimate to " << fMaxWeight << "\n"
212  << fCurDk2Nu->AsString() << "\n" << fCurNuChoice->AsString();
213  std::cout << std::flush;
214  } else {
215  LOG("Flux", pERROR)
216  << "** Fractional weight = " << f
217  << " > 1 !! Leave fMaxWeight frozen at " << fMaxWeight << " "
218  << " not bumped to " << bumped << ", max so far: "
219  << fMaxWeightMax << "\n"
220  << fCurDk2Nu->AsString() << "\n" << fCurNuChoice->AsString();
221  std::cout << std::flush;
222  if ( fMaxWgtFailModel >= 2 ) {
223  PrintConfig();
224  std::cout << std::flush;
225  LOG("Flux", pFATAL)
226  << "User requested not to continue";
227  abort();
228  }
229  }
230  }
231  double r = rnd->RndFlux().Rndm(); // (f < 1.) ? rnd->RndFlux().Rndm() : 0;
232  bool accept = ( r < f );
233  if ( accept ) {
234 
235 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
236  LOG("Flux", pNOTICE)
237  << "Generated beam neutrino: "
238  << "\n pdg-code: " << fCurNuChoice->pdgNu
239  << "\n p4: " << utils::print::P4AsShortString(&(fCurNuChoice->p4NuBeam))
240  << "\n x4: " << utils::print::X4AsString(&(fCurNuChoice->x4NuBeam))
241  << "\n p4: " << utils::print::P4AsShortString(&(fCurNuChoice->p4NuUser))
242  << "\n x4: " << utils::print::X4AsString(&(fCurNuChoice->x4NuUser));
243 #endif
244 
245  fWeight = 1.;
246  return true;
247  }
248 
249  //LOG("Flux", pNOTICE)
250  // << "** Rejecting current flux neutrino based on the flux weight only";
251  }
252  return false;
253 }
254 //___________________________________________________________________________
256 {
257 // Get next (weighted) flux ntuple entry on the specified detector location
258 //
259 
260  // Check whether a flux ntuple has been loaded
261  if ( ! fNuFluxTree ) {
262  LOG("Flux", pERROR)
263  << "The flux driver has not been properly configured";
264  return false;
265  }
266 
267  // Reuse an entry?
268  //std::cout << " ***** iuse " << fIUse << " nuse " << fNUse
269  // << " ientry " << fIEntry << " nentry " << fNEntries
270  // << " icycle " << fICycle << " ncycle " << fNCycles << std::endl;
271  if ( fIUse < fNUse && fIEntry >= 0 ) {
272  // Reuse this entry
273  fIUse++;
274  } else {
275  // Reset previously generated neutrino code / 4-p / 4-x
276  this->ResetCurrent();
277  // Move on, read next flux ntuple entry
278  fIEntry++;
279  if ( fIEntry >= fNEntries ) {
280  // Ran out of entries @ the current cycle of this flux file
281  // Check whether more (or infinite) number of cycles is requested
282  if (fICycle < fNCycles || fNCycles == 0 ) {
283  fICycle++;
284  fIEntry=0;
285  } else {
286  LOG("Flux", pWARN)
287  << "No more entries in input flux neutrino ntuple, cycle "
288  << fICycle << " of " << fNCycles;
289  fEnd = true;
290  //assert(0);
291  return false;
292  }
293  }
294 
295  fNuFluxTree->GetEntry(fIEntry);
296 
297 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
298  LOG("Flux",pDEBUG)
299  << "got " << fNNeutrinos << " new fIEntry " << fIEntry
300  << " pot# " << fCurDk2Nu->potnum;
301 #endif
302 
303  fIUse = 1;
304 
305  // here we might want to do flavor oscillations or simple mappings
306  fCurNuChoice->pdgNu = fCurDk2Nu->decay.ntype;
307  fCurNuChoice->impWgt = fCurDk2Nu->decay.nimpwt;
308  }
309 
310  // update the # POTs & number of neutrinos
311  // Do this HERE (before rejecting flavors that users might be weeding out)
312  // in order to keep the POT accounting correct. This allows one to get
313  // the right normalization for generating only events from the intrinsic
314  // nu_e entries.
315  fAccumPOTs += fEffPOTsPerNu / fMaxWeight;
316  fNNeutrinos++;
317 
318  // Check neutrino pdg against declared list of neutrino species declared
319  // by the current instance of the neutrino flux driver.
320  // No undeclared neutrino species will be accepted at this point as GENIE
321  // has already been configured to handle the specified list.
322  // Make sure that the appropriate list of flux neutrino species was set at
323  // initialization via GDk2NuFlux::SetFluxParticles(const PDGCodeList &)
324 
325  if ( ! fPdgCList->ExistsInPDGCodeList(fCurNuChoice->pdgNu) ) {
326  /// user might modify list via SetFluxParticles() in order to reject certain
327  /// flavors, even if they're found in the file. So don't make a big fuss.
328  /// Spit out a single message and then stop reporting that flavor as problematic.
329  int badpdg = fCurNuChoice->pdgNu;
330  if ( ! fPdgCListRej->ExistsInPDGCodeList(badpdg) ) {
331  fPdgCListRej->push_back(badpdg);
332  LOG("Flux", pWARN)
333  << "Encountered neutrino specie (" << badpdg
334  << " that wasn't in SetFluxParticles() list, "
335  << "\nDeclared list of neutrino species: " << *fPdgCList;
336  }
337  return false;
338  }
339 
340  // Update the curr neutrino weight and energy
341 
342  // Check current neutrino energy against the maximum flux neutrino energy
343  // declared by the current instance of the neutrino flux driver.
344  // No flux neutrino exceeding that maximum energy will be accepted at this
345  // point as that maximum energy has already been used for normalizing the
346  // interaction probabilities.
347  // Make sure that the appropriate maximum flux neutrino energy was set at
348  // initialization via GDk2NuFlux::SetMaxEnergy(double Ev)
349 
350  fCurNuChoice->x4NuBeam = fFluxWindowBase;
351 
352  double Ev = 0;
353  double& wgt_xy = fCurNuChoice->xyWgt;
354  // recalculate on x-y window
356  fCurNuChoice->x4NuBeam += ( rnd->RndFlux().Rndm()*fFluxWindowDir1 +
357  rnd->RndFlux().Rndm()*fFluxWindowDir2 );
358  bsim::calcEnuWgt(fCurDk2Nu->decay,fCurNuChoice->x4NuBeam.Vect(),Ev,wgt_xy);
359 
360  if (Ev > fMaxEv) {
361  LOG("Flux", pWARN)
362  << "Flux neutrino energy exceeds declared maximum neutrino energy"
363  << "\nEv = " << Ev << "(> Ev{max} = " << fMaxEv << ")";
364  }
365 
366  // Set the current flux neutrino 4-momentum
367  // this is in *beam* coordinates
368  fgX4dkvtx = TLorentzVector( fCurDk2Nu->decay.vx,
369  fCurDk2Nu->decay.vy,
370  fCurDk2Nu->decay.vz, 0.);
371  // don't use TLorentzVector here for Mag() due to - on metric
372  TVector3 dirNu = (fCurNuChoice->x4NuBeam.Vect() - fgX4dkvtx.Vect()).Unit();
373  fCurNuChoice->p4NuBeam.SetPxPyPzE( Ev*dirNu.X(),
374  Ev*dirNu.Y(),
375  Ev*dirNu.Z(), Ev);
376 
377  // Set the current flux neutrino 4-position, direction in user coord
378  Beam2UserP4(fCurNuChoice->p4NuBeam,fCurNuChoice->p4NuUser);
379  Beam2UserPos(fCurNuChoice->x4NuBeam,fCurNuChoice->x4NuUser);
380 
381  fWeight = fCurNuChoice->impWgt * fCurNuChoice->xyWgt; // full weight
382  if ( fApplyTiltWeight ) {
383  // additional weight due to window tilt relative to flux direction
384  double tiltwgt = dirNu.Dot( FluxWindowNormal() );
385  fWeight *= TMath::Abs( tiltwgt );
386  }
387 
388  // set the time component of the Lorentz vectors
389  double dist_dk2start = GetDecayDist();
390  size_t inu = fCurDk2Nu->indxnu();
391  double t_dk = 0;
392  if ( ! fCurDk2Nu->overflow() && ! fCurDk2Nu->ancestor.empty() ) {
393  t_dk = fCurDk2Nu->ancestor[inu].startt; // units? seconds, hopefully
394  }
395  const double c_mbys = 299792458;
396  const double c_cmbys = c_mbys * 100.;
397  double tstart = t_dk + ( dist_dk2start / c_cmbys );
398  fCurNuChoice->x4NuBeam.SetT(tstart);
399  fCurNuChoice->x4NuUser.SetT(tstart);
400  if ( t_dk == 0 ) {
401  // probably wasn't set
402  static int nmsg = 5;
403  if ( nmsg > 0 ) {
404  --nmsg;
405  LOG("Flux", pNOTICE)
406  << "Setting time at flux window, \n"
407  << "noticed that t_dk from ancestor list was 0, "
408  << "this probably means that the calculated time is wrong";
409  }
410  }
411 
412  // if desired, move to user specified user coord z
413  if ( TMath::Abs(fZ0) < 1.0e30 ) this->MoveToZ0(fZ0);
414 
415 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
416  LOG("Flux", pINFO)
417  << "Generated neutrino: " << fIEntry << " " << fCurDk2Nu->potnum
418  << "\n pdg-code: " << fCurNuChoice->pdgNu
419  << "\n p4 beam: " << utils::print::P4AsShortString(&fCurNuChoice->p4NuBeam)
420  << "\n x4 beam: " << utils::print::X4AsString(&fCurNuChoice->x4NuBeam)
421  << "\n p4 user: " << utils::print::P4AsShortString(&(fCurNuChoice->p4NuUser))
422  << "\n x4 user: " << utils::print::X4AsString(&(fCurNuChoice->x4NuUser));
423 #endif
424  if ( Ev > fMaxEv ) {
425  LOG("Flux", pFATAL)
426  << "Generated neutrino had E_nu = " << Ev << " > " << fMaxEv
427  << " maximum ";
428  assert(0);
429  }
430 
431  // update sum of weights
432  fSumWeight += this->Weight();
433 
434  return true;
435 }
436 //___________________________________________________________________________
437 double GDk2NuFlux::GetDecayDist() const
438 {
439  // return distance (user units) between dk point and start position
440  // these are in beam units
441  TVector3 x3diff = fCurNuChoice->x4NuBeam.Vect() - fgX4dkvtx.Vect();
442  return x3diff.Mag() * fLengthScaleB2U;
443 }
444 //___________________________________________________________________________
445 void GDk2NuFlux::MoveToZ0(double z0usr)
446 {
447  // move ray origin to specified user z0
448  // move beam coord entry correspondingly
449 
450 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
451  LOG("Flux", pDEBUG)
452  << "MoveToZ0 (z0usr=" << z0usr << ") before:"
453  << "\n p4 user: " << utils::print::P4AsShortString(&(fCurNuChoice->p4NuUser))
454  << "\n x4 user: " << utils::print::X4AsString(&(fCurNuChoice->x4NuUser));
455 #endif
456 
457  double pzusr = fCurNuChoice->p4NuUser.Pz();
458  if ( TMath::Abs(pzusr) < 1.0e-30 ) {
459  // neutrino is moving almost entirely in x-y plane
460  LOG("Flux", pWARN)
461  << "MoveToZ0(" << z0usr << ") not possible due to pz_usr (" << pzusr << ")";
462  return;
463  }
464 
465  double scale = (z0usr - fCurNuChoice->x4NuUser.Z()) / pzusr;
466  fCurNuChoice->x4NuUser += (scale*fCurNuChoice->p4NuUser);
467  fCurNuChoice->x4NuBeam += ((fLengthScaleU2B*scale)*fCurNuChoice->p4NuBeam);
468  // this scaling works for distances, but not the time component
469  fCurNuChoice->x4NuBeam.SetT(0);
470  fCurNuChoice->x4NuUser.SetT(0);
471 
472 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
473  LOG("Flux", pDEBUG)
474  << "MoveToZ0 (" << z0usr << ") after:"
475  << "\n x4 user: " << utils::print::X4AsString(&(fCurNuChoice->x4NuUser));
476 #endif
477 
478 }
479 
480 //___________________________________________________________________________
482 {
483  // do this if flux window changes or # of files changes
484 
485  if (!fNuFluxTree) return; // not yet fully configured
486 
487  // effpots = mc_pots * (wgtfunction-area) / window-area / wgt-max-est
488  // wgtfunction-area = pi * radius-det-element^2 = pi * (100.cm)^2
489 
490  // this should match what is used in the CalcEnuWgt()
491  const double kRDET = 100.0; // set to flux per 100 cm radius
492  const double kRDET2 = kRDET * kRDET;
493  double flux_area = fFluxWindowDir1.Vect().Cross(fFluxWindowDir2.Vect()).Mag();
494  LOG("Flux",pNOTICE) << "in CalcEffPOTsPerNu, area = " << flux_area;
495 
496  if ( flux_area < 1.0e-30 ) {
497  LOG("Flux", pWARN)
498  << "CalcEffPOTsPerNu called with flux window area effectively zero";
499  flux_area = 1;
500  }
501  double area_ratio = TMath::Pi() * kRDET2 / flux_area;
502  fEffPOTsPerNu = area_ratio * ( (double)fFilePOTs / (double)fNEntries );
503 }
504 
505 //___________________________________________________________________________
506 void GDk2NuFlux::LoadDkMeta(void)
507 {
508  // load the matching bsim::dkmeta entry that goes w/ bsim::dk2nu entry
509 
510  if ( fJobToMetaIndex.empty() ) {
511  int nmeta = fNuMetaTree->GetEntries();
512  for (int imeta =0; imeta < nmeta; ++imeta ) {
513  fNuMetaTree->GetEntry(imeta);
514  int mjob = fCurDkMeta->job;
515  // there shouldn't already be an entry in the map
516  // complain if there is
517  std::map<int,int>::const_iterator mitr = fJobToMetaIndex.find(mjob);
518  if ( mitr == fJobToMetaIndex.end() ) {
519  fJobToMetaIndex[mjob] = imeta; // make an entry
520  } else {
521  LOG("Flux", pERROR) << "LoadDkMeta already had an entry for job "
522  << mjob << " at " << mitr->second
523  << " which conflicts with new entry at "
524  << imeta;
525  }
526  }
527  }
528 
529  int job = fCurDk2Nu->job;
530  if ( fCurDkMeta->job == job ) return; // already loaded
531 
532  std::map<int,int>::const_iterator lookat = fJobToMetaIndex.find(job);
533  if ( lookat != fJobToMetaIndex.end() ) {
534  int indx = lookat->second;
535  fNuMetaTree->GetEntry(indx);
536  if ( fCurDkMeta->job != job ) {
537  LOG("Flux", pFATAL) << "Failed to get right metadata " << job
538  << " => " << fCurDkMeta->job
539  << " indx " << indx;
540  assert(0);
541  }
542  } else {
543  // wasn't already indexed
544  LOG("Flux", pFATAL) << "Failed index metadata for job " << job
545  << " indx " << lookat->second;
546  assert(0);
547  }
548 
549 }
550 
551 //___________________________________________________________________________
552 double GDk2NuFlux::UsedPOTs(void) const
553 {
554 // Compute current number of flux POTs
555 
556  if (!fNuFluxTree) {
557  LOG("Flux", pWARN)
558  << "The flux driver has not been properly configured";
559  return 0;
560  }
561  return fAccumPOTs;
562 }
563 
564 //___________________________________________________________________________
565 double GDk2NuFlux::POT_curr(void) {
566  // RWH: Not sure what POT_curr is supposed to represent I'll guess for
567  // now that that it means what I mean by UsedPOTs().
568  return UsedPOTs();
569 }
570 
571 //___________________________________________________________________________
572 void GDk2NuFlux::LoadBeamSimData(const std::vector<std::string>& patterns,
573  const std::string& config )
574 {
575 // Loads in a beam simulation root file into the GDk2NuFlux driver.
576 
577  bool found_cfg = this->LoadConfig(config);
578  if ( ! found_cfg ) {
579  LOG("Flux", pFATAL)
580  << "LoadBeamSimData could not find XML config \"" << config << "\"\n";
581  exit(1);
582  }
583 
584  fNuFluxFilePatterns = patterns;
585  std::vector<int> nfiles_from_pattern;
586 
587  // create a (sorted) set of file names
588  // this also helps ensure that the same file isn't listed multiple times
589  std::set<std::string> fnames;
590 
591  LOG("Flux", pINFO) << "LoadBeamSimData was passed " << patterns.size()
592  << " patterns";
593 
594  for (size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
595  std::string pattern = patterns[ipatt];
596  nfiles_from_pattern.push_back(0);
597 
598  LOG("Flux", pNOTICE)
599  << "Loading dk2nu flux tree from ROOT file pattern ["
600  << std::setw(3) << ipatt << "] \"" << pattern << "\"";
601 
602  // !WILDCARD only works for file name ... NOT directory
603  std::string dirname = gSystem->UnixPathName(gSystem->WorkingDirectory());
604  size_t slashpos = pattern.find_last_of("/");
605  size_t fbegin;
606  if ( slashpos != std::string::npos ) {
607  dirname = pattern.substr(0,slashpos);
608  LOG("Flux", pINFO) << "Look for flux using directory " << dirname;
609  fbegin = slashpos + 1;
610  } else { fbegin = 0; }
611 
612  void* dirp = gSystem->OpenDirectory(gSystem->ExpandPathName(dirname.c_str()));
613  if ( dirp ) {
614  std::string basename =
615  pattern.substr(fbegin,pattern.size()-fbegin);
616  TRegexp re(basename.c_str(),kTRUE);
617  const char* onefile;
618  while ( ( onefile = gSystem->GetDirEntry(dirp) ) ) {
619  TString afile = onefile;
620  if ( afile=="." || afile==".." ) continue;
621  if ( basename!=afile && afile.Index(re) == kNPOS ) continue;
622  std::string fullname = dirname + "/" + afile.Data();
623  fnames.insert(fullname);
624  nfiles_from_pattern[ipatt]++;
625  }
626  gSystem->FreeDirectory(dirp);
627  } // legal directory
628  } // loop over patterns
629 
630  // std::cout << "RWH start adding files" << std::endl << std::flush;
631 
632  size_t indx = 0;
633  std::set<std::string>::const_iterator sitr = fnames.begin();
634  for ( ; sitr != fnames.end(); ++sitr, ++indx ) {
635  std::string filename = *sitr;
636  // std::cout << "RWH [" << std::setw(3) << indx << "] \""
637  // << filename << "\"" << std::endl << std::flush;
638  bool isok = true;
639  // this next test only works for local files, so we can't do that
640  // if we want to stream via xrootd
641  // ! (gSystem->AccessPathName(filename.c_str()));
642  if ( isok ) {
643  TFile* tf = TFile::Open(filename.c_str(),"READ");
644  TTree* ftree = (TTree*)tf->Get(fTreeNames[0].c_str());
645  TTree* mtree = (TTree*)tf->Get(fTreeNames[1].c_str());
646  if ( ftree && mtree ) {
647  // found both trees
648  this->AddFile(ftree,mtree,filename);
649  // std::cout << "RWH AddFile(" << ftree << ","
650  // << mtree << "," << filename << ")" << std::endl << std::flush;
651  } else {
652  // std::cout << "RWH no AddFile" << std::endl << std::flush;
653  LOG("Flux", pNOTICE) << "File " << filename << " lacked a tree: "
654  << " \"" << fTreeNames[0] << "\" " << ftree
655  << " \"" << fTreeNames[1] << "\" " << mtree
656  << " and was not added to the file list";
657  }
658  tf->Close();
659  delete tf;
660  } // loop over tree type
661  } // loop over sorted file names
662 
663  if ( ! fNuFluxTree || ! fNuMetaTree ) {
664  LOG("Flux", pERROR)
665  << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
666  LOG("Flux", pERROR)
667  << "Missing dk2nu (" << fNuFluxTree << ") or dkmeta ("
668  << fNuMetaTree << ") tree";
669  LOG("Flux", pERROR)
670  << "please check that input files are valid";
671  LOG("Flux", pERROR)
672  << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
673  }
674 
675  // this will open all files and read header!!
676  fNEntries = fNuFluxTree->GetEntries();
677 
678  if ( fNEntries == 0 ) {
679  LOG("Flux", pERROR)
680  << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
681  LOG("Flux", pERROR)
682  << "Loaded flux tree contains " << fNEntries << " entries";
683  LOG("Flux", pERROR)
684  << "Was passed the file patterns: ";
685  for (size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
686  std::string pattern = patterns[ipatt];
687  LOG("Flux", pERROR)
688  << " [" << std::setw(3) << ipatt <<"] " << pattern;
689  }
690  LOG("Flux", pERROR)
691  << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
692  } else {
693  LOG("Flux", pNOTICE)
694  << "Loaded flux tree contains " << fNEntries << " entries"
695  << " from " << fnames.size() << " unique files";
696  for (size_t ipatt = 0; ipatt < patterns.size(); ++ipatt ) {
697  std::string pattern = patterns[ipatt];
698  LOG("Flux", pINFO)
699  << " pattern: " << pattern << " yielded "
700  << nfiles_from_pattern[ipatt] << " files";
701  }
702  }
703 
704  TBranch* bdk2nu = fNuFluxTree->GetBranch("dk2nu");
705  if (bdk2nu) bdk2nu->SetAutoDelete(false);
706  TBranch* bdkmeta = fNuMetaTree->GetBranch("dkmeta");
707  if (bdkmeta) bdkmeta->SetAutoDelete(false);
708 
709  // we have a file we can work with
710  if (!fDetLocIsSet) {
711  LOG("Flux", pERROR)
712  << "LoadBeamSimData left detector location unset";
713  }
714  if (fMaxWeight<=0) {
715  LOG("Flux", pINFO)
716  << "Run ScanForMaxWeight() as part of LoadBeamSimData";
717  this->ScanForMaxWeight();
718  }
719 
720  // current ntuple cycle # (flux ntuples may be recycled)
721  fICycle = 0;
722  // pick a starting entry index [0:fNEntries-1]
723  // pretend we just used up the the previous one
725  fIUse = 9999999;
726  fIEntry = rnd->RndFlux().Integer(fNEntries) - 1;
727 
728  // don't count things we used to estimate max weight
729  fSumWeight = 0;
730  fNNeutrinos = 0;
731  fAccumPOTs = 0;
732 
733  LOG("Flux",pNOTICE) << "about to CalcEffPOTsPerNu";
734  this->CalcEffPOTsPerNu();
735 
736 }
737 
738 // #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,9,0)
739 //___________________________________________________________________________
740 void GDk2NuFlux::GetBranchInfo(std::vector<std::string>& branchNames,
741  std::vector<std::string>& branchClassNames,
742  std::vector<void**>& branchObjPointers)
743 {
744  // allow flux driver to report back current status and/or ntuple entry
745  // info for possible recording in the output file by supplying
746  // the class name, and a pointer to the object that will be filled
747  // as well as a suggested name for the branch.
748 
749  branchNames.push_back("dk2nu");
750  branchClassNames.push_back("bsim::Dk2Nu");
751  branchObjPointers.push_back((void**)&fCurDk2Nu);
752 
753  branchNames.push_back("nuchoice");
754  branchClassNames.push_back("bsim::NuChoice");
755  branchObjPointers.push_back((void**)&fCurNuChoice);
756 
757 }
758 TTree* GDk2NuFlux::GetMetaDataTree() { return fNuMetaTree; }
759 
760 // #else
761 // //___________________________________________________________________________
762 // // migrated to GFluxFileConfigI
763 // void GDk2NuFlux::SetUpstreamZ(double z0) { fZ0 = z0; }
764 // void GDk2NuFlux::SetNumOfCycles(long int ncycle) { fNCycles = TMath::Max(0L, ncycle); }
765 // void GDk2NuFlux::LoadBeamSimData(std::set<std::string> fileset, std::string config)
766 // {
767 // // have a set<> want a vector<>
768 // std::vector<std::string> filevec;
769 // std::copy(fileset.begin(),fileset.end(),std::back_inserter(filevec));
770 // LoadBeamSimData(filevec,config); // call the one that takes a vector
771 // }
772 // void GDk2NuFlux::LoadBeamSimData(std::string filename, std::string config)
773 // {
774 // // Loads a beam simulation root file into the GDk2NuFlux driver.
775 // std::vector<std::string> filevec;
776 // filevec.push_back(filename);
777 // LoadBeamSimData(filevec,config); // call the one that takes a vector
778 // }
779 // void GDk2NuFlux::SetFluxParticles(const PDGCodeList & particles)
780 // {
781 // if (!fPdgCList) {
782 // fPdgCList = new PDGCodeList;
783 // }
784 // fPdgCList->Copy(particles);
785 //
786 // LOG("Flux", pINFO)
787 // << "Declared list of neutrino species: " << *fPdgCList;
788 // }
789 // #endif
790 
791 //___________________________________________________________________________
793 {
794  if (!fDetLocIsSet) {
795  LOG("Flux", pERROR)
796  << "Specify a detector location before scanning for max weight";
797  return;
798  }
799 
800  // scan for the maximum weight
801 
802  double wgtgenmx = 0, enumx = 0;
803  TStopwatch t;
804  t.Start();
805  for (int itry=0; itry < fMaxWgtEntries; ++itry) {
806  this->GenerateNext_weighted();
807  double wgt = this->Weight();
808  if ( wgt > wgtgenmx ) wgtgenmx = wgt;
809  double enu = fCurNuChoice->p4NuBeam.Energy();
810  if ( enu > enumx ) enumx = enu;
811  }
812  t.Stop();
813  t.Print("u");
814  LOG("Flux", pNOTICE) << "Maximum flux weight for spin = "
815  << wgtgenmx << ", energy = " << enumx
816  << " (" << fMaxWgtEntries << ")";
817 
818  if (wgtgenmx > fMaxWeight ) fMaxWeight = wgtgenmx;
819  // apply a fudge factor to estimated weight
820  fMaxWeight *= fMaxWgtFudge;
821  if ( fMaxWgtFudge < 1 ) fMaxWgtFudge=1; // avoid run-away downward spiral
822 
823  fMaxWeightScan = fMaxWeight;
824  // apply possible user supplied minimum (floor)
825  if (fMinMaxWeight > fMaxWeight) {
826  LOG("Flux", pNOTICE) << "Maximum flux weight set by user = "
827  << fMinMaxWeight << " exceeds estimate, use that";
828  fMaxWeight = fMinMaxWeight;
829  }
830 
831  fMaxWeightInit = fMaxWeight; // what we started with before generation
832  fMaxWeightMax = fMaxWeight;
833  fMaxWgtExceeded = 0; // not yet exceeded
834 
835  // adjust max energy?
836  if ( enumx*fMaxEFudge > fMaxEv ) {
837  LOG("Flux", pNOTICE) << "Adjust max: was=" << fMaxEv
838  << " now " << enumx << "*" << fMaxEFudge
839  << " = " << enumx*fMaxEFudge;
840  fMaxEv = enumx * fMaxEFudge;
841  }
842 
843  LOG("Flux", pNOTICE) << "Maximum flux weight = " << fMaxWeight
844  << ", energy = " << fMaxEv;
845 
846 }
847 //___________________________________________________________________________
848 void GDk2NuFlux::SetMaxEnergy(double Ev)
849 {
850  fMaxEv = TMath::Max(0.,Ev);
851 
852  LOG("Flux", pINFO)
853  << "Declared maximum flux neutrino energy: " << fMaxEv;
854 }
855 //___________________________________________________________________________
856 void GDk2NuFlux::SetEntryReuse(long int nuse)
857 {
858 // With nuse > 1 then the same entry in the file is used "nuse" times
859 // before moving on to the next entry in the ntuple
860 
861  fNUse = TMath::Max(1L, nuse);
862 }
863 
864 //___________________________________________________________________________
865 void GDk2NuFlux::SetFluxWindow(TVector3 p0, TVector3 p1, TVector3 p2)
866  // bool inDetCoord) future extension
867 {
868  // set flux window
869  // NOTE: internally these are in "cm", but user might have set a preference
870  fDetLocIsSet = true;
871  fIsSphere = false;
872 
873  fFluxWindowPtUser[0] = p0;
874  fFluxWindowPtUser[1] = p1;
875  fFluxWindowPtUser[2] = p2;
876 
877  // convert from user to beam coord and from 3 points to base + 2 directions
878  // apply units conversion
879  TLorentzVector ptbm0, ptbm1, ptbm2;
880  User2BeamPos(TLorentzVector(fFluxWindowPtUser[0],0),ptbm0);
881  User2BeamPos(TLorentzVector(fFluxWindowPtUser[1],0),ptbm1);
882  User2BeamPos(TLorentzVector(fFluxWindowPtUser[2],0),ptbm2);
883 
884  fFluxWindowBase = ptbm0;
885  fFluxWindowDir1 = ptbm1 - ptbm0;
886  fFluxWindowDir2 = ptbm2 - ptbm0;
887 
888  fFluxWindowLen1 = fFluxWindowDir1.Mag();
889  fFluxWindowLen2 = fFluxWindowDir2.Mag();
890  fFluxWindowNormal = fFluxWindowDir1.Vect().Cross(fFluxWindowDir2.Vect()).Unit();
891 
892  double dot = fFluxWindowDir1.Dot(fFluxWindowDir2);
893  if ( TMath::Abs(dot) > 1.0e-8 )
894  LOG("Flux",pWARN) << "Dot product between window direction vectors was "
895  << dot << "; please check for orthoganality";
896 
897  LOG("Flux",pNOTICE) << "about to CalcEffPOTsPerNu";
898  this->CalcEffPOTsPerNu();
899 }
900 
901 //___________________________________________________________________________
902 void GDk2NuFlux::SetFluxSphere(TVector3 center, double radius,
903  bool inDetCoord)
904 {
905  fIsSphere = true;
906  fFluxSphereRadius = radius;
907  if ( inDetCoord ) {
908  fFluxSphereCenterUser = center;
909  TLorentzVector v4beam;
910  User2BeamPos(TLorentzVector(center,0),v4beam);
911  fFluxSphereCenterBeam = v4beam.Vect();
912  } else {
913  fFluxSphereCenterBeam = center;
914  TLorentzVector v4user;
915  Beam2UserPos(TLorentzVector(center,0),v4user);
916  fFluxSphereCenterBeam = v4user.Vect();
917  }
918 
919  // not yet supported
920  LOG("Flux", pFATAL) << "use of spherical flux window not yet fully supported";
921  assert(0);
922 
923 }
924 
925 //___________________________________________________________________________
926 void GDk2NuFlux::GetFluxSphere(TVector3& center, double& radius,
927  bool inDetCoord) const
928 {
929  center = ( inDetCoord ) ? fFluxSphereCenterUser : fFluxSphereCenterBeam;
930  radius = fFluxSphereRadius;
931 }
932 
933 //___________________________________________________________________________
934 void GDk2NuFlux::GetFluxWindow(TVector3& p0, TVector3& p1, TVector3& p2) const
935 {
936  // return flux window points
937  p0 = fFluxWindowPtUser[0];
938  p1 = fFluxWindowPtUser[1];
939  p2 = fFluxWindowPtUser[2];
940 
941 }
942 //___________________________________________________________________________
943 void GDk2NuFlux::SetBeamRotation(TRotation beamrot)
944 {
945  // rotation is really only 3-d vector, but we'll be operating on LorentzV's
946  fBeamRot = TLorentzRotation(beamrot);
947  fBeamRotInv = fBeamRot.Inverse();
948 }
949 
950 void GDk2NuFlux::SetBeamCenter(TVector3 beam0)
951 {
952  // set coord transform between detector and beam
953  // NOTE: internally these are in "cm", but user might have set a preference
954  fBeamZero = TLorentzVector(beam0,0); // no time shift
955 }
956 
957 //___________________________________________________________________________
958 TRotation GDk2NuFlux::GetBeamRotation() const
959 {
960  // rotation is really only 3-d vector, but we'll be operating on LorentzV's
961  // give people back the original TRotation ... not pretty
962  // ... it think this is right
963  TRotation rot3;
964  const TLorentzRotation& rot4 = fBeamRot;
965  TVector3 newX(rot4.XX(),rot4.XY(),rot4.XZ());
966  TVector3 newY(rot4.YX(),rot4.YY(),rot4.YZ());
967  TVector3 newZ(rot4.ZX(),rot4.ZY(),rot4.ZZ());
968  rot3.RotateAxes(newX,newY,newZ);
969  return rot3.Inverse();
970 }
971 TVector3 GDk2NuFlux::GetBeamCenter() const
972 {
973  TVector3 beam0 = fBeamZero.Vect();
974  return beam0;
975 }
976 
977 //___________________________________________________________________________
978 TVector3 GDk2NuFlux::SphereNormal() const
979 {
980  // return the normal for this current location on the sphere
981  TVector3 dir = (fCurNuChoice->x4NuBeam.Vect() - fFluxSphereCenterBeam);
982  return dir.Unit();
983 }
984 
985 //___________________________________________________________________________
986 //void GDk2NuFlux::SetCoordTransform(TVector3 beam0, TRotation beamrot)
987 //{
988 // // set coord transform between detector and beam
989 // // NOTE: internally these are in "cm", but user might have set a preference
990 //
991 // beam0 *= (1./fLengthScaleB2U);
992 // fDetectorZero = TLorentzVector(beam0,0); // no time shift
993 // fDetectorRot = TLorentzRotation(beamrot);
994 //
995 //}
996 //___________________________________________________________________________
997 //void GDk2NuFlux::GetDetectorCoord(TLorentzVector& det0, TLorentzRotation& detrot) const
998 //{
999 // // get coord transform between detector and beam
1000 // // NOTE: internally these are in "cm", but user might have set a preference
1001 //
1002 // det0 = fDetectorZero;
1003 // det0 *= fLengthScaleB2U;
1004 // detrot = fDetectorRot;
1005 //
1006 //}
1007 //___________________________________________________________________________
1008 
1009 void GDk2NuFlux::Beam2UserPos(const TLorentzVector& beamxyz,
1010  TLorentzVector& usrxyz) const
1011 {
1012  usrxyz = fLengthScaleB2U*(fBeamRot*beamxyz) + fBeamZero;
1013 }
1014 void GDk2NuFlux::Beam2UserDir(const TLorentzVector& beamdir,
1015  TLorentzVector& usrdir) const
1016 {
1017  usrdir = fLengthScaleB2U*(fBeamRot*beamdir);
1018 }
1019 void GDk2NuFlux::Beam2UserP4 (const TLorentzVector& beamp4,
1020  TLorentzVector& usrp4 ) const
1021 {
1022  usrp4 = fBeamRot*beamp4;
1023 }
1024 
1025 void GDk2NuFlux::User2BeamPos(const TLorentzVector& usrxyz,
1026  TLorentzVector& beamxyz) const
1027 {
1028  beamxyz = fLengthScaleU2B*(fBeamRotInv*(usrxyz-fBeamZero));
1029 }
1030 void GDk2NuFlux::User2BeamDir(const TLorentzVector& usrdir,
1031  TLorentzVector& beamdir) const
1032 {
1033  beamdir = fLengthScaleU2B*(fBeamRotInv*usrdir);
1034 }
1035 void GDk2NuFlux::User2BeamP4 (const TLorentzVector& usrp4,
1036  TLorentzVector& beamp4) const
1037 {
1038  beamp4 = fBeamRotInv*usrp4;
1039 }
1040 
1041 //___________________________________________________________________________
1042 void GDk2NuFlux::PrintCurrent(void)
1043 {
1044  LOG("Flux", pNOTICE) << "CurrentEntry:\n"
1045  << fCurDk2Nu->AsString() << "\n"
1046  << fCurNuChoice->AsString();
1047 }
1048 //___________________________________________________________________________
1049 void GDk2NuFlux::Clear(Option_t * opt)
1050 {
1051  // Clear the driver state
1052  //
1053  LOG("Flux", pWARN) << "GSimpleNtpFlux::Clear(" << opt << ") called";
1054  // do it in all cases, but EVGDriver/GMCJDriver will pass "CycleHistory"
1055 
1056  fICycle = 0;
1057 
1058  fSumWeight = 0;
1059  fNNeutrinos = 0;
1060  fAccumPOTs = 0;
1061 }
1062 //___________________________________________________________________________
1063 void GDk2NuFlux::GenerateWeighted(bool gen_weighted)
1064 {
1065  // Set whether to generate weighted rays
1066  //
1067  fGenWeighted = gen_weighted;
1068 }
1069 //___________________________________________________________________________
1070 void GDk2NuFlux::Initialize(void)
1071 {
1072  LOG("Flux", pNOTICE) << "Initializing GDk2NuFlux driver";
1073 
1074  fMaxEv = 0;
1075  fEnd = false;
1076 
1077 // #if __GENIE_RELEASE_CODE__ <= GRELCODE(2,9,0)
1078 // fPdgCList = new PDGCodeList;
1079 // fPdgCListRej = new PDGCodeList;
1080 // #endif
1081 
1082  fTreeNames[0] = "dk2nuTree";
1083  fTreeNames[1] = "dkmetaTree";
1084  fNuFluxTree = 0;
1085  fNuMetaTree = 0;
1086  fCurDk2Nu = 0;
1087  fCurDkMeta = 0;
1088  fCurNuChoice = 0;
1089  fNFiles = 0;
1090 
1091  fNEntries = 0;
1092  fIEntry = -1;
1093  fICycle = 0;
1094  fNUse = 1;
1095  fIUse = 999999;
1096 
1097  fNuTot = 0;
1098  fFilePOTs = 0;
1099 
1100  fMaxWeight = -1;
1101  fMinMaxWeight = 0;
1102  fMaxWgtFudge = 1.05;
1103  fMaxWgtEntries = 2500000;
1104  fMaxEFudge = 0;
1105 
1106  fMaxWgtExceeded = 0;
1107  fMaxWgtFailModel = 0; // default: 0=bump, 1=leave frozen, 2=abort
1108 
1109  fSumWeight = 0;
1110  fNNeutrinos = 0;
1111  fEffPOTsPerNu = 0;
1112  fAccumPOTs = 0;
1113 
1114  fGenWeighted = false;
1115  fApplyTiltWeight = true;
1116  fIsSphere = false;
1117  fDetLocIsSet = false;
1118  // by default assume user length is m
1119  SetLengthUnits(genie::utils::units::UnitFromString("m"));
1120 
1121 // #if __GENIE_RELEASE_CODE__ < GRELCODE(2,9,0)
1122 // // migrated to GFluxFileConfigI
1123 // fXMLbasename = "GNuMIFlux.xml";
1124 // fNCycles = 0;
1125 // fZ0 = -3.4e38;
1126 // #endif
1127 
1128  this->SetDefaults();
1129  this->ResetCurrent();
1130 }
1131 //___________________________________________________________________________
1132 void GDk2NuFlux::SetDefaults(void)
1133 {
1134 // - Set default neutrino species list (nue, nuebar, numu, numubar) and
1135 // maximum energy (120 GeV).
1136 // These defaults can be overwritten by user calls (at the driver init) to
1137 // GDk2NuFlux::SetMaxEnergy(double Ev) and
1138 // GDk2NuFlux::SetFluxParticles(const PDGCodeList & particles)
1139 // - Set the default file normalization to 1E+21 POT
1140 // - Set the default flux neutrino start z position at -5m (z=0 is the
1141 // detector centre).
1142 // - Set number of cycles to 1
1143 
1144  LOG("Flux", pNOTICE) << "Setting default GDk2NuFlux driver options";
1145 
1146  PDGCodeList particles;
1147  particles.push_back(kPdgNuMu);
1148  particles.push_back(kPdgAntiNuMu);
1149  particles.push_back(kPdgNuE);
1150  particles.push_back(kPdgAntiNuE);
1151 
1152  this->SetFluxParticles (particles);
1153  this->SetMaxEnergy (120./*GeV*/); // was 200, but that would be wasteful
1154  this->SetUpstreamZ (-3.4e38); // way upstream ==> use flux window
1155  this->SetNumOfCycles (0);
1156  this->SetEntryReuse (1);
1157 
1158  std::string xmlfile = "GNuMIFlux.xml";
1159  const char* altxml = gSystem->Getenv("GDK2NUFLUXXML");
1160  if ( altxml ) xmlfile = altxml;
1161 // #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,9,0)
1162  this->SetXMLFileBase(xmlfile);
1163 // #else
1164 // this->SetXMLFile(xmlfile);
1165 // #endif
1166 
1167 }
1168 //___________________________________________________________________________
1169 void GDk2NuFlux::ResetCurrent(void)
1170 {
1171 // reset running values of neutrino pdg-code, 4-position & 4-momentum
1172 // and the input ntuple leaves
1173 
1174  if ( fCurDk2Nu ) fCurDk2Nu->clear();
1175  if ( fCurNuChoice ) fCurNuChoice->clear();
1176 }
1177 //___________________________________________________________________________
1178 void GDk2NuFlux::CleanUp(void)
1179 {
1180  LOG("Flux", pNOTICE) << "Cleaning up...";
1181 
1182  if ( fPdgCList ) delete fPdgCList;
1183  if ( fPdgCListRej ) delete fPdgCListRej;
1184  if ( fCurNuChoice ) delete fCurNuChoice;
1185 
1186  LOG("Flux", pNOTICE)
1187  << " flux file cycles: " << fICycle << " of " << fNCycles
1188  << ", entry " << fIEntry << " use: " << fIUse << " of " << fNUse;
1189 }
1190 
1191 //___________________________________________________________________________
1192 void GDk2NuFlux::AddFile(TTree* ftree, TTree* mtree, std::string fname)
1193 {
1194  // Add a file to the chain
1195 
1196  ULong64_t nentries = ftree->GetEntries();
1197 
1198  // generally files will only have one meta data entry, but if they've
1199  // been combined (i.e. "hadd") there might be more than one.
1200  int nmeta = mtree->GetEntries();
1201  double potsum = 0;
1203  mtree->SetBranchAddress("dkmeta",&dkmeta);
1204  for (int imeta = 0; imeta < nmeta; ++imeta ) {
1205  mtree->GetEntry(imeta);
1206  double potentry = dkmeta->pots;
1207  potsum += potentry;
1208  }
1209  delete dkmeta;
1210 
1211  // don't need these anymore
1212  delete ftree;
1213  delete mtree;
1214 
1215  // make sure the chains are defined and a branch object attached
1216  if ( ! fNuFluxTree ) {
1217  fNuFluxTree = new TChain(fTreeNames[0].c_str());
1218  fNuMetaTree = new TChain(fTreeNames[1].c_str());
1219  fCurDk2Nu = new bsim::Dk2Nu;
1220  fCurDkMeta = new bsim::DkMeta;
1221  fCurNuChoice = new bsim::NuChoice;
1222  fNuFluxTree->SetBranchAddress("dk2nu",&fCurDk2Nu);
1223  fNuMetaTree->SetBranchAddress("dkmeta",&fCurDkMeta);
1224  }
1225 
1226  // add the file to the chains
1227  int stat0 = fNuFluxTree->AddFile(fname.c_str());
1228  int stat1 = fNuMetaTree->AddFile(fname.c_str());
1229 
1230  LOG("Flux",pINFO)
1231  << "flux->AddFile() of " << nentries
1232  << " " << ((mtree)?"[+meta]":"[no-meta]")
1233  << " [status=" << stat0 << "," << stat1 << "]"
1234  << nentries << " (" << nmeta << ")"
1235  << " entries in file: " << fname;
1236 
1237  if ( stat0 != 1 || stat1 != 1 ) {
1238  SLOG("GDk2NuFlux", pFATAL) << "Add: \"" << fname << "\" failed";
1239  }
1240 
1241  fNuTot += nentries;
1242  fFilePOTs += potsum;
1243  fNFiles++;
1244 
1245 }
1246 
1247 //___________________________________________________________________________
1248 void GDk2NuFlux::SetLengthUnits(double user_units)
1249 {
1250  // Set the scale factor for lengths going from beam (cm) to user coordinates
1251 
1252  // GDk2NuFlux uses "cm" as the length unit consistently internally (this is
1253  // the length units used by both the g3 and g4 ntuples). User interactions
1254  // setting the beam-to-detector coordinate transform, flux window, and the
1255  // returned position might need to be in other units. Use:
1256  // double scale = genie::utils::units::UnitFromString("cm");
1257  // ( #include "Utils/UnitUtils.h for declaration )
1258  // to get the correct scale factor to pass in.
1259 
1260  double rescale = fLengthUnits / user_units;
1261  fLengthUnits = user_units;
1262  double cm = genie::utils::units::UnitFromString("cm");
1263  fLengthScaleB2U = cm / user_units;
1264  fLengthScaleU2B = user_units / cm;
1265 
1266  // in case we're changing units without resetting transform/window
1267  // not recommended, but should work
1268  if (fCurNuChoice) fCurNuChoice->x4NuUser *= rescale;
1269  fBeamZero *= rescale;
1270  fFluxWindowPtUser[0] *= rescale;
1271  fFluxWindowPtUser[1] *= rescale;
1272  fFluxWindowPtUser[2] *= rescale;
1273 
1274  // case GDk2NuFlux::kmeter: fLengthScaleB2U = 0.01 ; break;
1275  // case GDk2NuFlux::kcm: fLengthScaleB2U = 1. ; break;
1276  // case GDk2NuFlux::kmm: fLengthScaleB2U = 10. ; break;
1277  // case GDk2NuFlux::kfm: fLengthScaleB2U = 1.e13 ; break; // 10e-2m -> 10e-15m
1278 
1279 }
1280 
1281 //___________________________________________________________________________
1282 double GDk2NuFlux::LengthUnits(void) const
1283 {
1284  // Return the scale factor for lengths the user is getting
1285  double cm = genie::utils::units::UnitFromString("cm");
1286  return fLengthScaleB2U * cm ;
1287 }
1288 
1289 //___________________________________________________________________________
1290 
1292 {
1293  genie::flux::GDk2NuFluxXMLHelper helper(this);
1294  return helper.LoadConfig(cfg);
1295 }
1296 
1297 //___________________________________________________________________________
1298 
1300 {
1301 
1302  std::ostringstream s;
1303  PDGCodeList::const_iterator itr = fPdgCList->begin();
1304  for ( ; itr != fPdgCList->end(); ++itr) s << (*itr) << " ";
1305  s << "[rejected: ";
1306  itr = fPdgCListRej->begin();
1307  for ( ; itr != fPdgCListRej->end(); ++itr) s << (*itr) << " ";
1308  s << " ] ";
1309 
1310  std::ostringstream fpattout;
1311  for (size_t i = 0; i < fNuFluxFilePatterns.size(); ++i)
1312  fpattout << "\n [" << std::setw(3) << i << "] " << fNuFluxFilePatterns[i];
1313 
1314  std::ostringstream flistout;
1315  std::vector<std::string> flist = GetFileList();
1316  for (size_t i = 0; i < flist.size(); ++i)
1317  flistout << "\n [" << std::setw(3) << i << "] " << flist[i];
1318 
1319  TLorentzVector usr0(0,0,0,0);
1320  TLorentzVector usr0asbeam;
1321  User2BeamPos(usr0,usr0asbeam);
1322 
1323  const int w=10, p=6;
1324  std::ostringstream beamrot_str, beamrotinv_str;
1325  beamrot_str
1326  << "fBeamRot: " << std::setprecision(p) << "\n"
1327  << " [ "
1328  << std::setw(w) << fBeamRot.XX() << " "
1329  << std::setw(w) << fBeamRot.XY() << " "
1330  << std::setw(w) << fBeamRot.XZ() << " ]\n"
1331  << " [ "
1332  << std::setw(w) << fBeamRot.YX() << " "
1333  << std::setw(w) << fBeamRot.YY() << " "
1334  << std::setw(w) << fBeamRot.YZ() << " ]\n"
1335  << " [ "
1336  << std::setw(w) << fBeamRot.ZX() << " "
1337  << std::setw(w) << fBeamRot.ZY() << " "
1338  << std::setw(w) << fBeamRot.ZZ() << " ]";
1339  beamrotinv_str
1340  << "fBeamRotInv: " << std::setprecision(p) << "\n"
1341  << " [ "
1342  << std::setw(w) << fBeamRotInv.XX() << " "
1343  << std::setw(w) << fBeamRotInv.XY() << " "
1344  << std::setw(w) << fBeamRotInv.XZ() << " ]\n"
1345  << " [ "
1346  << std::setw(w) << fBeamRotInv.YX() << " "
1347  << std::setw(w) << fBeamRotInv.YY() << " "
1348  << std::setw(w) << fBeamRotInv.YZ() << " ]\n"
1349  << " [ "
1350  << std::setw(w) << fBeamRotInv.ZX() << " "
1351  << std::setw(w) << fBeamRotInv.ZY() << " "
1352  << std::setw(w) << fBeamRotInv.ZZ() << " ]";
1353 
1354  std::ostringstream config_str;
1355  config_str
1356  << "GDk2NuFlux Config:"
1357  << "\n Enu_max " << fMaxEv
1358  << "\n pdg-codes: " << s.str() << "\n "
1359  << fNEntries << " entries"
1360  << " (FilePOTs " << fFilePOTs << ") "
1361  << "in " << fNFiles << " files: "
1362  << flistout.str()
1363  << "\n from file patterns:"
1364  << fpattout.str()
1365  << "\n wgt max=" << fMaxWeight
1366  << " min=" << fMinMaxWeight << " scan=" << fMaxWeightScan
1367  << " ( fudge=" << fMaxWgtFudge << " using scan of "
1368  << fMaxWgtEntries << " entries )"
1369  << "\n exceeded init=" << fMaxWeightInit << ", N=" << fMaxWgtExceeded
1370  << " times (fail model " << fMaxWgtFailModel << ") max = " << fMaxWeightMax
1371  << "\n Z0 pushback " << fZ0
1372  << "\n used entry " << fIEntry << " " << fIUse << "/" << fNUse
1373  << " times, in " << fICycle << "/" << fNCycles << " cycles"
1374  << "\n SumWeight " << fSumWeight << " for " << fNNeutrinos << " neutrino entries"
1375  << "\n EffPOTsPerNu " << fEffPOTsPerNu << " AccumPOTs " << fAccumPOTs
1376  << "\n GenWeighted: \"" << (fGenWeighted?"true":"false") << "\", "
1377  << "ApplyTiltWeight: \"" << (fApplyTiltWeight?"true":"false") << "\", "
1378  << "Detector location set: \"" << (fDetLocIsSet?"true":"false") << "\", "
1379  << "\n LengthUnits " << fLengthUnits << ", scale b2u " << fLengthScaleB2U
1380  << ", scale u2b " << fLengthScaleU2B;
1381 
1382  if ( IsFluxSphere() ) {
1383  config_str
1384  << "\n Flux Sphere radius: " << fFluxSphereRadius
1385  << "\n User Center: " << utils::print::Vec3AsString(&fFluxSphereCenterUser)
1386  << "\n Beam Center: " << utils::print::Vec3AsString(&fFluxSphereCenterBeam);
1387  } else {
1388  config_str
1389  << "\n User Flux Window: "
1390  << "\n " << utils::print::Vec3AsString(&(fFluxWindowPtUser[0]))
1391  << "\n " << utils::print::Vec3AsString(&(fFluxWindowPtUser[1]))
1392  << "\n " << utils::print::Vec3AsString(&(fFluxWindowPtUser[2]))
1393  << "\n Flux Window (cm, beam coord): "
1394  << "\n base " << utils::print::X4AsString(&fFluxWindowBase)
1395  << "\n dir1 " << utils::print::X4AsString(&fFluxWindowDir1) << " len " << fFluxWindowLen1
1396  << "\n dir2 " << utils::print::X4AsString(&fFluxWindowDir2) << " len " << fFluxWindowLen2;
1397  }
1398  config_str
1399  << "\n User Beam Origin: "
1400  << "\n base " << utils::print::X4AsString(&fBeamZero)
1401  << "\n " << beamrot_str.str() << " "
1402  << "\n Detector Origin (beam coord): "
1403  << "\n base " << utils::print::X4AsString(&usr0asbeam)
1404  << "\n " << beamrotinv_str.str();
1405 
1406  LOG("Flux", pNOTICE) << config_str.str();
1407 }
1408 
1409 //___________________________________________________________________________
1410 std::vector<std::string> GDk2NuFlux::GetFileList()
1411 {
1412  std::vector<std::string> flist;
1413  if ( ! fNuFluxTree ) return flist;
1414  TObjArray *fileElements=fNuFluxTree->GetListOfFiles();
1415  if ( ! fileElements ) return flist;
1416  TIter next(fileElements);
1417  TChainElement *chEl=0;
1418  while (( chEl=(TChainElement*)next() )) {
1419  flist.push_back(chEl->GetTitle());
1420  }
1421  return flist;
1422 }
1423 
1424 //___________________________________________________________________________
1425 
1427 {
1428  // turn string into vector<double>
1429  // be liberal about separators, users might punctuate for clarity
1430  std::vector<std::string> strtokens = genie::utils::str::Split(str," ,;:()[]=");
1431  std::vector<double> vect;
1432  size_t ntok = strtokens.size();
1433 
1434  if ( fVerbose > 2 )
1435  std::cout << "GetDoubleVector \"" << str << "\"" << std::endl;
1436 
1437  for (size_t i=0; i < ntok; ++i) {
1438  std::string trimmed = utils::str::TrimSpaces(strtokens[i]);
1439  if ( " " == trimmed || "" == trimmed ) continue; // skip empty strings
1440  double val = strtod(trimmed.c_str(), (char**)NULL);
1441  if ( fVerbose > 2 )
1442  std::cout << "(" << vect.size() << ") = " << val << std::endl;
1443  vect.push_back(val);
1444  }
1445 
1446  return vect;
1447 }
1448 
1449 std::vector<long int> GDk2NuFluxXMLHelper::GetIntVector(std::string str)
1450 {
1451  // turn string into vector<long int>
1452  // be liberal about separators, users might punctuate for clarity
1453  std::vector<std::string> strtokens = genie::utils::str::Split(str," ,;:()[]=");
1454  std::vector<long int> vect;
1455  size_t ntok = strtokens.size();
1456 
1457  if ( fVerbose > 2 )
1458  std::cout << "GetIntVector \"" << str << "\"" << std::endl;
1459 
1460  for (size_t i=0; i < ntok; ++i) {
1461  std::string trimmed = utils::str::TrimSpaces(strtokens[i]);
1462  if ( " " == trimmed || "" == trimmed ) continue; // skip empty strings
1463  long int val = strtol(trimmed.c_str(),(char**)NULL,10);
1464  if ( fVerbose > 2 )
1465  std::cout << "(" << vect.size() << ") = " << val << std::endl;
1466  vect.push_back(val);
1467  }
1468  return vect;
1469 }
1470 
1472 {
1473 // #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,9,0)
1474  std::string fname = utils::xml::GetXMLFilePath(fGDk2NuFlux->GetXMLFileBase());
1475 // #else
1476 // string fname = utils::xml::GetXMLFilePath(fGDk2NuFlux->GetXMLFile());
1477 // #endif
1478 
1479  bool is_accessible = ! (gSystem->AccessPathName(fname.c_str()));
1480  if (!is_accessible) {
1481  SLOG("GDk2NuFlux", pERROR)
1482  << "The XML doc doesn't exist! (filename: " << fname << ")";
1483  return false;
1484  }
1485 
1486  xmlDocPtr xml_doc = xmlParseFile( fname.c_str() );
1487  if ( xml_doc == NULL) {
1488  SLOG("GDk2NuFlux", pERROR)
1489  << "The XML doc can't be parsed! (filename: " << fname << ")";
1490  return false;
1491  }
1492 
1493  xmlNodePtr xml_root = xmlDocGetRootElement( xml_doc );
1494  if ( xml_root == NULL ) {
1495  SLOG("GDk2NuFlux", pERROR)
1496  << "The XML doc is empty! (filename: " << fname << ")";
1497  return false;
1498  }
1499 
1500  std::string rootele = "gdk2nu_config";
1501  std::string rooteleAlt = "gnumi_config"; // read older GNuMIFlux files
1502  if ( xmlStrcmp(xml_root->name, (const xmlChar*)rootele.c_str() ) &&
1503  xmlStrcmp(xml_root->name, (const xmlChar*)rooteleAlt.c_str() ) ) {
1504  SLOG("GDk2NuFlux", pERROR)
1505  << "The XML doc has invalid root element! (filename: " << fname << ")"
1506  << " expected \"" << rootele << "\", saw \"" << xml_root->name << "\"";
1507  return false;
1508  }
1509 
1510  SLOG("GDk2NuFlux", pINFO) << "Attempt to load config \"" << cfg
1511  << "\" from file: " << fname;
1512 
1513  bool found = this->LoadParamSet(xml_doc,cfg);
1514 
1515  xmlFree(xml_doc);
1516  return found;
1517 
1518 }
1519 
1520 bool GDk2NuFluxXMLHelper::LoadParamSet(xmlDocPtr& xml_doc, std::string cfg)
1521 {
1522 
1523  xmlNodePtr xml_root = xmlDocGetRootElement( xml_doc );
1524 
1525  // loop over all xml tree nodes that are children of the root node
1526  // read the entries looking for "param_set" of the right name
1527 
1528  // loop looking for particular config
1529  bool found = false;
1530  xmlNodePtr xml_pset = xml_root->xmlChildrenNode;
1531  for ( ; xml_pset != NULL ; xml_pset = xml_pset->next ) {
1532  if ( ! xmlStrEqual(xml_pset->name, (const xmlChar*)"param_set") ) continue;
1533  // every time there is a 'param_set' tag
1534  std::string param_set_name =
1536 
1537  if ( param_set_name != cfg ) continue;
1538 
1539  SLOG("GDk2NuFlux", pINFO) << "Found config \"" << cfg;
1540 
1541  this->ParseParamSet(xml_doc,xml_pset);
1542  found = true;
1543 
1544  } // loop over elements of root
1545  xmlFree(xml_pset);
1546 
1547  return found;
1548 }
1549 
1550 void GDk2NuFluxXMLHelper::ParseParamSet(xmlDocPtr& xml_doc, xmlNodePtr& xml_pset)
1551 {
1552  xmlNodePtr xml_child = xml_pset->xmlChildrenNode;
1553  for ( ; xml_child != NULL ; xml_child = xml_child->next ) {
1554  // handle basic gnumi_config/param_set
1555  // bad cast away const on next line, but function sig requires it
1556  std::string pname =
1557  utils::xml::TrimSpaces(const_cast<xmlChar*>(xml_child->name));
1558  if ( pname == "text" || pname == "comment" ) continue;
1559  std::string pval =
1561  xmlNodeListGetString(xml_doc, xml_child->xmlChildrenNode, 1));
1562 
1563  if ( fVerbose > 1 )
1564  SLOG("GDk2NuFlux", pINFO)
1565  << " pname \"" << pname << "\", string value \"" << pval << "\"";
1566 
1567  if ( pname == "verbose" ) {
1568  fVerbose = atoi(pval.c_str());
1569 
1570  } else if ( pname == "using_param_set" ) {
1571  SLOG("GDk2NuFlux", pWARN) << "start using_param_set: \"" << pval << "\"";
1572  bool found = this->LoadParamSet(xml_doc,pval); // recurse
1573  if ( ! found ) {
1574  SLOG("GDk2NuFlux", pFATAL) << "using_param_set: \"" << pval << "\" NOT FOUND";
1575  assert(found);
1576  }
1577  SLOG("GDk2NuFlux", pWARN) << "done using_param_set: \"" << pval << "\"";
1578  } else if ( pname == "units" ) {
1580  fGDk2NuFlux->SetLengthUnits(scale);
1581  SLOG("GDk2NuFlux", pINFO) << "set user units to \"" << pval << "\"";
1582 
1583  } else if ( pname == "beamdir" ) {
1584  ParseBeamDir(xml_doc,xml_child);
1585  fGDk2NuFlux->SetBeamRotation(fBeamRotXML);
1586 
1587  } else if ( pname == "beampos" ) {
1588  ParseBeamPos(pval);
1589  fGDk2NuFlux->SetBeamCenter(fBeamPosXML);
1590 
1591  } else if ( pname == "window" ) {
1592  ParseWindowSeries(xml_doc,xml_child);
1593  // RWH !!!! MEMORY LEAK!!!!
1594  //std::cout << " flux window " << std::endl
1595  // << " [0] " << utils::print::X4AsString(new TLorentzVector(fFluxWindowPt[0],0)) << std::endl
1596  // << " [1] " << utils::print::X4AsString(new TLorentzVector(fFluxWindowPt[1],0)) << std::endl
1597  // << " [2] " << utils::print::X4AsString(new TLorentzVector(fFluxWindowPt[2],0)) << std::endl;
1598 
1599  fGDk2NuFlux->SetFluxWindow(fFluxWindowPtXML[0],
1600  fFluxWindowPtXML[1],
1601  fFluxWindowPtXML[2]);
1602 
1603  } else if ( pname == "enumax" ) {
1604  ParseEnuMax(pval);
1605 
1606  } else if ( pname == "minmaxwgt" ) {
1607  double minmaxwgt = 0;
1608  std::vector<double> v = GetDoubleVector(pval);
1609  if ( v.size() > 0 ) minmaxwgt = v[0];
1610  fGDk2NuFlux->SetMinMaxWeight(minmaxwgt);
1611  SLOG("GDk2NuFlux", pINFO) << "set minmaxwgt = " << minmaxwgt;
1612 
1613  } else if ( pname == "maxwgtfail" ) {
1614  ParseMaxWgtFail(pval);
1615 
1616  } else if ( pname == "upstreamz" ) {
1617  double z0usr = -3.4e38;
1618  std::vector<double> v = GetDoubleVector(pval);
1619  if ( v.size() > 0 ) z0usr = v[0];
1620  fGDk2NuFlux->SetUpstreamZ(z0usr);
1621  SLOG("GDk2NuFlux", pINFO) << "set upstreamz = " << z0usr;
1622 
1623  } else if ( pname == "reuse" ) {
1624  long int nreuse = 1;
1625  std::vector<long int> v = GetIntVector(pval);
1626  if ( v.size() > 0 ) nreuse = v[0];
1627  fGDk2NuFlux->SetEntryReuse(nreuse);
1628  SLOG("GDk2NuFlux", pINFO) << "set entry reuse = " << nreuse;
1629 
1630  } else {
1631  SLOG("GDk2NuFlux", pWARN)
1632  << " NOT HANDLED: pname \"" << pname
1633  << "\", string value \"" << pval << "\"";
1634 
1635  }
1636 
1637  } // loop over param_set contents
1638  xmlFree(xml_child);
1639 }
1640 
1641 void GDk2NuFluxXMLHelper::ParseBeamDir(xmlDocPtr& xml_doc, xmlNodePtr& xml_beamdir)
1642 {
1643  fBeamRotXML.SetToIdentity(); // start fresh
1644 
1645  std::string dirtype =
1647  utils::xml::GetAttribute(xml_beamdir,"type"));
1648 
1649  std::string pval =
1651  xmlNodeListGetString(xml_doc, xml_beamdir->xmlChildrenNode, 1));
1652 
1653  if ( dirtype == "series" ) {
1654  // series of rotations around an axis
1655  ParseRotSeries(xml_doc,xml_beamdir);
1656 
1657  } else if ( dirtype == "thetaphi3") {
1658  // G3 style triplet of (theta,phi) pairs
1659  std::vector<double> thetaphi3 = GetDoubleVector(pval);
1660  std::string units =
1661  utils::str::TrimSpaces(utils::xml::GetAttribute(xml_beamdir,"units"));
1662  if ( thetaphi3.size() == 6 ) {
1663  TRotation fTempRot;
1664  TVector3 newX = AnglesToAxis(thetaphi3[0],thetaphi3[1],units);
1665  TVector3 newY = AnglesToAxis(thetaphi3[2],thetaphi3[3],units);
1666  TVector3 newZ = AnglesToAxis(thetaphi3[4],thetaphi3[5],units);
1667  fTempRot.RotateAxes(newX,newY,newZ);
1668  fBeamRotXML = fTempRot; //.Inverse();
1669  } else {
1670  SLOG("GDk2NuFlux", pWARN)
1671  << " type=\"" << dirtype << "\" within <beamdir> needs 6 values";
1672  }
1673 
1674  } else if ( dirtype == "newxyz" ) {
1675  // G4 style new axis values
1676  std::vector<double> newdir = GetDoubleVector(pval);
1677  if ( newdir.size() == 9 ) {
1678  TRotation fTempRot;
1679  TVector3 newX = TVector3(newdir[0],newdir[1],newdir[2]).Unit();
1680  TVector3 newY = TVector3(newdir[3],newdir[4],newdir[5]).Unit();
1681  TVector3 newZ = TVector3(newdir[6],newdir[7],newdir[8]).Unit();
1682  fTempRot.RotateAxes(newX,newY,newZ);
1683  fBeamRotXML = fTempRot.Inverse(); // weirdly necessary: frame vs. obj rot
1684  } else {
1685  SLOG("GDk2NuFlux", pWARN)
1686  << " type=\"" << dirtype << "\" within <beamdir> needs 9 values";
1687  }
1688 
1689  } else {
1690  // yet something else ... what? 3 choices weren't sufficient?
1691  SLOG("GDk2NuFlux", pWARN)
1692  << " UNHANDLED type=\"" << dirtype << "\" within <beamdir>";
1693  }
1694 
1695  if ( fVerbose > 1 ) {
1696  int w=10, p=6;
1697  std::cout << " fBeamRotXML: " << std::setprecision(p) << std::endl;
1698  std::cout << " [ "
1699  << std::setw(w) << fBeamRotXML.XX() << " "
1700  << std::setw(w) << fBeamRotXML.XY() << " "
1701  << std::setw(w) << fBeamRotXML.XZ() << std::endl
1702  << " "
1703  << std::setw(w) << fBeamRotXML.YX() << " "
1704  << std::setw(w) << fBeamRotXML.YY() << " "
1705  << std::setw(w) << fBeamRotXML.YZ() << std::endl
1706  << " "
1707  << std::setw(w) << fBeamRotXML.ZX() << " "
1708  << std::setw(w) << fBeamRotXML.ZY() << " "
1709  << std::setw(w) << fBeamRotXML.ZZ() << " ] " << std::endl;
1710  std::cout << std::endl;
1711  }
1712 
1713 }
1714 
1716 {
1717  std::vector<double> xyz = GetDoubleVector(str);
1718  if ( xyz.size() == 3 ) {
1719  fBeamPosXML = TVector3(xyz[0],xyz[1],xyz[2]);
1720  } else if ( xyz.size() == 6 ) {
1721  // should check for '=' between triplets but we won't be so pedantic
1722  // ( userx, usery, userz ) = ( beamx, beamy, beamz )
1723  TVector3 userpos(xyz[0],xyz[1],xyz[2]);
1724  TVector3 beampos(xyz[3],xyz[4],xyz[5]);
1725  fBeamPosXML = userpos - fBeamRotXML*beampos;
1726  } else {
1727  SLOG("GDk2NuFlux", pWARN)
1728  << "Unable to parse " << xyz.size() << " values in <beampos>";
1729  return;
1730  }
1731  if ( fVerbose > 1 ) {
1732  int w=16, p=10;
1733  std::cout << " fBeamPosXML: [ " << std::setprecision(p)
1734  << std::setw(w) << fBeamPosXML.X() << " , "
1735  << std::setw(w) << fBeamPosXML.Y() << " , "
1736  << std::setw(w) << fBeamPosXML.Z() << " ] "
1737  << std::endl;
1738  }
1739 }
1740 
1742 {
1743  std::vector<double> v = GetDoubleVector(str);
1744  size_t n = v.size();
1745  if ( n > 0 ) {
1746  fGDk2NuFlux->SetMaxEnergy(v[0]);
1747  if ( fVerbose > 1 )
1748  std::cout << "ParseEnuMax SetMaxEnergy(" << v[0] << ") " << std::endl;
1749  }
1750  if ( n > 1 ) {
1751  fGDk2NuFlux->SetMaxEFudge(v[1]);
1752  if ( fVerbose > 1 )
1753  std::cout << "ParseEnuMax SetMaxEFudge(" << v[1] << ")" << std::endl;
1754  }
1755  if ( n > 2 ) {
1756  if ( n == 3 ) {
1757  fGDk2NuFlux->SetMaxWgtScan(v[2]);
1758  if ( fVerbose > 1 )
1759  std::cout << "ParseEnuMax SetMaxWgtScan(" << v[2] << ")" << std::endl;
1760  } else {
1761  long int nentries = (long int)v[3];
1762  fGDk2NuFlux->SetMaxWgtScan(v[2],nentries);
1763  if ( fVerbose > 1 )
1764  std::cout << "ParseEnuMax SetMaxWgtScan(" << v[2] << "," << nentries << ")" << std::endl;
1765  }
1766  }
1767 }
1768 
1770 {
1771  // what to do if we exceed estimate maximum weight
1772  // 0 = "bump"
1773  // 1 = "frozen" treat as frozen: accept and move on
1774  // 2 = "abort"
1775  int model = 0; // default
1776  if ( str == "0" || str == "bump" ) model = 0;
1777  else if ( str == "1" || str == "frozen" ) model = 1;
1778  else if ( str == "2" || str == "abort" ) model = 2;
1779  else {
1780  std::cout << "ParseMaxWeightFail \"" << str << "\" unrecognized "
1781  << std::endl;
1782  }
1783  if ( fVerbose > 1 ) {
1784  std::cout << "ParseMaxWeightFail \"" << str << "\" treat as "
1785  << model << std::endl;
1786  }
1787  fGDk2NuFlux->SetMaxWeightFailModel(model);
1788 }
1789 
1790 void GDk2NuFluxXMLHelper::ParseRotSeries(xmlDocPtr& xml_doc, xmlNodePtr& xml_pset)
1791 {
1792  TRotation fTempRot; // reset matrix
1793 
1794  xmlNodePtr xml_child = xml_pset->xmlChildrenNode;
1795  for ( ; xml_child != NULL ; xml_child = xml_child->next ) {
1796  // in a <beamdir> of type "series"
1797  // should be a sequence of <rotation> entries
1798  std::string name =
1799  utils::xml::TrimSpaces(const_cast<xmlChar*>(xml_child->name));
1800  if ( name == "text" || name == "comment" ) continue;
1801 
1802  if ( name == "rotation" ) {
1804  xmlNodeListGetString(xml_doc, xml_child->xmlChildrenNode, 1));
1805  std::string axis =
1807 
1808  std::string units =
1810 
1811  double rot = atof(val.c_str());
1812  // assume radians unless given a hint that it's degrees
1813  if ( 'd' == units[0] || 'D' == units[0] ) rot *= TMath::DegToRad();
1814 
1815  if ( fVerbose > 0 )
1816  SLOG("GDk2NuFlux", pINFO)
1817  << " rotate " << rot << " radians around " << axis << " axis";
1818 
1819  if ( axis[0] == 'x' || axis[0] == 'X' ) fTempRot.RotateX(rot);
1820  else if ( axis[0] == 'y' || axis[0] == 'Y' ) fTempRot.RotateY(rot);
1821  else if ( axis[0] == 'z' || axis[0] == 'Z' ) fTempRot.RotateZ(rot);
1822  else {
1823  SLOG("GDk2NuFlux", pINFO)
1824  << " no " << axis << " to rotate around";
1825  }
1826 
1827  } else {
1828  SLOG("GDk2NuFlux", pWARN)
1829  << " found <" << name << "> within <beamdir type=\"series\">";
1830  }
1831  }
1832  // TRotation rotates objects not frames, so we want the inverse
1833  fBeamRotXML = fTempRot.Inverse();
1834  xmlFree(xml_child);
1835 }
1836 
1837 void GDk2NuFluxXMLHelper::ParseWindowSeries(xmlDocPtr& xml_doc, xmlNodePtr& xml_pset)
1838 {
1839  int ientry = -1;
1840 
1841  xmlNodePtr xml_child = xml_pset->xmlChildrenNode;
1842  for ( ; xml_child != NULL ; xml_child = xml_child->next ) {
1843  // in a <windowr> element
1844  // should be a sequence of <point> entries
1845  std::string name =
1846  utils::xml::TrimSpaces(const_cast<xmlChar*>(xml_child->name));
1847  if ( name == "text" || name == "comment" ) continue;
1848 
1849  if ( name == "point" ) {
1850  std::string val =
1852  xmlNodeListGetString(xml_doc, xml_child->xmlChildrenNode, 1));
1853  std::string coord =
1855 
1856  std::vector<double> xyz = GetDoubleVector(val);
1857  if ( xyz.size() != 3 || coord != "det" ) {
1858  SLOG("GDk2NuFlux", pWARN)
1859  << "parsing <window> found <point> but size=" << xyz.size()
1860  << " (expect 3) and coord=\"" << coord << "\" (expect \"det\")"
1861  << " IGNORE problem";
1862  }
1863  ++ientry;
1864  if ( ientry < 3 && ientry >= 0 ) {
1865  TVector3 pt(xyz[0],xyz[1],xyz[2]);
1866  if ( fVerbose > 0 ) {
1867  int w=16, p=10;
1868  std::cout << " point[" << ientry <<"] = [ " << std::setprecision(p)
1869  << std::setw(w) << pt.X() << " , "
1870  << std::setw(w) << pt.Y() << " , "
1871  << std::setw(w) << pt.Z() << " ] "
1872  << std::endl;
1873  }
1874  fFluxWindowPtXML[ientry] = pt; // save the point
1875  } else {
1876  SLOG("GDk2NuFlux", pWARN)
1877  << " <window><point> ientry " << ientry << " out of range (0-2)";
1878  }
1879 
1880  } else {
1881  SLOG("GDk2NuFlux", pWARN)
1882  << " found <" << name << "> within <window>";
1883  }
1884  }
1885  xmlFree(xml_child);
1886 }
1887 
1888 TVector3 GDk2NuFluxXMLHelper::AnglesToAxis(double theta, double phi, std::string units)
1889 {
1890  double xyz[3];
1891  // assume radians unless given a hint that it's degrees
1892  double scale = ('d'==units[0]||'D'==units[0]) ? TMath::DegToRad() : 1.0 ;
1893 
1894  xyz[0] = TMath::Cos(scale*phi)*TMath::Sin(scale*theta);
1895  xyz[1] = TMath::Sin(scale*phi)*TMath::Sin(scale*theta);
1896  xyz[2] = TMath::Cos(scale*theta);
1897  // condition vector to eliminate most floating point errors
1898  for (int i=0; i<3; ++i) {
1899  const double eps = 1.0e-15;
1900  if (TMath::Abs(xyz[i]) < eps ) xyz[i] = 0;
1901  if (TMath::Abs(xyz[i]-1) < eps ) xyz[i] = 1;
1902  if (TMath::Abs(xyz[i]+1) < eps ) xyz[i] = -1;
1903  }
1904  return TVector3(xyz[0],xyz[1],xyz[2]);
1905 }
1906 
1907 TVector3 GDk2NuFluxXMLHelper::ParseTV3(const std::string& str)
1908 {
1909  std::vector<double> xyz = GetDoubleVector(str);
1910  if ( xyz.size() != 3 ) {
1911  return TVector3();
1912  SLOG("GDk2NuFlux", pWARN)
1913  << " ParseTV3 \"" << str << "\" had " << xyz.size() << " elements ";
1914  }
1915  return TVector3(xyz[0],xyz[1],xyz[2]);
1916 
1917 }
1918 //___________________________________________________________________________
void Beam2UserDir(const TLorentzVector &beamdir, TLorentzVector &usrdir) const
const XML_Char * name
Definition: expat.h:151
std::vector< double > GetDoubleVector(std::string str)
A class for generating concrete GFluxI derived classes based on the factory pattern. This code supplies a CPP macro which allows the classes to self-register and thus no modification of this class is needed in order to expand the list of classes it knows about.
double UsedPOTs(void) const
of protons-on-target used
Definition: GDk2NuFlux.cxx:552
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:52
bool LoadConfig(std::string cfg)
load a named configuration
const int kPdgNuE
Definition: PDGCodes.h:28
void PrintConfig()
print the current configuration
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
virtual void LoadBeamSimData(const std::vector< std::string > &filenames, const std::string &det_loc)
Definition: GDk2NuFlux.cxx:572
#define pERROR
Definition: Messenger.h:60
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
Definition: GDk2NuFlux.cxx:166
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
void GetFluxSphere(TVector3 &center, double &radius, bool inDetCoord=true) const
specification of a sphere
Definition: GDk2NuFlux.cxx:926
def AddFile(fout, deets, f)
string TrimSpaces(xmlChar *xmls)
void SetMaxEnergy(double Ev)
specify max flx nu energy
Definition: GDk2NuFlux.cxx:848
const int kPdgAntiNuE
Definition: PDGCodes.h:29
void ParseBeamDir(xmlDocPtr &, xmlNodePtr &)
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
const int kPdgNuMu
Definition: PDGCodes.h:30
void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
Definition: config.py:1
string filename
Definition: shutoffs.py:106
void Beam2UserP4(const TLorentzVector &beamp4, TLorentzVector &usrp4) const
double GetDecayDist() const
dist (user units) from dk to current pos
Definition: GDk2NuFlux.cxx:437
Timing fit.
void SetLengthUnits(double user_units)
Set units assumed by user.
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.C:184
Loaders::FluxType flux
void MoveToZ0(double z0)
move ray origin to user coord Z0
Definition: GDk2NuFlux.cxx:445
TVector3 SphereNormal() const
Definition: GDk2NuFlux.cxx:978
A list of PDG codes.
Definition: PDGCodeList.h:33
virtual void GetBranchInfo(std::vector< std::string > &branchNames, std::vector< std::string > &branchClassNames, std::vector< void ** > &branchObjPointers)
Definition: GDk2NuFlux.cxx:740
Double_t pots
protons-on-target
Definition: dkmeta.h:87
void SetEntryReuse(long int nuse=1)
of times to use entry before moving to next
Definition: GDk2NuFlux.cxx:856
void SetFluxSphere(TVector3 center, double radius, bool inDetCoord=true)
specification of a sphere
Definition: GDk2NuFlux.cxx:902
void User2BeamPos(const TLorentzVector &usrxyz, TLorentzVector &beamxyz) const
void User2BeamP4(const TLorentzVector &usrp4, TLorentzVector &beamp4) const
GENIE flux drivers.
TVector3 ParseTV3(const std::string &)
const XML_Char * s
Definition: expat.h:262
void SetBeamRotation(TRotation beamrot)
< beam (0,0,0) relative to user frame, beam direction in user frame
Definition: GDk2NuFlux.cxx:943
Double_t scale
Definition: plot.C:25
double enu
Definition: runWimpSim.h:113
void ParseParamSet(xmlDocPtr &, xmlNodePtr &)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
static constexpr double L
int PdgCode(void)
returns the flux neutrino pdg code
Definition: GDk2NuFlux.cxx:151
Long64_t nentries
void ParseWindowSeries(xmlDocPtr &, xmlNodePtr &)
bool GenerateNext_weighted(void)
Definition: GDk2NuFlux.cxx:255
void Beam2UserPos(const TLorentzVector &beamxyz, TLorentzVector &usrxyz) const
double UnitFromString(string u)
Definition: UnitUtils.cxx:26
void PrintCurrent(void)
print current entry from leaves
GENIE interface for uniform flux exposure iterface.
string GetXMLFilePath(string basename)
void GetFluxWindow(TVector3 &p1, TVector3 &p2, TVector3 &p3) const
3 points define a plane in beam coordinate
Definition: GDk2NuFlux.cxx:934
void AddFile(TTree *fluxtree, TTree *metatree, std::string fname)
virtual double GetTotalExposure() const
Definition: GDk2NuFlux.cxx:144
double POT_curr(void)
current average POT (RWH?)
Definition: GDk2NuFlux.cxx:565
Double_t radius
#define pINFO
Definition: Messenger.h:63
std::vector< std::string > GetFileList()
list of files currently part of chain
bool LoadParamSet(xmlDocPtr &, std::string cfg)
TVector3 Unit() const
string pname
Definition: eplot.py:33
const ana::Var wgt
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
int calcEnuWgt(const bsim::Decay &decay, const TVector3 &xyz, double &enu, double &wgt_xy)
workhorse routine
#define pWARN
Definition: Messenger.h:61
double LengthUnits(void) const
Return user units.
OStream cout
Definition: OStream.cxx:6
void User2BeamDir(const TLorentzVector &usrdir, TLorentzVector &beamdir) const
string TrimSpaces(string input)
Definition: StringUtils.cxx:24
void ParseRotSeries(xmlDocPtr &, xmlNodePtr &)
double dot(const std::vector< double > &x, const std::vector< double > &y)
Definition: dot.hpp:10
void SetFluxWindow(TVector3 p1, TVector3 p2, TVector3 p3)
3 points define a plane (by default in user coordinates)
Definition: GDk2NuFlux.cxx:865
bsim::DkMeta * dkmeta
const TLorentzVector & Position(void)
returns the flux neutrino 4-position (note: expect SI rather than physical units) ...
Definition: GDk2NuFlux.cxx:161
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
TDirectory * dir
Definition: macro.C:5
void End(void)
Definition: gXSecComp.cxx:210
bool LoadConfig(std::string cfg)
void SetBeamCenter(TVector3 beam0)
Definition: GDk2NuFlux.cxx:950
An implementation of the GENIE GFluxI interface ("flux driver") encapsulating reading/processing the ...
Definition: GDk2NuFlux.h:69
std::vector< long int > GetIntVector(std::string str)
float Mag() const
exit(0)
void Clear(Option_t *opt)
reset state variables based on opt
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
GDk2NuFluxXMLHelper(GDk2NuFlux *dk2nuFlux)
Definition: GDk2NuFlux.cxx:96
TVector3 GetBeamCenter() const
beam origin in user frame
Definition: GDk2NuFlux.cxx:971
#define pNOTICE
Definition: Messenger.h:62
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:72
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:64
Float_t e
Definition: plot.C:35
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
string GetAttribute(xmlNodePtr xml_cur, string attr_name)
void ScanForMaxWeight(void)
scan for max flux weight (before generating unweighted flux neutrinos)
Definition: GDk2NuFlux.cxx:792
TRotation GetBeamRotation() const
rotation to apply from beam->user
Definition: GDk2NuFlux.cxx:958
Float_t w
Definition: plot.C:20
const XML_Char XML_Content * model
Definition: expat.h:151
static constexpr Double_t cm
Definition: Munits.h:140
void next()
Definition: show_event.C:84
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:87
virtual TTree * GetMetaDataTree()
Definition: GDk2NuFlux.cxx:758
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
#define pDEBUG
Definition: Messenger.h:64
TVector3 AnglesToAxis(double theta, double phi, std::string units="deg")
const TLorentzVector & Momentum(void)
returns the flux neutrino 4-momentum
Definition: GDk2NuFlux.cxx:156
enum BeamMode string