runWimpSim.h

Go to the documentation of this file.

#ifndef RUN_WIMP_SIM_H
#define RUN_WIMP_SIM_H

#include <iostream>
#include <ctime>
#include <sstream>
#include <cmath>
#include <cstdlib>
#include <set>
#include <string>

#include "GENIE/Framework/Numerical/RandomGen.h"
#include "GENIE/Tools/Flux/GNuMIFlux.h"
#include "GENIE/Tools/Flux/GSimpleNtpFlux.h"
#include "GENIE/Framework/Utils/UnitUtils.h"
#include "GENIE/Framework/Utils/AppInit.h"
#include "TSystem.h"
#include "TStopwatch.h"
#include "TLorentzVector.h"
#include "TNtuple.h"
#include "TFile.h"
#include "TRandom3.h"

const bool COMMAND_LINE_INPUT = false;
const double PI = M_PI;

// clockwise angle from true North to NOvA +z axis
// taken from MoonShadowTrigger, which cites an email from Virgil Bocean
// to Alec Habig.
const double NOVA_AZIMUTH = 332.066111; // degrees

struct complex{ double r, i; };

const unsigned nFluxEntries = 100000;
//const unsigned nFluxEntries = 100;
const double enumin = 1.0; // GeV
const double T_MAX = 5.0e5, T_MIN = 0,
             X_MAX = 16000,   X_MIN = -16000,
             Y_MAX = -9200, Y_MIN = -9200,
             Z_MAX = 72000,  Z_MIN = -48000;

// wimpann defaults
const unsigned int ANN_CHAN      = 11, // 5 -> b bbar
                   N_ANN         = 1e7;
const double       CM_ENERGY     = 100, // 50 GeV WIMPs
                   THETA_12      = 33.58 * PI / 180.0,
                   THETA_13      = 9.12 * PI / 180.0,
                   THETA_23      = 40.40 * PI / 180.0,
                   DELTA_CP      = 0,
                   DM2_21        = 7.58e-5,
                   DM2_31        = 2.35e-5;
const bool         WIMPANN_OUTP  = false,
                   USE_EARTH_ANN = false,
                   WIMPANN_SUM   = false;
const std::string  FILE_S        = "wimpann_sun.out",
                   FILE_E        = "wimpann_earth.out",
                   FILE_INFO     = "wimpann.info",
                   FILE_SUM_S    = "wimpann_sun.sum",
                   FILE_SUM_E    = "wimpann_earth.sum";

// wimpevent defaults
const double       WELAT         = 48.4171297,
                   WELONG        = -92.7890442;
const unsigned int WETF          = 2,
                   WE_INT        = 3,
                   WE_MED        = 1;

extern"C" {
  void dsinit_();
  void wimpinit_();
  void nusetup_();
  void pygive_(char *, size_t);
  void nuoscinit_();
  void modpyt_();
  void genevt_(unsigned int *channel, unsigned int *nev);
  void weinit_();
  double astxt2mjd_(char *);
  void nupropsu_();
  double pyr_(int *);

  struct {
    double enumax;
    bool output, earth;
    double mwimp;
  } neutrino_;

  struct {
    complex nusc[700000][3];
    double nusr[700000][3], nuer[700000];
    unsigned int nusi[700000][7], nuei[700000][4], nnus, nnue, nuunierr;
  } nus_;

  extern struct {
    int kf1[14], kf2[14], ach, achmax;
  } annich_;
 
  struct {
    double th12, th13, th23, delta, dm212, dm312;
  } nuosc_;

  struct {
    unsigned int seed, nevfiles, nannperfile, fileno;
  } wimppara_;

  extern struct {
    double welat, welong, wet1, wet2, mjd1, mjd2;
    unsigned int wemed, weint, wetf, wesunacc, weradec, weseed;
    char wefile[128], wesumdir[128], d1txt[26], d2txt[26];
  } wesim_;

  struct {
    complex nustate[3];
    double enu, lang, le, ledet, langdet, hang, he, tev, dfluxnu,
           dfluxvert, dfluxdet, nuang, dist, nuaz, laz, haz, lazdet,
           sunnadir, mjd, mjdfr, sunra, sundec, sunaz, sunel, nura,
           nudec, nuel, lepra, lepdec, lepaz, lepel, hadra, haddec,
           hadaz, hadel, dpsun[3], lst;
    int nut, nui, iev, lcv;
  } weevent_;
}

inline unsigned int nutype(unsigned int u) { if (u==1||u==3||u==5) return 1;
                                             else return 2;
                                       }

inline int signum(double d) { if (d<0) return -1;
                              else if (d>0) return 1;
                              else return 0;
                            }
inline double zamp(complex z) { return sqrt(z.i*z.i+z.r*z.r); }
inline double zarg(complex z) { if (z.r==0 && z.i==0) return 0;
                                else if (z.r==0) return PI/2*signum(z.i);
                                else return atan(z.i/z.r);
                              }
inline int get_pdg (unsigned wimpevent_nui) {
  switch (wimpevent_nui) {
    case 1:
      return 11;// nu_e
    case 2:
      return -11;// ~nu_e
    case 3:
      return 14;// nu_mu
    case 4:
      return -14;//  ~nu_mu
    case 5:
      return 15;// nu_tau
    case 6:
      return -15;// ~nu_tau
  }
  return 0;
}

#endif