---

ReactorFsim.cpp

Go to the documentation of this file.

#include "ReactorFortran.hh"
#include "ReactorConstants.hh"
#include "ReactorEvent.hh"
#include "ReactorDetector.hh"
#include "ReactorXC.hh"
#include "ReactorNtuple.hh"
#include <iostream.h>
#include <fstream.h>   // file I/O
#include <iomanip.h>   // format manipulation
#include <math.h>
#include <string.h>
#include <cstdio>

int main(){

  //  ReactorNeutrons Neutrons;
  //
  //  Set number of events, paw file name, and ntuple name
  //
  // defaults
  int nevents = 100;
  double r0=260;
  double r1=260;
  double rmaxgen=270;
  double r2=350;
  int fast=0;
  char* fileName="test.ntuple";
  char* ntName="event data";
  //
  // setup the event
  fstream fin("src/event_setup.txt", ios::in);
  //cout << "setup file: src/event_setup.txt" << endl;

  char paramLine[256];
  std::string dummy;
  bool done = false;
  std::string fastneutron;

  while(!done){

    // '#' specifies a comment line -- ignore it
    if(fin.peek() == '#'){

      fin.getline(paramLine, 256);
      continue;
    }

    // 'E'ND or 'e'nd specifies end of file -- set done flag to true
    if(fin.peek() == 'e' || fin.peek() == 'E'){

      done = true;
      continue;
    }

    // 'N'EVENT or 'n'event specifies the number of events
    if(fin.peek() == 'n' || fin.peek() == 'N'){

      fin >> dummy;
      while(fin.peek() == ' ' || fin.peek() == '=') 
        fin.ignore(1);

      fin >> nevents;

      if(nevents <= 0){
        cout << "  SETUP ERROR: number of events <= 0!" << endl;
        nevents = 1;
      }

      continue;
    }

    // 'O'UTPUTFILE or 'o'utputfile specifies the output filename
    if(fin.peek() == 'o' || fin.peek() == 'O'){

      fin >> dummy;
      while(fin.peek() == ' ' || fin.peek() == '=') 
        fin.ignore(1);

      char newFile[256];
      fin >> newFile;
      fileName = newFile;

      continue;
    }

    // 'R'ADII or 'r'adii specifies the three detector radii: R0 R1 R2
    if(fin.peek() == 'r' || fin.peek() == 'R'){

      fin >> dummy;
      while(fin.peek() == ' ' || fin.peek() == '=') 
        fin.ignore(1);

      fin >> r0 >> r1 >> r2 >> rmaxgen;

      continue;
    }

    // 'F'ASTNEUTRON or 'f'astneutron turns on/off fast neutrons
    if(fin.peek() == 'f' || fin.peek() == 'F'){

      fin >> dummy;
      while(fin.peek() == ' ' || fin.peek() == '=') 
        fin.ignore(1);

      fin >> fastneutron;

      if(fastneutron == "off"){
        fast = 0;
      }
      else{
        fast = 1;
      }

      continue;
    }

    // skip T, P. and D
    if(fin.peek() == 'T' || fin.peek() == 't' ||
       fin.peek() == 'P' || fin.peek() == 'p' ||
       fin.peek() == 'D' || fin.peek() == 'd'){

      fin.getline(paramLine, 256);
      continue;
    }

    // ignore extra whitespace
    if(fin.peek() == ' ' || fin.peek() == '\n'){

      fin.ignore(1);
      continue;
    }
  }
  //
  cout << "ReactorFsim v0_0_2 by T. Bolton" << endl;
  cout << "===============================" << endl;
  cout << "  " << nevents << " events" << endl;
  cout << "  Output file: " << fileName << endl;
  cout << "  R0 = " << r0 << ", R1 = " << r1 << ", R2 = " << r2 
       << " (Rmaxgen = " << rmaxgen << ") cm" << endl;
  cout << "  Fast neutrons are " << fastneutron << endl;
  cout << "===============================" << endl;
  //
  // init random number seed
  //
  int luxury = 2;
  int seed = 8796896;
  int k1=0;
  int k2=0;
  rluxgo_(luxury,seed,k1,k2);
  //
  // Paw setup.  
  //
  hbook_init__();
  //
  int iostat;
  int ntLun=99;
  int ntBuf=1024;
  char* dirName = "TEST";
  char* chopt = "N";
  hropen_(ntLun,dirName,fileName,chopt,ntBuf,iostat,
          strlen(dirName),strlen(fileName),strlen(chopt));
  //
  // Set up simple detector models
  // Set up simple Ntuple and histogranms 
  // See ReactorNtuple.hh for ntuple/histogram details.
  //
  ReactorDetector DetectorA;
  DetectorA.SetParam_FastNeutronOption(fast);
  DetectorA.SetParam_R0(r0);
  DetectorA.SetParam_R1(r1);
  DetectorA.SetParam_R2(r2);
  ReactorNtuple ModelA(1,"model A",dirName);
  //
  // event loop
  //
  for (int events = 0;events<nevents;events++){
    //
    // create neutrino event
    //
    ReactorEvent Event(rmaxgen);
    //
    // create fast detector simulation
    //
    DetectorA.LightsOut(Event);
    ModelA.Fill(Event,DetectorA);
  }
  //
  // close up paw
  //
  hldir_(" "," ",strlen(" "),strlen(" "));
  char dirName2[strlen(dirName)+2];
  *dirName2='\0';
  strcat(dirName2,"//");
  strcat(dirName2,dirName);
  hcdir_(dirName2," ",strlen(dirName2),strlen(" "));
  int outopt = 0;
  int icycle;
  hrout_(outopt,icycle," ",strlen(" "));
  hrend_(dirName,strlen(dirName));
}

---