---

ReactorFsim.cpp

Go to the documentation of this file.

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

#define PAWC_SIZE 50000

// CERN recommended way of interfacing to HBOOK
typedef struct { float PAW[PAWC_SIZE]; } PAWC_DEF;
#define PAWC COMMON_BLOCK(PAWC,pawc)
COMMON_BLOCK_DEF(PAWC_DEF,PAWC);
PAWC_DEF PAWC;

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=290;
  double r2=350;
  int fast=0;
  int PMTdebug = 0;
  double MergeTime = 1;
  int trigger=0;
  char* fileName="test.ntuple";
  char* rootName="test.root";
  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;
  std::string trigsim;
  std::string PMTdebugoption;
  char pawFile[256];
  char rootFile[256];

  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);

      fin >> pawFile;
      fileName = pawFile;

      continue;
    }

    // 'A'ROOTFILE or 'a'rootfile specifies the output filename
    if(fin.peek() == 'a' || fin.peek() == 'A'){

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

      fin >> rootFile;
      rootName = rootFile;

      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;
    }

    // 'B'PMTOUTPUTFILE 'b'PMTOUTPUTFILE turns on/off the PMT Debug file
    if(fin.peek() == 'b' || fin.peek() == 'B'){

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

      fin >> PMTdebugoption;

      if(PMTdebugoption == "on"){
        PMTdebug = 1;
      }
      else{
        PMTdebug = 0;
      }
      continue;
    }

    // 'M'ERGETIME or 'm'ergetime specifies the PMTOutput Resolution (MergeTime)
    if(fin.peek() == 'm' || fin.peek() == 'M'){

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

      fin >> MergeTime;

      continue;
    }

    // 'T'RIGGER or 't'rigger turns on/off trigger simulation
    if(fin.peek() == 't' || fin.peek() == 'T'){

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

      fin >> trigsim;

      if(trigsim == "off"){
        trigger = 0;
      }
      else{
        trigger = 1;
      }

      continue;
    }

    // skip K, P, and D
    if(fin.peek() == 'K' || fin.peek() == 'k' ||
       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_5 by T. Bolton" << endl;
  cout << "===============================" << endl;
  cout << "  " << nevents << " events" << endl;
  cout << "  Paw output file: " << fileName << endl;
  cout << "  Root output file: " << rootName << endl;
  cout << "  R0 = " << r0 << ", R1 = " << r1 << ", R2 = " << r2
       << " (Rmaxgen = " << rmaxgen << ") cm" << endl;
  cout << "  Fast neutrons are " << fastneutron << endl;
  cout << "  Trigger simulation is " << trigsim << endl;
  cout << "  PMT output debugging is " << PMTdebugoption << endl;
  cout << "  PMT MergeTime is " << MergeTime << endl;
  cout << "===============================" << endl;
  //
  // init random number seed
  //
  int luxury = 2;
  int seed = 8796896;
  int k1=0;
  int k2=0;
  rluxgo_(luxury,seed,k1,k2);
  //
  // setup root output file and create tree
  //
  TFile *rfile= new TFile(rootName,"RECREATE","ReactorTree");
  ReactorTree rtree("RTree");
  //
  // Paw setup.
  //
  hlimit_(PAWC_SIZE);
  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_PMTDebugOption(PMTdebug);
  DetectorA.SetParam_PMTMergeTime(MergeTime);
  DetectorA.SetParam_R0(r0);
  DetectorA.SetParam_R1(r1);
  DetectorA.SetParam_R2(r2);
  ReactorNtuple ModelA(1,"model A",dirName);
  ReactorTrigger TriggerA;
  //
  // event loop
  //
  for (int events = 0;events<nevents;events++){
    //
    // create neutrino event
    //
    ReactorEvent Event(nevents,rmaxgen);
    //
    // create fast detector simulation
    //
    DetectorA.LightsOut(Event, events);
    ModelA.Fill(Event,DetectorA);
    rtree.Fill(Event,DetectorA);
    if(trigger) TriggerA.Trigger(Event,DetectorA);
  }
  // close up root
  //
     rfile->Write();
     rtree.CloseOut();
  // close up paw
  //
  hldir_(" "," ",strlen(" "),strlen(" "));
  char *dirName2 = new char[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));
}

---