runWimpSim.cpp

Go to the documentation of this file.

#include "runWimpSim.h"


void setup_wimpann();
void get_wimpann_info(unsigned int&, double&, unsigned int&, 
                      bool&, bool&, std::string&,
                      std::string&, std::string&, bool&, std::string&,
                      std::string&);
void get_wimpevent_info();

int main(int argc, char* argv[]) {

  if (argc < 2) {
    std::cout << "usage: runWimpSim flux_out_filename" << std::endl;
    return 1;
  } 

  if (!getenv("WIMPSIM_DIR")) {
    std::cout << "wimpsim must be set up before running this program." << std::endl;
    std::cout << "do: setup wimpsim v3_05_00" << std::endl;
    return 1;
  }

  std::string fnameout(argv[1]);

//*********************** WIMPANN *************************
  setup_wimpann();
  // variables
  unsigned int ann_chan;
  double cm_energy, mwimp;
  unsigned int nann;
  bool wimpann_outp, use_earth_ann, wimpann_sum;
  std::string file_s, file_e, file_info, file_sum_s, file_sum_e;
  get_wimpann_info(ann_chan, cm_energy, nann, wimpann_outp, 
                   use_earth_ann, file_s, file_e, file_info,
                   wimpann_sum, file_sum_s, file_sum_e); 
  mwimp = cm_energy / 2.0;

  // setup pythia seed
  std::stringstream ss;
  ss << "MRPY(1)=" << wimppara_.seed;
  std::cout << "Using Pythia seed " << ss.str() << std::endl;
  //std::cout << ss.str().length() << std::endl;
  char cseed[20];
  for (size_t i=0; i<ss.str().length() && i<20; i++) {
     cseed[i] = ss.str().c_str()[i];
  }
  for (size_t i=ss.str().length(); i<20; i++) {
     cseed[i] = ' ';
  }
  //std::cout << cseed << std::endl;
  //cseed[20] = ' ';
  size_t len = 20;
  //pass hidden array length variable
  pygive_((char*) cseed, len);
 
  // Initialize neutrino oscillation routines
  nuoscinit_();

  // Modify Pythia
  modpyt_();

  //int nnustot=0, nnuetot=0;
  neutrino_.enumax = mwimp;
  neutrino_.mwimp = mwimp;

//*********************** WIMPEVENT *************************
  weinit_();
  std::cout << "wimpevent setup complete." << std::endl;
  get_wimpevent_info();

//*********************** GENIE *****************************

  //bool doaux = false;

  long int gseed = (long int) wimppara_.seed;
  genie::utils::app_init::RandGen(gseed);

  TFile* file = TFile::Open(fnameout.c_str(),"RECREATE");
  TTree* fluxntp = new TTree("flux","a simple flux n-tuple");
  TTree* metantp = new TTree("meta","metadata for flux n-tuple");
  genie::flux::GSimpleNtpEntry* fentry = new genie::flux::GSimpleNtpEntry;
  genie::flux::GSimpleNtpNuMI*  fnumi  = new genie::flux::GSimpleNtpNuMI;
  genie::flux::GSimpleNtpAux*   faux   = new genie::flux::GSimpleNtpAux;
  genie::flux::GSimpleNtpMeta*  fmeta  = new genie::flux::GSimpleNtpMeta;
  fluxntp->Branch("entry",&fentry);
  fluxntp->Branch("numi",&fnumi);
  fluxntp->Branch("aux",&faux);
  metantp->Branch("meta",&fmeta);

  //long int ngen = 0;

  std::set<int> pdglist;
  double maxe = 0;
  double minwgt = +1.0e10;
  double maxwgt = -1.0e10;

  UInt_t metakey = TString::Hash(fnameout.c_str(),strlen(fnameout.c_str()));
  std::cout << "metakey " << metakey << std::endl;
  // ensure that we don't get smashed by UInt_t vs Int_t
  metakey &= 0x7FFFFFFF;
  std::cout << "GENIE metakey " << metakey << " after 0x7FFFFFFF" << std::endl;

  TRandom3 *myRand = new TRandom3(0);

//********************** Simulation *************************
  std::cout << "=========================== Start " << std::endl;
  
  TStopwatch sw;
  sw.Start();

  long int nann_sofar = 0;

  unsigned nok = 0, nwrite = 0;
  // simulate in batches of 100,000
  for (unsigned int i=0; i<=nann/100000; i++) { // removed nok < nFluxEntries
    unsigned int todo = (nann-i*100000<100000) ? nann-i*100000 : 100000;
    genevt_(&ann_chan, &todo);
    nann_sofar += todo; 

    // now neutrino information is given in nus_
    // run the current set of neutrinos through wimpevent
    for (unsigned int j=0; j<nus_.nnus; j++) {
      if (!!nus_.nusi[j][6]&0x1) { // not absorbed

         nwrite++;

         double enu = nus_.nusr[j][1];
         //std::cout << "\n" << enu << std::endl;
         weevent_.enu = enu;
         if (enu < enumin) continue; // energy too low

         unsigned nut = nutype(nus_.nusi[j][4]);
         //std::cout << std::endl << nus_.nusi[j][4] << std::endl;
         //std::cout << nut << std::endl;
         weevent_.nut = nut;

         // see wimpann/writeevs
         double a_e, a_mu, a_tau, phi_e, phi_mu, phi_tau;
         a_e = zamp(nus_.nusc[j][0]);
         phi_e = zarg(nus_.nusc[j][0]);
         a_mu = zamp(nus_.nusc[j][1]);
         phi_mu = zarg(nus_.nusc[j][1]);
         a_tau = zamp(nus_.nusc[j][2]);
         phi_tau = zarg(nus_.nusc[j][2]);

         phi_mu = phi_mu-phi_e;
         phi_tau = phi_tau-phi_e;

         //std::cout << "a_e: " << a_e << std::endl;
         //std::cout << "a_mu: " << a_mu << std::endl;
         //std::cout << "phi_mu: " << phi_mu << std::endl;

         // now setup event info for wimpevent
         weevent_.nustate[0].r = a_e;
         weevent_.nustate[0].i = 0;
         weevent_.nustate[1].r = a_mu*cos(phi_mu);
         weevent_.nustate[1].i = a_mu*sin(phi_mu);
         weevent_.nustate[2].r = a_tau*cos(phi_tau);
         weevent_.nustate[2].i = a_tau*sin(phi_tau);

         //std::cout << "mu.r before: " << weevent_.nustate[1].r << std::endl;
         //std::cout << "a_mu before: " << zamp(weevent_.nustate[1]) << std::endl;

         // propagate neutrino to detector
         nupropsu_();

         //std::cout << "a_mu after: " << zamp(weevent_.nustate[1]) << std::endl;
         // project onto a neutrino state at the detector 
         // see wimpann/nuproject
         double pe  = weevent_.nustate[0].r*weevent_.nustate[0].r
                    + weevent_.nustate[0].i*weevent_.nustate[0].i,
                pmu = weevent_.nustate[1].r*weevent_.nustate[1].r
                    + weevent_.nustate[1].i*weevent_.nustate[1].i;
         //std::cout << "pe: " << pe << std::endl;
         //std::cout << "pmu: " << pmu << std::endl;
         int fam;
         //int temp = 0;
         double r = myRand->Rndm();
         //std::cout << "r: " << r << std::endl;
         if (r < pe) fam = 1;
         else if (r < (pe+pmu)) fam = 2;
         else fam = 3;
         weevent_.nui = fam*2-2+weevent_.nut;
         //std::cout << "nut: " << weevent_.nut << std::endl;
         //std::cout << "nui: " << weevent_.nui << std::endl;
         if (weevent_.nui != 3 && weevent_.nui != 4)
           continue; // only want nu_mu or nu_mu~
                  

         // save event into gnumiflux
         fentry->Reset();
         fnumi->Reset();
         faux->Reset();

         int pdg = get_pdg(weevent_.nui);
         double wgt, x, y, z, px, py, pz;//, t;
         wgt = 1.0/(4*PI*weevent_.dist*weevent_.dist)/nann;
         x   = (X_MAX - X_MIN)*myRand->Rndm()+X_MIN;
         y   = (Y_MAX - Y_MIN)*myRand->Rndm()+Y_MIN;
         z   = (Z_MAX - Z_MIN)*myRand->Rndm()+Z_MIN;
         //t   = (T_MAX - T_MIN)*myRand->Rndm()+T_MIN;
         // GENIE wants units in meters
         x /= 100;
         y /= 100;
         z /= 100;

         double part_az = weevent_.nuaz;
         double part_el = weevent_.nuel;
         double mjd     = weevent_.mjd;

//Original code
/*
         if (part_el*180/PI < 10) continue; // only want upward-going

         //coordinate transformations taken from MoonShadowTrigger
         double phi   = part_az-(NOVA_AZIMUTH*PI/180.0); // 332o03'58.07169"
         px = -1.0 * enu * cos(part_el) * sin(phi);
         py = enu * sin(part_el);
         pz = enu * cos(part_el) * cos(phi);
*/

//New code
         //coordinate transformations taken from MoonShadowTrigger
         double phi   = part_az-(NOVA_AZIMUTH*PI/180.0); // 332o03'58.07169"

         //part_el and phi are the angles where the particles came from!
         /*         px = -1.0 * enu * cos(-part_el) * sin(phi+180);
         py = enu * sin(-part_el);          
         pz = enu * cos(-part_el) * cos(phi+180);*/

        
         px =   enu * cos(part_el) * sin(phi);
         py = - enu * sin(part_el);          
         pz = - enu * cos(part_el) * cos(phi);
         
         //      if (py < 0) continue; //only up mu

         //x ... west
         //y ... up
         //z ... north

         fentry->metakey = metakey;
         fentry->pdg     = pdg;
         fentry->wgt     = wgt;
         fentry->vtxx    = x;
         fentry->vtxy    = y;
         fentry->vtxz    = z;
         fentry->dist    = 0;//gnumi->GetDecayDist();
         fentry->px      = px;
         fentry->py      = py;
         fentry->pz      = pz;
         fentry->E       = enu;

         fnumi->run      = 1;//gnumi->PassThroughInfo().run;
         fnumi->evtno    = nok;//gnumi->PassThroughInfo().evtno;
         fnumi->entryno  = nok;//gnumi->GetEntryNumber();

         fnumi->tpx      = mjd;//gnumi->PassThroughInfo().tpx;
         fnumi->tpy      = 0;//gnumi->PassThroughInfo().tpy;
         fnumi->tpz      = 0;//gnumi->PassThroughInfo().tpz;
         fnumi->vx       = 0;//gnumi->PassThroughInfo().vx;
         fnumi->vy       = 0;//gnumi->PassThroughInfo().vy;
         fnumi->vz       = 0;//gnumi->PassThroughInfo().vz;

         fnumi->pdpx     = 0;//gnumi->PassThroughInfo().pdpx;
         fnumi->pdpy     = 0;//gnumi->PassThroughInfo().pdpy;
         fnumi->pdpz     = 0;//gnumi->PassThroughInfo().pdpz;

         //double apppz = 0;//gnumi->PassThroughInfo().pppz;
         fnumi->pppx     = 0;//gnumi->PassThroughInfo().ppdxdz * apppz;
         fnumi->pppy     = 0;//gnumi->PassThroughInfo().ppdydz * apppz;
         fnumi->pppz     = 0;//apppz;

         fnumi->ndecay   = 0;//gnumi->PassThroughInfo().ndecay;
         fnumi->ptype    = 0;//gnumi->PassThroughInfo().ptype;
         fnumi->ppmedium = 0;//gnumi->PassThroughInfo().ppmedium;
         fnumi->tptype   = 0;//gnumi->PassThroughInfo().tptype;

         faux->auxdbl.push_back(wgt);
         faux->auxint.push_back(nann);

         ++nok;
         fluxntp->Fill();

         pdglist.insert(fentry->pdg);
         minwgt = TMath::Min(minwgt, fentry->wgt);
         maxwgt = TMath::Max(maxwgt, fentry->wgt);
         maxe   = TMath::Max(maxe, fentry->E);
      } // if (not absorbed)
    } // loop on nu's in batch
  } // loop on annihilation batches
  std::cout << "=========================== Complete " << std::endl;

  fmeta->pdglist.clear();
  std::set<int>::const_iterator setitr = pdglist.begin();
  for ( ; setitr != pdglist.end(); ++setitr)  fmeta->pdglist.push_back(*setitr);

  fmeta->maxEnergy = maxe;
  fmeta->minWgt    = minwgt;
  fmeta->maxWgt    = maxwgt;
  fmeta->protons   = nann; // don't use nann_sofar because neutrino weighting is
                           // based on total requested annihilations
  fmeta->metakey   = metakey;
  metantp->Fill();

  sw.Stop();
  std::cout << "Generated " << nwrite << " (" << nok << " w/ E >" 
       << enumin << " ) " << " ( request " << nFluxEntries << " ) "
       << std::endl
       << nann_sofar << " nAnn " << " ( request " << nann << " )"
       << std::endl
       << "Time to generate: " << std::endl;
  sw.Print();

  std::cout << "=================================================" << std::endl;

  std::cout << "Last GSimpleNtpEntry: " << std::endl
       << *fentry
       << std::endl
       << "Last GSimpleNtpMeta: " << std::endl
       << *fmeta
       << std::endl;

  file->cd();
  // write ntuples out
  fluxntp->Write();
  metantp->Write();
  file->Close();

  return 0;
}


void setup_wimpann() {
  dsinit_();
  wimpinit_();
  nusetup_();
  std::cout << "wimpann setup complete." << std::endl;
}


void get_wimpann_info(unsigned int&  ann_chan, 
                      double&       cm_energy, 
                      unsigned int&      nann,
                      bool&      wimpann_outp,        
                      bool&     use_earth_ann, 
                      std::string&     file_s,
                      std::string&     file_e, 
                      std::string&  file_info, 
                      bool&       wimpann_sum, 
                      std::string& file_sum_s,
                      std::string& file_sum_e) {
  // command line input not yet implemented
  //if (COMMAND_LINE_INPUT) {
     
  //}
  //else {
      ann_chan = ANN_CHAN;
      cm_energy = CM_ENERGY; // == 10 GeV WIMPs
      nann = N_ANN;
      wimpann_outp = WIMPANN_OUTP;
      use_earth_ann = USE_EARTH_ANN;
      file_s = FILE_S;
      file_e = FILE_E;
      file_info = FILE_INFO;
      wimpann_sum = WIMPANN_SUM;
      file_sum_s = FILE_SUM_S;
      file_sum_e = FILE_SUM_E;

      nuosc_.th12  = THETA_12;
      nuosc_.th13  = THETA_13;
      nuosc_.th23  = THETA_23;
      nuosc_.delta = DELTA_CP;
      nuosc_.dm212 = DM2_21;
      nuosc_.dm312 = DM2_31;

      wimppara_.seed = (unsigned int) std::time(NULL);
      wimppara_.nevfiles = 1;
      wimppara_.nannperfile = nann;
      wimppara_.fileno = 1;

      std::cout << "Using defaults for all wimpann parameters." << std::endl;
  //}
}

void get_wimpevent_info() {
//  if (COMMAND_LINE_INPUT) {

//  }
//  else {
    //wesim_.wafile = "null.txt"
    wesim_.welat = WELAT;
    wesim_.welong = WELONG;

    wesim_.wetf = WETF;
    //wesim_.d1txt = "D0V 2013\4"
    //wesim_.d2txt = "D0V 2014\4"
    char time_text[27] = "D0V 2013                  ";
    wesim_.mjd1 = astxt2mjd_((char *) time_text);
    time_text[7] = '4';
    wesim_.mjd2 = astxt2mjd_((char *) time_text);

    wesim_.weint = WE_INT;
    wesim_.wemed = WE_MED;

    //wesim_.wefile = "none";
    //wesim_.wesumdir = "none";
    wesim_.weseed = wimppara_.seed;

    std::cout << "Using defaults for all wimpevent parameters." << std::endl;
//  } 
}