runWimpSim.cpp
Go to the documentation of this file.
1 
2 #include "runWimpSim.h"
3 
4 
5 void setup_wimpann();
6 void get_wimpann_info(unsigned int&, double&, unsigned int&,
7  bool&, bool&, std::string&,
9  std::string&);
10 void get_wimpevent_info();
11 
12 int main(int argc, char* argv[]) {
13 
14  if (argc < 2) {
15  std::cout << "usage: runWimpSim flux_out_filename" << std::endl;
16  return 1;
17  }
18 
19  if (!getenv("WIMPSIM_DIR")) {
20  std::cout << "wimpsim must be set up before running this program." << std::endl;
21  std::cout << "do: setup wimpsim v3_05_00" << std::endl;
22  return 1;
23  }
24 
25  std::string fnameout(argv[1]);
26 
27 //*********************** WIMPANN *************************
28  setup_wimpann();
29  // variables
30  unsigned int ann_chan;
31  double cm_energy, mwimp;
32  unsigned int nann;
33  bool wimpann_outp, use_earth_ann, wimpann_sum;
34  std::string file_s, file_e, file_info, file_sum_s, file_sum_e;
35  get_wimpann_info(ann_chan, cm_energy, nann, wimpann_outp,
36  use_earth_ann, file_s, file_e, file_info,
37  wimpann_sum, file_sum_s, file_sum_e);
38  mwimp = cm_energy / 2.0;
39 
40  // setup pythia seed
41  std::stringstream ss;
42  ss << "MRPY(1)=" << wimppara_.seed;
43  std::cout << "Using Pythia seed " << ss.str() << std::endl;
44  //std::cout << ss.str().length() << std::endl;
45  char cseed[20];
46  for (size_t i=0; i<ss.str().length() && i<20; i++) {
47  cseed[i] = ss.str().c_str()[i];
48  }
49  for (size_t i=ss.str().length(); i<20; i++) {
50  cseed[i] = ' ';
51  }
52  //std::cout << cseed << std::endl;
53  //cseed[20] = ' ';
54  size_t len = 20;
55  //pass hidden array length variable
56  pygive_((char*) cseed, len);
57 
58  // Initialize neutrino oscillation routines
59  nuoscinit_();
60 
61  // Modify Pythia
62  modpyt_();
63 
64  //int nnustot=0, nnuetot=0;
65  neutrino_.enumax = mwimp;
66  neutrino_.mwimp = mwimp;
67 
68 //*********************** WIMPEVENT *************************
69  weinit_();
70  std::cout << "wimpevent setup complete." << std::endl;
72 
73 //*********************** GENIE *****************************
74 
75  //bool doaux = false;
76 
77  long int gseed = (long int) wimppara_.seed;
79 
80  TFile* file = TFile::Open(fnameout.c_str(),"RECREATE");
81  TTree* fluxntp = new TTree("flux","a simple flux n-tuple");
82  TTree* metantp = new TTree("meta","metadata for flux n-tuple");
87  fluxntp->Branch("entry",&fentry);
88  fluxntp->Branch("numi",&fnumi);
89  fluxntp->Branch("aux",&faux);
90  metantp->Branch("meta",&fmeta);
91 
92  //long int ngen = 0;
93 
94  std::set<int> pdglist;
95  double maxe = 0;
96  double minwgt = +1.0e10;
97  double maxwgt = -1.0e10;
98 
99  UInt_t metakey = TString::Hash(fnameout.c_str(),strlen(fnameout.c_str()));
100  std::cout << "metakey " << metakey << std::endl;
101  // ensure that we don't get smashed by UInt_t vs Int_t
102  metakey &= 0x7FFFFFFF;
103  std::cout << "GENIE metakey " << metakey << " after 0x7FFFFFFF" << std::endl;
104 
105  TRandom3 *myRand = new TRandom3(0);
106 
107 //********************** Simulation *************************
108  std::cout << "=========================== Start " << std::endl;
109 
110  TStopwatch sw;
111  sw.Start();
112 
113  long int nann_sofar = 0;
114 
115  unsigned nok = 0, nwrite = 0;
116  // simulate in batches of 100,000
117  for (unsigned int i=0; i<=nann/100000; i++) { // removed nok < nFluxEntries
118  unsigned int todo = (nann-i*100000<100000) ? nann-i*100000 : 100000;
119  genevt_(&ann_chan, &todo);
120  nann_sofar += todo;
121 
122  // now neutrino information is given in nus_
123  // run the current set of neutrinos through wimpevent
124  for (unsigned int j=0; j<nus_.nnus; j++) {
125  if (!!nus_.nusi[j][6]&0x1) { // not absorbed
126 
127  nwrite++;
128 
129  double enu = nus_.nusr[j][1];
130  //std::cout << "\n" << enu << std::endl;
131  weevent_.enu = enu;
132  if (enu < enumin) continue; // energy too low
133 
134  unsigned nut = nutype(nus_.nusi[j][4]);
135  //std::cout << std::endl << nus_.nusi[j][4] << std::endl;
136  //std::cout << nut << std::endl;
137  weevent_.nut = nut;
138 
139  // see wimpann/writeevs
140  double a_e, a_mu, a_tau, phi_e, phi_mu, phi_tau;
141  a_e = zamp(nus_.nusc[j][0]);
142  phi_e = zarg(nus_.nusc[j][0]);
143  a_mu = zamp(nus_.nusc[j][1]);
144  phi_mu = zarg(nus_.nusc[j][1]);
145  a_tau = zamp(nus_.nusc[j][2]);
146  phi_tau = zarg(nus_.nusc[j][2]);
147 
148  phi_mu = phi_mu-phi_e;
149  phi_tau = phi_tau-phi_e;
150 
151  //std::cout << "a_e: " << a_e << std::endl;
152  //std::cout << "a_mu: " << a_mu << std::endl;
153  //std::cout << "phi_mu: " << phi_mu << std::endl;
154 
155  // now setup event info for wimpevent
156  weevent_.nustate[0].r = a_e;
157  weevent_.nustate[0].i = 0;
158  weevent_.nustate[1].r = a_mu*cos(phi_mu);
159  weevent_.nustate[1].i = a_mu*sin(phi_mu);
160  weevent_.nustate[2].r = a_tau*cos(phi_tau);
161  weevent_.nustate[2].i = a_tau*sin(phi_tau);
162 
163  //std::cout << "mu.r before: " << weevent_.nustate[1].r << std::endl;
164  //std::cout << "a_mu before: " << zamp(weevent_.nustate[1]) << std::endl;
165 
166  // propagate neutrino to detector
167  nupropsu_();
168 
169  //std::cout << "a_mu after: " << zamp(weevent_.nustate[1]) << std::endl;
170  // project onto a neutrino state at the detector
171  // see wimpann/nuproject
172  double pe = weevent_.nustate[0].r*weevent_.nustate[0].r
173  + weevent_.nustate[0].i*weevent_.nustate[0].i,
174  pmu = weevent_.nustate[1].r*weevent_.nustate[1].r
175  + weevent_.nustate[1].i*weevent_.nustate[1].i;
176  //std::cout << "pe: " << pe << std::endl;
177  //std::cout << "pmu: " << pmu << std::endl;
178  int fam;
179  //int temp = 0;
180  double r = myRand->Rndm();
181  //std::cout << "r: " << r << std::endl;
182  if (r < pe) fam = 1;
183  else if (r < (pe+pmu)) fam = 2;
184  else fam = 3;
185  weevent_.nui = fam*2-2+weevent_.nut;
186  //std::cout << "nut: " << weevent_.nut << std::endl;
187  //std::cout << "nui: " << weevent_.nui << std::endl;
188  if (weevent_.nui != 3 && weevent_.nui != 4)
189  continue; // only want nu_mu or nu_mu~
190 
191 
192  // save event into gnumiflux
193  fentry->Reset();
194  fnumi->Reset();
195  faux->Reset();
196 
197  int pdg = get_pdg(weevent_.nui);
198  double wgt, x, y, z, px, py, pz;//, t;
199  wgt = 1.0/(4*PI*weevent_.dist*weevent_.dist)/nann;
200  x = (X_MAX - X_MIN)*myRand->Rndm()+X_MIN;
201  y = (Y_MAX - Y_MIN)*myRand->Rndm()+Y_MIN;
202  z = (Z_MAX - Z_MIN)*myRand->Rndm()+Z_MIN;
203  //t = (T_MAX - T_MIN)*myRand->Rndm()+T_MIN;
204  // GENIE wants units in meters
205  x /= 100;
206  y /= 100;
207  z /= 100;
208 
209  double part_az = weevent_.nuaz;
210  double part_el = weevent_.nuel;
211  double mjd = weevent_.mjd;
212 
213 //Original code
214 /*
215  if (part_el*180/PI < 10) continue; // only want upward-going
216 
217  //coordinate transformations taken from MoonShadowTrigger
218  double phi = part_az-(NOVA_AZIMUTH*PI/180.0); // 332o03'58.07169"
219  px = -1.0 * enu * cos(part_el) * sin(phi);
220  py = enu * sin(part_el);
221  pz = enu * cos(part_el) * cos(phi);
222 */
223 
224 //New code
225  //coordinate transformations taken from MoonShadowTrigger
226  double phi = part_az-(NOVA_AZIMUTH*PI/180.0); // 332o03'58.07169"
227 
228  //part_el and phi are the angles where the particles came from!
229  /* px = -1.0 * enu * cos(-part_el) * sin(phi+180);
230  py = enu * sin(-part_el);
231  pz = enu * cos(-part_el) * cos(phi+180);*/
232 
233 
234  px = enu * cos(part_el) * sin(phi);
235  py = - enu * sin(part_el);
236  pz = - enu * cos(part_el) * cos(phi);
237 
238  // if (py < 0) continue; //only up mu
239 
240  //x ... west
241  //y ... up
242  //z ... north
243 
244  fentry->metakey = metakey;
245  fentry->pdg = pdg;
246  fentry->wgt = wgt;
247  fentry->vtxx = x;
248  fentry->vtxy = y;
249  fentry->vtxz = z;
250  fentry->dist = 0;//gnumi->GetDecayDist();
251  fentry->px = px;
252  fentry->py = py;
253  fentry->pz = pz;
254  fentry->E = enu;
255 
256  fnumi->run = 1;//gnumi->PassThroughInfo().run;
257  fnumi->evtno = nok;//gnumi->PassThroughInfo().evtno;
258  fnumi->entryno = nok;//gnumi->GetEntryNumber();
259 
260  fnumi->tpx = mjd;//gnumi->PassThroughInfo().tpx;
261  fnumi->tpy = 0;//gnumi->PassThroughInfo().tpy;
262  fnumi->tpz = 0;//gnumi->PassThroughInfo().tpz;
263  fnumi->vx = 0;//gnumi->PassThroughInfo().vx;
264  fnumi->vy = 0;//gnumi->PassThroughInfo().vy;
265  fnumi->vz = 0;//gnumi->PassThroughInfo().vz;
266 
267  fnumi->pdpx = 0;//gnumi->PassThroughInfo().pdpx;
268  fnumi->pdpy = 0;//gnumi->PassThroughInfo().pdpy;
269  fnumi->pdpz = 0;//gnumi->PassThroughInfo().pdpz;
270 
271  //double apppz = 0;//gnumi->PassThroughInfo().pppz;
272  fnumi->pppx = 0;//gnumi->PassThroughInfo().ppdxdz * apppz;
273  fnumi->pppy = 0;//gnumi->PassThroughInfo().ppdydz * apppz;
274  fnumi->pppz = 0;//apppz;
275 
276  fnumi->ndecay = 0;//gnumi->PassThroughInfo().ndecay;
277  fnumi->ptype = 0;//gnumi->PassThroughInfo().ptype;
278  fnumi->ppmedium = 0;//gnumi->PassThroughInfo().ppmedium;
279  fnumi->tptype = 0;//gnumi->PassThroughInfo().tptype;
280 
281  faux->auxdbl.push_back(wgt);
282  faux->auxint.push_back(nann);
283 
284  ++nok;
285  fluxntp->Fill();
286 
287  pdglist.insert(fentry->pdg);
288  minwgt = TMath::Min(minwgt, fentry->wgt);
289  maxwgt = TMath::Max(maxwgt, fentry->wgt);
290  maxe = TMath::Max(maxe, fentry->E);
291  } // if (not absorbed)
292  } // loop on nu's in batch
293  } // loop on annihilation batches
294  std::cout << "=========================== Complete " << std::endl;
295 
296  fmeta->pdglist.clear();
297  std::set<int>::const_iterator setitr = pdglist.begin();
298  for ( ; setitr != pdglist.end(); ++setitr) fmeta->pdglist.push_back(*setitr);
299 
300  fmeta->maxEnergy = maxe;
301  fmeta->minWgt = minwgt;
302  fmeta->maxWgt = maxwgt;
303  fmeta->protons = nann; // don't use nann_sofar because neutrino weighting is
304  // based on total requested annihilations
305  fmeta->metakey = metakey;
306  metantp->Fill();
307 
308  sw.Stop();
309  std::cout << "Generated " << nwrite << " (" << nok << " w/ E >"
310  << enumin << " ) " << " ( request " << nFluxEntries << " ) "
311  << std::endl
312  << nann_sofar << " nAnn " << " ( request " << nann << " )"
313  << std::endl
314  << "Time to generate: " << std::endl;
315  sw.Print();
316 
317  std::cout << "=================================================" << std::endl;
318 
319  std::cout << "Last GSimpleNtpEntry: " << std::endl
320  << *fentry
321  << std::endl
322  << "Last GSimpleNtpMeta: " << std::endl
323  << *fmeta
324  << std::endl;
325 
326  file->cd();
327  // write ntuples out
328  fluxntp->Write();
329  metantp->Write();
330  file->Close();
331 
332  return 0;
333 }
334 
335 
337  dsinit_();
338  wimpinit_();
339  nusetup_();
340  std::cout << "wimpann setup complete." << std::endl;
341 }
342 
343 
344 void get_wimpann_info(unsigned int& ann_chan,
345  double& cm_energy,
346  unsigned int& nann,
347  bool& wimpann_outp,
348  bool& use_earth_ann,
349  std::string& file_s,
350  std::string& file_e,
351  std::string& file_info,
352  bool& wimpann_sum,
353  std::string& file_sum_s,
354  std::string& file_sum_e) {
355  // command line input not yet implemented
356  //if (COMMAND_LINE_INPUT) {
357 
358  //}
359  //else {
360  ann_chan = ANN_CHAN;
361  cm_energy = CM_ENERGY; // == 10 GeV WIMPs
362  nann = N_ANN;
363  wimpann_outp = WIMPANN_OUTP;
364  use_earth_ann = USE_EARTH_ANN;
365  file_s = FILE_S;
366  file_e = FILE_E;
367  file_info = FILE_INFO;
368  wimpann_sum = WIMPANN_SUM;
369  file_sum_s = FILE_SUM_S;
370  file_sum_e = FILE_SUM_E;
371 
372  nuosc_.th12 = THETA_12;
373  nuosc_.th13 = THETA_13;
374  nuosc_.th23 = THETA_23;
375  nuosc_.delta = DELTA_CP;
376  nuosc_.dm212 = DM2_21;
377  nuosc_.dm312 = DM2_31;
378 
379  wimppara_.seed = (unsigned int) std::time(NULL);
380  wimppara_.nevfiles = 1;
381  wimppara_.nannperfile = nann;
382  wimppara_.fileno = 1;
383 
384  std::cout << "Using defaults for all wimpann parameters." << std::endl;
385  //}
386 }
387 
389 // if (COMMAND_LINE_INPUT) {
390 
391 // }
392 // else {
393  //wesim_.wafile = "null.txt"
394  wesim_.welat = WELAT;
395  wesim_.welong = WELONG;
396 
397  wesim_.wetf = WETF;
398  //wesim_.d1txt = "D0V 2013\4"
399  //wesim_.d2txt = "D0V 2014\4"
400  char time_text[27] = "D0V 2013 ";
401  wesim_.mjd1 = astxt2mjd_((char *) time_text);
402  time_text[7] = '4';
403  wesim_.mjd2 = astxt2mjd_((char *) time_text);
404 
405  wesim_.weint = WE_INT;
406  wesim_.wemed = WE_MED;
407 
408  //wesim_.wefile = "none";
409  //wesim_.wesumdir = "none";
410  wesim_.weseed = wimppara_.seed;
411 
412  std::cout << "Using defaults for all wimpevent parameters." << std::endl;
413 // }
414 }
const XML_Char int len
Definition: expat.h:262
void wimpinit_()
Double_t E
energy in lab frame
void RandGen(long int seed)
Definition: AppInit.cxx:31
const unsigned int N_ANN
Definition: runWimpSim.h:44
const unsigned int WE_INT
Definition: runWimpSim.h:65
const double THETA_13
Definition: runWimpSim.h:47
int main(int argc, char *argv[])
Definition: runWimpSim.cpp:12
void setup_wimpann()
Definition: runWimpSim.cpp:336
const double X_MIN
Definition: runWimpSim.h:38
void dsinit_()
const std::string FILE_SUM_E
Definition: runWimpSim.h:59
struct @20 nus_
Double_t pdpx
nu parent momentum at time of decay
Double_t px
x momentum in lab frame
Double_t vtxy
y position in lab frame
Float_t x1[n_points_granero]
Definition: compare.C:5
Double_t vx
vertex position of hadron/muon decay
Double_t maxEnergy
maximum energy
const double WELAT
Definition: runWimpSim.h:62
const unsigned int WE_MED
Definition: runWimpSim.h:66
double mwimp
Definition: runWimpSim.h:84
Float_t ss
Definition: plot.C:24
const unsigned int ANN_CHAN
Definition: runWimpSim.h:43
const unsigned int WETF
Definition: runWimpSim.h:64
const double WELONG
Definition: runWimpSim.h:63
const double DM2_31
Definition: runWimpSim.h:51
const std::string FILE_E
Definition: runWimpSim.h:56
const bool WIMPANN_OUTP
Definition: runWimpSim.h:52
Int_t tptype
parent particle type at target exit
std::vector< Double_t > auxdbl
additional doubles associated w/ entry
const double Y_MAX
Definition: runWimpSim.h:39
int nut
Definition: runWimpSim.h:118
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
std::vector< Int_t > auxint
additional ints associated w/ entry
double astxt2mjd_(char *)
const bool USE_EARTH_ANN
Definition: runWimpSim.h:53
const double Y_MIN
Definition: runWimpSim.h:39
const double NOVA_AZIMUTH
Definition: runWimpSim.h:30
void get_wimpann_info(unsigned int &, double &, unsigned int &, bool &, bool &, std::string &, std::string &, std::string &, bool &, std::string &, std::string &)
Definition: runWimpSim.cpp:344
const double DM2_21
Definition: runWimpSim.h:50
Double_t protons
represented number of protons-on-target
struct @24 wesim_
const double THETA_23
Definition: runWimpSim.h:48
void nusetup_()
void weinit_()
unsigned int nutype(unsigned int u)
Definition: runWimpSim.h:122
Double_t pppx
nu parent momentum at production point
double enu
Definition: runWimpSim.h:113
std::vector< Int_t > pdglist
list of neutrino flavors
void pygive_(char *, size_t)
const double X_MAX
Definition: runWimpSim.h:38
const double Z_MIN
Definition: runWimpSim.h:40
std::string getenv(std::string const &name)
Double_t vtxz
z position in lab frame
const double PI
Definition: runWimpSim.h:25
const unsigned nFluxEntries
Definition: runWimpSim.h:34
const double THETA_12
Definition: runWimpSim.h:46
void genevt_(unsigned int *channel, unsigned int *nev)
Double_t minWgt
minimum weight
Int_t ppmedium
tracking medium where parent was produced
const double j
Definition: BetheBloch.cxx:29
struct @19 neutrino_
const double DELTA_CP
Definition: runWimpSim.h:49
double zamp(complex z)
Definition: runWimpSim.h:130
Double_t tpx
parent particle px at target exit
struct @22 nuosc_
const ana::Var wgt
z
Definition: test.py:28
struct @23 wimppara_
const std::string FILE_SUM_S
Definition: runWimpSim.h:58
void nuoscinit_()
OStream cout
Definition: OStream.cxx:6
double zarg(complex z)
Definition: runWimpSim.h:131
const double enumin
Definition: runWimpSim.h:36
Double_t maxWgt
maximum weight
Float_t sw
Definition: plot.C:20
std::size_t Hash(const std::vector< std::size_t > &vals)
Generate a unique hash from a sequence of integers.
Definition: Utilities.cxx:384
T sin(T number)
Definition: d0nt_math.hpp:132
UInt_t metakey
index key to tie to individual entries
struct @25 weevent_
int get_pdg(unsigned wimpevent_nui)
Definition: runWimpSim.h:135
Double_t vtxx
x position in lab frame
Double_t pz
z momentum in lab frame
TFile * file
Definition: cellShifts.C:17
T cos(T number)
Definition: d0nt_math.hpp:78
TRandom3 r(0)
const double CM_ENERGY
Definition: runWimpSim.h:45
const std::string FILE_S
Definition: runWimpSim.h:55
Double_t dist
distance from hadron decay
const bool WIMPANN_SUM
Definition: runWimpSim.h:54
void get_wimpevent_info()
Definition: runWimpSim.cpp:388
void modpyt_()
const double Z_MAX
Definition: runWimpSim.h:40
UInt_t metakey
key to meta data
Double_t py
y momentum in lab frame
const std::string FILE_INFO
Definition: runWimpSim.h:57
void nupropsu_()
double mjd
Definition: runWimpSim.h:113
Int_t ptype
parent type (PDG)
enum BeamMode string