Project

General

Profile

myfakefluxgen3.C

R-3_X_Y & ROOT6 - Robert Hatcher, 04/03/2018 04:59 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
// only for R-3_
7
#include "Framework/Conventions/GVersion.h"
8
#include "Framework/Numerical/RandomGen.h"
9
#include "Tools/Flux/GSimpleNtpFlux.h"
10
#include "Framework/Utils/UnitUtils.h"
11
#include "Framework/Utils/AppInit.h"
12

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

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

    
26
// my own laziness:
27
using namespace std;
28
using namespace genie;
29
using namespace genie::flux;
30

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

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

    
43
  genie::utils::app_init::RandGen(seed);
44

    
45
  /// create the file, trees and branches
46
  TFile* file = TFile::Open(fnameout.c_str(),"RECREATE");
47
  TTree* fluxntp = new TTree("flux","a simple flux n-tuple");
48
  TTree* metantp = new TTree("meta","metadata for flux n-tuple");
49
  genie::flux::GSimpleNtpEntry* fentry = new genie::flux::GSimpleNtpEntry;
50
  genie::flux::GSimpleNtpMeta*  fmeta  = new genie::flux::GSimpleNtpMeta;
51
  fluxntp->Branch("entry",&fentry);
52
  metantp->Branch("meta",&fmeta);
53

    
54
  TLorentzVector p4u, x4u;
55
  double wgt, decayDist;
56
  int pdg;
57

    
58
  double POTs=0, minwgt=1.0e30, maxwgt=0, maxe=0;
59

    
60
  /// use metakey to associate entries w/ metadata
61
  UInt_t metakey = TString::Hash(fnameout.c_str(),strlen(fnameout.c_str()));
62

    
63
  /// generate & store entries
64
  long int ngen = 0;
65
  while ( ngen < nentries ) {
66
    fentry->Reset();
67
    ngen++;
68

    
69
    // make up an entry ...
70
    // p4u and x4u are the 4-momentum and 4-position in the user
71
    // (i.e. detector) frame
72
    getInfo(p4u,x4u,pdg,wgt,decayDist);
73

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

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

    
83
    // flux positions are expected to be in SI units (i.e. meters)
84
    fentry->vtxx    = x4u.X();
85
    fentry->vtxy    = x4u.Y();
86
    fentry->vtxz    = x4u.Z();
87
    fentry->dist    = decayDist;
88

    
89
    fentry->px      = p4u.Px();
90
    fentry->py      = p4u.Py();
91
    fentry->pz      = p4u.Pz();
92
    fentry->E       = p4u.E();
93

    
94
    fluxntp->Fill();
95

    
96
    // accumulate meta data
97
    fmeta->AddFlavor(fentry->pdg);
98
    minwgt = TMath::Min(minwgt,fentry->wgt);
99
    maxwgt = TMath::Max(maxwgt,fentry->wgt);
100
    maxe   = TMath::Max(maxe,fentry->E);
101

    
102
  }
103

    
104
  fmeta->maxEnergy = maxe;
105
  fmeta->minWgt    = minwgt;
106
  fmeta->maxWgt    = maxwgt;
107
  fmeta->protons   = POTs;
108

    
109
  fmeta->seed    = seed;
110
  fmeta->metakey = metakey;
111
  metantp->Fill();
112

    
113
  // flush/write ntuples out
114
  fluxntp->Write();
115
  metantp->Write();
116
  file->Close();
117

    
118
}
119

    
120
// forward declaration
121
void getInfo(TLorentzVector& p4u, TLorentzVector& x4u,
122
             int& pdg, double& wgt, double& decayDist)
123
{
124

    
125
  // GENIE's pseudo-random number central station
126
  RandomGen* rg  = RandomGen::Instance();
127

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

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

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

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

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

    
170
  // unweighted rays are easiest
171
  wgt = 1.0;
172

    
173
}