---

ReactorFsim.cpp

#include "ReactorFortran.hh"
#include "ReactorConstants.hh"
#include "ReactorEvent.hh"
#include "ReactorDetector.hh"
#include "ReactorXC.hh"
#include "ReactorNtuple.hh"
#include <iostream.h>
#include <math.h>
#include <string.h>

int main(){

  //  ReactorNeutrons Neutrons;
  //
  //  Set number of events, paw file name, and ntuple name
  //
  int nevents = 40000;
  double rmaxgen=225;
  int fast=0;
  char* fileName="th13.ntuple";
  char* ntName="event data";
  //
  // 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(175);
  DetectorA.SetParam_R1(200);
  DetectorA.SetParam_R2(275);
  ReactorNtuple ModelA(1,"model A",dirName);
  //
  ReactorDetector DetectorB;
  DetectorB.SetParam_FastNeutronOption(fast);
  DetectorB.SetParam_R0(175);
  DetectorB.SetParam_R1(200);
  DetectorB.SetParam_R2(275);
  DetectorB.SetParam_GdConcentration(0.11);
  ReactorNtuple ModelB(2,"10%higher f-Gd",dirName);
  //
  ReactorDetector DetectorC;
  DetectorC.SetParam_FastNeutronOption(fast);
  DetectorC.SetParam_R0(175);
  DetectorC.SetParam_R1(200);
  DetectorC.SetParam_R2(275);
  DetectorC.SetParam_attenlGd(440);
  ReactorNtuple ModelC(3,"10% higher l-Gd",dirName);
  //
  ReactorDetector DetectorD;
  DetectorD.SetParam_FastNeutronOption(fast);
  DetectorD.SetParam_R0(175);
  DetectorD.SetParam_R1(200);
  DetectorD.SetParam_R2(275);
  DetectorD.SetParam_attenlSc(1100);
  ReactorNtuple ModelD(4,"10% higher l-Sc",dirName);
  //
  ReactorDetector DetectorE;
  DetectorE.SetParam_FastNeutronOption(fast);
  DetectorE.SetParam_R0(175);
  DetectorE.SetParam_R1(200);
  DetectorE.SetParam_R2(275);
  ReactorNtuple ModelE(5,"-20% f-Gd vs t",dirName);
  //
  ReactorDetector DetectorF;
  DetectorF.SetParam_FastNeutronOption(fast);
  DetectorF.SetParam_R0(175);
  DetectorF.SetParam_R1(200);
  DetectorF.SetParam_R2(275);
  ReactorNtuple ModelF(6,"-20% l-Gd vs t",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);
    //
    DetectorB.LightsOut(Event);
    ModelB.Fill(Event,DetectorB);
    //
    DetectorC.LightsOut(Event);
    ModelC.Fill(Event,DetectorC);
    //
    DetectorD.LightsOut(Event);
    ModelD.Fill(Event,DetectorD);
    //
    DetectorE.SetParam_GdConcentration(0.01-0.002*events/nevents);
    DetectorE.LightsOut(Event);
    ModelE.Fill(Event,DetectorE);
    //
    DetectorF.SetParam_attenlGd(400-80*events/nevents);
    DetectorF.LightsOut(Event);
    ModelF.Fill(Event,DetectorF);
  }
  //
  // 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));
}

---