Project

General

Profile

gnumi2simple_basic.C

Robert Hatcher, 09/20/2011 01:47 PM

 
1
#include "Numerical/RandomGen.h"
2
#include "FluxDrivers/GNuMIFlux.h"
3
#include "FluxDrivers/GSimpleNtpFlux.h"
4
#include "Utils/UnitUtils.h"
5

    
6
#include "TSystem.h"
7
#include "TStopwatch.h"
8
#include "TLorentzVector.h"
9
#include "TNtuple.h"
10
#include "TFile.h"
11

    
12
#include <iostream>
13
#include <iomanip>
14
#include <string>
15
#include <sstream>
16
#include <set>
17

    
18
using namespace std;
19
using namespace genie;
20
using namespace genie::flux;
21

    
22
// forward declaration
23
void encode_transform(flux::GNuMIFlux* gnumi, string& rot, string& pos);
24

    
25
// main routine
26
void gnumi2simple_basic(string fnameout="fluxntgenie.root",
27
                        string fluxpatt="",
28
                        string cfg="MINOS-Near", 
29
                        long int nentries=0, 
30
                        double pots=1.0e30,
31
                        double enumin=0.0,
32
                        bool doaux=true,
33
                        string srcpatt="")
34
{
35
  
36
  cout << "Creating:    " << fnameout << endl;
37
  cout << "Input files: " << fluxpatt << endl;
38
  cout << "Config:      " << cfg << endl;
39
  cout << "NEntries:    " << nentries << endl;
40
  cout << "POTs:        " << pots << endl;
41
  if ( enumin > 0.0 ) 
42
    cout << "EnuMin:      " << enumin << endl;
43
  cout << "AuxTree:     " << doaux << endl;
44
  if ( srcpatt != "" )
45
    cout << "SrcPatt:     " << srcpatt << endl;
46
  
47
  string fluxfname(gSystem->ExpandPathName(fluxpatt.c_str()));
48
  flux::GNuMIFlux* gnumi = new GNuMIFlux();
49
  gnumi->LoadBeamSimData(fluxfname,cfg);
50
  gnumi->SetEntryReuse(1); // don't reuse entries when reformatting
51

    
52
  //gnumi->PrintConfig();
53
  //cout << " change to cm:" << endl;
54
  //double cm_unit = genie::utils::units::UnitFromString("cm");
55
  //gnumi->SetLengthUnits(cm_unit);
56
  //gnumi->PrintConfig();
57
  //cout << " change to m:" << endl;
58
  ///////RWH: the following 2 lines should work ...
59
  //double m_unit = genie::utils::units::UnitFromString("m");
60
  //gnumi->SetLengthUnits(m_unit);
61
  gnumi->PrintConfig();
62
  //double scaleusr = 100.;
63

    
64
  gnumi->SetEntryReuse(1);     // reuse entries
65
  gnumi->SetUpstreamZ(-3e38);  // leave ray on flux window
66
  
67
  GFluxI* fdriver = dynamic_cast<GFluxI*>(gnumi);
68

    
69
  if (nentries == 0) nentries = 2147483647;
70

    
71
  // so as not to include scan time generate 1 nu
72
  fdriver->GenerateNext();
73

    
74
  TFile* file = TFile::Open(fnameout.c_str(),"RECREATE");
75
  TTree* fluxntp = new TTree("flux","a simple flux n-tuple");
76
  TTree* metantp = new TTree("meta","metadata for flux n-tuple");
77
  genie::flux::GSimpleNtpEntry* fentry = new genie::flux::GSimpleNtpEntry;
78
  genie::flux::GSimpleNtpNuMI*  fnumi  = new genie::flux::GSimpleNtpNuMI;
79
  genie::flux::GSimpleNtpAux*   faux   = new genie::flux::GSimpleNtpAux;
80
  genie::flux::GSimpleNtpMeta*  fmeta  = new genie::flux::GSimpleNtpMeta;
81
  fluxntp->Branch("entry",&fentry);
82
  fluxntp->Branch("numi",&fnumi);
83
  if (doaux) fluxntp->Branch("aux",&faux);
84
  metantp->Branch("meta",&fmeta);
85

    
86
  TLorentzVector p4u, x4u;
87
  long int ngen = 0;
88

    
89
  set<int> pdglist;
90
  double maxe = 0;
91
  double minwgt = +1.0e10;
92
  double maxwgt = -1.0e10;
93

    
94
  UInt_t metakey = TString::Hash(fnameout.c_str(),strlen(fnameout.c_str()));
95
  cout << "metakey " << metakey << endl;
96
  // ensure that we don't get smashed by UInt_t vs Int_t
97
  metakey &= 0x7FFFFFFF;
98
  cout << "metakey " << metakey << " after 0x7FFFFFFF" << endl;
99
  cout << "=========================== Start " << endl;
100

    
101
  TStopwatch sw;
102
  sw.Start();
103
  int nok = 0, nwrite = 0;
104
  while ( nok < nentries && gnumi->UsedPOTs() < pots ) {
105
    fdriver->GenerateNext();
106
    ngen++;
107
    fentry->Reset();
108
    fnumi->Reset();
109
    faux->Reset();
110

    
111
    fentry->metakey = metakey;
112
    fentry->pdg     = fdriver->PdgCode();
113
    fentry->wgt     = gnumi->Weight();
114
    x4u     = fdriver->Position();
115
    fentry->vtxx    = x4u.X();
116
    fentry->vtxy    = x4u.Y();
117
    fentry->vtxz    = x4u.Z();
118
    fentry->dist    = gnumi->GetDecayDist();
119
    p4u     = fdriver->Momentum();
120
    fentry->px      = p4u.Px();
121
    fentry->py      = p4u.Py();
122
    fentry->pz      = p4u.Pz();
123
    fentry->E       = p4u.E();
124

    
125
    fnumi->run      = gnumi->PassThroughInfo().run;
126
    fnumi->evtno    = gnumi->PassThroughInfo().evtno;
127
    fnumi->entryno  = gnumi->GetEntryNumber();
128

    
129
    fnumi->tpx      = gnumi->PassThroughInfo().tpx;
130
    fnumi->tpy      = gnumi->PassThroughInfo().tpy;
131
    fnumi->tpz      = gnumi->PassThroughInfo().tpz;
132
    fnumi->vx       = gnumi->PassThroughInfo().vx;
133
    fnumi->vy       = gnumi->PassThroughInfo().vy;
134
    fnumi->vz       = gnumi->PassThroughInfo().vz;
135

    
136
    fnumi->ndecay   = gnumi->PassThroughInfo().ndecay;
137
    fnumi->ppmedium = gnumi->PassThroughInfo().ppmedium;
138
    fnumi->tptype   = gnumi->PassThroughInfo().tptype;
139

    
140

    
141
    if ( doaux ) {
142
      faux->auxint.push_back(gnumi->PassThroughInfo().ndecay);
143
      faux->auxint.push_back(gnumi->PassThroughInfo().ptype);
144
      faux->auxint.push_back(gnumi->PassThroughInfo().tgen);
145
      faux->auxdbl.push_back(gnumi->PassThroughInfo().fgXYWgt);
146
      faux->auxdbl.push_back(gnumi->PassThroughInfo().nimpwt);
147
    }
148

    
149
    if ( fentry->E > enumin ) ++nok;
150
    fluxntp->Fill();
151
    ++nwrite;
152

    
153
    // accumulate meta data
154
    pdglist.insert(fentry->pdg);
155
    minwgt = TMath::Min(minwgt,fentry->wgt);
156
    maxwgt = TMath::Max(maxwgt,fentry->wgt);
157
    maxe   = TMath::Max(maxe,fentry->E);
158

    
159
  }
160
  cout << "=========================== Complete " << endl;
161

    
162
  fmeta->pdglist.clear();
163
  set<int>::const_iterator setitr = pdglist.begin();
164
  for ( ; setitr != pdglist.end(); ++setitr)  fmeta->pdglist.push_back(*setitr);
165
 
166
  fmeta->maxEnergy = maxe;
167
  fmeta->minWgt    = minwgt;
168
  fmeta->maxWgt    = maxwgt;
169
  fmeta->protons   = gnumi->UsedPOTs();
170
  TVector3 p0, p1, p2;
171
  gnumi->GetFluxWindow(p0,p1,p2);
172
  TVector3 d1 = p1 - p0;
173
  TVector3 d2 = p2 - p0;
174
  fmeta->windowBase[0] = p0.X();
175
  fmeta->windowBase[1] = p0.Y();
176
  fmeta->windowBase[2] = p0.Z();
177
  fmeta->windowDir1[0] = d1.X();
178
  fmeta->windowDir1[1] = d1.Y();
179
  fmeta->windowDir1[2] = d1.Z();
180
  fmeta->windowDir2[0] = d2.X();
181
  fmeta->windowDir2[1] = d2.Y();
182
  fmeta->windowDir2[2] = d2.Z();
183
  if ( doaux ) {
184
    fmeta->auxintname.push_back("ndecay");
185
    fmeta->auxintname.push_back("ptype");
186
    fmeta->auxintname.push_back("tgen");
187
    fmeta->auxdblname.push_back("fgXYWgt");
188
    fmeta->auxdblname.push_back("nimpwt");
189
  }
190

    
191
  std::string rotinfo, posinfo;
192
  encode_transform(gnumi,rotinfo,posinfo);
193
  fmeta->infiles.push_back(rotinfo);
194
  fmeta->infiles.push_back(posinfo);
195

    
196
  if ( srcpatt != "" ) {
197
    string srcstring = "srcpatt=" + srcpatt;
198
    fmeta->infiles.push_back(srcstring);
199
  }
200
  string fname = "flux_pattern=" + fluxfname;
201
  fmeta->infiles.push_back(fname); // should get this expanded from gnumi
202
  vector<string> flist = gnumi->GetFileList();
203
  for (size_t i = 0; i < flist.size(); ++i) {
204
    fname = "flux_infile=" + flist[i];
205
    fmeta->infiles.push_back(fname);
206
  }
207

    
208
  fmeta->seed    = RandomGen::Instance()->GetSeed();
209
  fmeta->metakey = metakey;
210

    
211
  metantp->Fill();
212

    
213
  sw.Stop();
214
  cout << "Generated " << nwrite << " (" << nok << " w/ E >" << enumin << " ) "
215
       << " ( request " << nentries << " ) "
216
       << endl
217
       << gnumi->UsedPOTs() << " POTs " << " ( request " << pots << " )"
218
       << endl
219
       << " pulled NFluxNeutrinos " << gnumi->NFluxNeutrinos()
220
       << endl
221
       << "Time to generate: " << endl;
222
  sw.Print();
223
  
224
  cout << "===================================================" << endl;
225

    
226
  cout << "Last GSimpleNtpEntry: " << endl
227
       << *fentry
228
       << endl
229
       << "Last GSimpleNtpMeta: " << endl
230
       << *fmeta
231
       << endl;
232

    
233
  file->cd();
234

    
235
  // write meta as simple object at top of file
236
  //fmeta->Write();
237

    
238
  // write ntuples out
239
  fluxntp->Write();
240
  metantp->Write();
241
  file->Close();
242

    
243
  cout << endl << endl;
244
  gnumi->PrintConfig();
245
}
246

    
247
void encode_transform(flux::GNuMIFlux* gnumi, string& rot, string& pos)
248
{
249
  TRotation rot3x3 = gnumi->GetBeamRotation();
250
  TVector3  posv3  = gnumi->GetBeamCenter();
251

    
252
  std::ostringstream rotstr;
253
  std::ostringstream posstr;
254

    
255
  const int p=15, w=20;
256
  
257
  rotstr << std::setprecision(p);
258
  rotstr << "<beamdir type=\"newxyz\"> \n"
259
         << "[ "
260
         << std::setw(w) << rot3x3.XX() << " "
261
         << std::setw(w) << rot3x3.XY() << " "
262
         << std::setw(w) << rot3x3.XZ() << " ]\n"
263
         << "[ "
264
         << std::setw(w) << rot3x3.YX() << " "
265
         << std::setw(w) << rot3x3.YY() << " "
266
         << std::setw(w) << rot3x3.YZ() << " ]\n"
267
         << "[ "
268
         << std::setw(w) << rot3x3.ZX() << " "
269
         << std::setw(w) << rot3x3.ZY() << " "
270
         << std::setw(w) << rot3x3.ZZ() << " ]\n"
271
         << "    </beamdir>" << std::endl;
272
  posstr << std::setprecision(p);
273
  posstr << "<beampos> ( "
274
         << std::setw(w) << posv3.X() << " "
275
         << std::setw(w) << posv3.Y() << " "
276
         << std::setw(w) << posv3.Z() << " ) "
277
         << "</beampos>" << std::endl;
278

    
279
  rot = rotstr.str();
280
  pos = posstr.str();
281
}
282