Project

General

Profile

myfakefluxgen.C

example for creating gsimple flux file from nothing - Robert Hatcher, 04/03/2018 04:33 PM

 
1
// a sample script for creating/filling a GSimpleNtpFlux file
2
// use:   $ genie -b -q myfakefluxgen.C+
3
// 
4
// R. Hatcher <rhatcher@fnal.gov> 2014-05-07
5

    
6
#if !defined(__CINT__) || defined(__MAKECINT__)
7
// hide from CINT, but not ACLiC
8
  #include "Numerical/RandomGen.h"
9
  #include "FluxDrivers/GSimpleNtpFlux.h"
10
  #include "Utils/UnitUtils.h"
11
  #include "Utils/AppInit.h"
12
  #include "Numerical/RandomGen.h"
13

    
14
  #include "TSystem.h"
15
  #include "TStopwatch.h"
16
  #include "TLorentzVector.h"
17
  #include "TNtuple.h"
18
  #include "TFile.h"
19

    
20
  #include <iostream>
21
  #include <iomanip>
22
  #include <string>
23
  #include <sstream>
24
  #include <set>
25
  #include <stdlib.h>  // for strtol, putenv, unsetenv
26
#endif 
27

    
28
// my own laziness:
29
using namespace std;
30
using namespace genie;
31
using namespace genie::flux;
32

    
33
// forward declaration
34
void getInfo(TLorentzVector& p4u, TLorentzVector& x4u,
35
             int& pdg, double& wgt, double& decayDist);
36

    
37
/// main routine
38
void myfakefluxgen(std::string fnameout="myfakeflux.gsimple.root",
39
                   long int nentries = 100000,
40
                   long int seed = 12345) 
41
{
42
  /// Generate a GSimpleNtpFlux file
43
  /// ...out of nothing
44

    
45
  genie::utils::app_init::RandGen(seed);
46

    
47
  /// create the file, trees and branches
48
  TFile* file = TFile::Open(fnameout.c_str(),"RECREATE");
49
  TTree* fluxntp = new TTree("flux","a simple flux n-tuple");
50
  TTree* metantp = new TTree("meta","metadata for flux n-tuple");
51
  genie::flux::GSimpleNtpEntry* fentry = new genie::flux::GSimpleNtpEntry;
52
  genie::flux::GSimpleNtpMeta*  fmeta  = new genie::flux::GSimpleNtpMeta;
53
  fluxntp->Branch("entry",&fentry);
54
  metantp->Branch("meta",&fmeta);
55
  
56
  TLorentzVector p4u, x4u;
57
  double wgt, decayDist;
58
  int pdg;
59

    
60
  double POTs=0, minwgt=1.0e30, maxwgt=0, maxe=0;
61

    
62
  /// use metakey to associate entries w/ metadata
63
  UInt_t metakey = TString::Hash(fnameout.c_str(),strlen(fnameout.c_str()));
64
  
65
  /// generate & store entries
66
  long int ngen = 0;
67
  while ( ngen < nentries ) {
68
    fentry->Reset();
69
    ngen++;
70
    
71
    // make up an entry ... 
72
    // p4u and x4u are the 4-momentum and 4-position in the user
73
    // (i.e. detector) frame
74
    getInfo(p4u,x4u,pdg,wgt,decayDist);
75

    
76
    // in our fake setup each entry represents 
77
    //  3.1415 protons-on-target on average
78
    POTs += 3.1415;
79

    
80
    // fill the entry in the tree
81
    fentry->metakey = metakey;
82
    fentry->pdg     = pdg;
83
    fentry->wgt     = wgt;
84

    
85
    // flux positions are expected to be in SI units (i.e. meters)
86
    fentry->vtxx    = x4u.X();
87
    fentry->vtxy    = x4u.Y();
88
    fentry->vtxz    = x4u.Z();
89
    fentry->dist    = decayDist;
90
    
91
    fentry->px      = p4u.Px();
92
    fentry->py      = p4u.Py();
93
    fentry->pz      = p4u.Pz();
94
    fentry->E       = p4u.E();
95
    
96
    fluxntp->Fill();
97
    
98
    // accumulate meta data
99
    fmeta->AddFlavor(fentry->pdg);
100
    minwgt = TMath::Min(minwgt,fentry->wgt);
101
    maxwgt = TMath::Max(maxwgt,fentry->wgt);
102
    maxe   = TMath::Max(maxe,fentry->E);
103
    
104
  }
105
  
106
  fmeta->maxEnergy = maxe;
107
  fmeta->minWgt    = minwgt;
108
  fmeta->maxWgt    = maxwgt;
109
  fmeta->protons   = POTs;
110
  
111
  fmeta->seed    = seed;
112
  fmeta->metakey = metakey;
113
  metantp->Fill();
114
  
115
  // flush/write ntuples out
116
  fluxntp->Write();
117
  metantp->Write();
118
  file->Close();
119
  
120
}
121

    
122
// forward declaration
123
void getInfo(TLorentzVector& p4u, TLorentzVector& x4u,
124
             int& pdg, double& wgt, double& decayDist)
125
{
126
  
127
  // GENIE's pseudo-random number central station
128
  RandomGen* rg  = RandomGen::Instance();
129

    
130
  // here's the pseudo-random # engine we'll use
131
  // make sure to use a _reference_ here
132
  // otherwise you're creating a copy and it will not reflect changing
133
  // internal state (and thus you'll get the same values each call to getInfo)
134
  // Ah, CINT barfs at this use of a reference ... don't know why 
135
  // ACLiC compiles it fine ... whatever!  just use this expression everywhere
136
  // TRandom3& rndm = rg->RndFlux();
137

    
138
  // pick a flavor
139
  double r = rg->RndFlux().Uniform(1.0);
140
  pdg = 14;  // nu_mu
141
  if ( r > 0.75 ) pdg = -14;  // nu_mu_bar
142
  if ( r > 0.95 ) pdg =  12;  // nu_e
143
  if ( r > 0.99 ) pdg = -12;  // nu_e_bar
144

    
145
  // energy spectrum ... how about Landau
146
  // make sure it sensible, though ... no zero or negative energies
147
  // nothing crazy high either (cut it at 100GeV)
148
  double e = 0;
149
  while ( e <= 0 || e > 100 ) {
150
    e = rg->RndFlux().Landau(2.0,1.0);
151
  }
152
  // just have it come face on +z direction
153
  p4u = TLorentzVector(0,0,e,e);
154

    
155
  // flux coming off a plane at z=0
156
  // guassian intensity in x-y, centered at x=0.5 w/ sigmax=4, sigmay=3
157
  // but cut so nothing |x|>15 or |y|>15
158
  bool okay_xy = false;
159
  double a, b, x, y, z=0, t=0;
160
  while ( ! okay_xy ) {
161
    rg->RndFlux().Rannor(a,b);
162
    x = 0.5 + a*4.0;
163
    y = b*3.0;
164
    if ( fabs(x) < 15 || fabs(y) < 15 ) okay_xy = true;
165
  }
166
  x4u = TLorentzVector(x,y,z,t);
167

    
168
  // decayDist is distance from decay of particle giving the neutrino
169
  // to the x4u position ... for this test just put 735e3
170
  decayDist = 735.e3;
171

    
172
  // unweighted rays are easiest
173
  wgt = 1.0;
174

    
175
}
176