10 #include "GENIE/Framework/Numerical/RandomGen.h" 11 #include "GENIE/Tools/Flux/GSimpleNtpFlux.h" 12 #include "GENIE/Framework/Utils/UnitUtils.h" 15 #include "TStopwatch.h" 16 #include "TLorentzVector.h" 32 using namespace genie;
61 TLorentzVector fmb4 = TLorentzVector(entry->
px, entry->
py, entry->
pz,
69 double k =
sqrt(1+(entry->
E*entry->
E-fmb4.E()*fmb4.E())/fmb4.Vect().Mag2());
75 static double trash = 0;
76 if ( numi ) trash = numi->
entryno;
77 if ( aux ) trash = aux->
auxint.size();
78 if ( trash == -1 )
cout <<
"trash was -1" <<
endl;
83 meta->
infiles.push_back(
"THIS IS MDC WARPED FLUX FROM:");
84 meta->
infiles.push_back(fnamein);
86 meta->
infiles.push_back(
"NO ACTUAL WARPING APPLIED");
88 string msg =
"USING WARP FUNCTION CODENAMED WHISKEY-TANGO-FOXTROT";
105 return fcn->GetRandom(emin, emax);
111 auto profile = [](
double y,
double*
p){
112 return p[0]+
p[1]*TMath::Power(y, -
p[2])+
p[3]*TMath::Power(y, -
p[4])+
p[5]*
y;
114 double thresholdE = 0.1;
115 if(energy > thresholdE)
116 return profile(energy, par);
118 return profile(100., par);
131 string fnamein=
"gsimple_input.root",
140 fcn->SetParameters(-2.43864
e+02, 4.42211
e+02, 3.00555
e+00, 1.23288
e+04, 8.95128
e-01, 2.38448
e+00);
144 string fnameinlong(gSystem->ExpandPathName(fnamein.c_str()));
146 TFile*
fin = TFile::Open(fnameinlong.c_str(),
"READONLY");
147 TTree* etreein = (TTree*)fin->Get(
"flux");
148 TTree* mtreein = (TTree*)fin->Get(
"meta");
154 long int nentries = etreein->GetEntries();
156 int sba_status[4] = { -999, -999, -999, -999 };
157 sba_status[0] = etreein->SetBranchAddress(
"entry",&entry_in);
158 sba_status[1] = etreein->SetBranchAddress(
"numi",&numi_in);
159 sba_status[2] = etreein->SetBranchAddress(
"aux",&aux_in);
160 sba_status[3] = mtreein->SetBranchAddress(
"meta",&meta_in);
161 cout <<
"SetBranchAddress results " 162 << sba_status[0] <<
"," 163 << sba_status[1] <<
"," 164 << sba_status[2] <<
"," 167 bool donumi = ( sba_status[1] == 0 );
168 bool doaux = ( sba_status[2] == 0 );
170 int nindices = mtreein->BuildIndex(
"metakey");
171 cout <<
"saw " << nindices <<
" metakey indices" <<
endl;
173 cout <<
"Creating: " << fnameout <<
endl;
174 cout <<
"Input file: " << fnamein <<
endl;
175 cout <<
"Branches: entry " 176 << ( (doaux) ?
"aux ":
"" )
177 << ( (donumi) ?
"numi ":
"" )
181 TFile*
fout = TFile::Open(fnameout.c_str(),
"RECREATE");
182 TTree* fluxntp =
new TTree(
"flux",
"a simple flux n-tuple");
183 TTree* metantp =
new TTree(
"meta",
"metadata for flux n-tuple");
189 fluxntp->Branch(
"entry",&fentry);
190 if ( doaux ) fluxntp->Branch(
"aux",&faux);
191 if ( donumi ) fluxntp->Branch(
"numi",&fnumi);
192 metantp->Branch(
"meta",&fmeta);
194 cout <<
"=========================== Start " <<
endl;
199 const double large = 1.0e300;
200 double minwgt = large, maxwgt = -large, maxenergy = -large;
201 std::set<int> gFlavors;
203 for (
long int ientry = 0; ientry <
nentries; ++ientry ) {
212 etreein->GetEntry(ientry);
213 UInt_t metakey = entry_in->
metakey;
214 if ( fmeta->
metakey != metakey ) {
216 int nbmeta = mtreein->GetEntryWithIndex(metakey);
217 cout <<
"on entry " << ientry <<
" fetched metadata w/ key " 218 << metakey <<
" read " << nbmeta <<
" bytes" <<
endl;
234 if ( fentry->
wgt < minwgt ) minwgt = fentry->
wgt;
235 if ( fentry->
wgt > maxwgt ) maxwgt = fentry->
wgt;
237 if ( fentry->
E > maxenergy ) maxenergy = fentry->
E;
239 gFlavors.insert(fentry->
pdg);
254 cout <<
"metantp->Fill() " << *fmeta <<
endl;
279 cout <<
"metantp->Fill() " << *fmeta <<
endl;
283 cout <<
"=========================== Complete " <<
endl;
286 cout <<
"Generated " << nentries
288 <<
"Time to generate: " <<
endl;
304 std::set<int>::const_iterator flitr = gFlavors.begin();
305 std::vector<Int_t>& flvlist = meta->
pdglist;
306 for ( ; flitr != gFlavors.end(); ++flitr ) {
307 int seen_pdg = *flitr;
308 std::vector<Int_t>::iterator knwitr =
309 std::find(flvlist.begin(),flvlist.end(),seen_pdg);
310 bool known_pdg = ( knwitr != flvlist.end() );
311 if ( ! known_pdg ) meta->
pdglist.push_back(seen_pdg);
Double_t E
energy in lab frame
void warp_gsimple(string fnameout="gsimple_output.root", string fnamein="gsimple_input.root", bool dowarp=true)
THE MAIN GENIE PROJECT NAMESPACE
Double_t px
x momentum in lab frame
std::vector< Int_t > auxint
additional ints associated w/ entry
void update_flavors(GSimpleNtpMeta *meta, std::set< int > &gFlavors)
void warp_entry(GSimpleNtpEntry *entry, GSimpleNtpAux *aux, GSimpleNtpNuMI *numi)
void warp_meta(GSimpleNtpMeta *meta, string fnamein)
double pick_energy(double, double)
double fit_fcn(double *x, double *par)
Double_t pz
z momentum in lab frame
UInt_t metakey
key to meta data
Double_t py
y momentum in lab frame