00001
00010 #include "ReactorFortran.hh"
00011 #include "ReactorConstants.hh"
00012 #include "ReactorEvent.hh"
00013 #include "ReactorDetector.hh"
00014 #include "ReactorXC.hh"
00015 #include "ReactorNtuple.hh"
00016 #include "ReactorTrigger.hh"
00017 #include <iostream.h>
00018 #include <fstream.h>
00019 #include <iomanip.h>
00020 #include <math.h>
00021 #include <string.h>
00022 #include <cstdio>
00023
00030 int main(){
00031
00032
00033
00034
00035
00036
00037 int nevents = 100;
00038 double r0=260;
00039 double r1=260;
00040 double rmaxgen=270;
00041 double r2=350;
00042 int fast=0;
00043 int trigger=0;
00044 char* fileName="test.ntuple";
00045 char* ntName="event data";
00046
00047
00048 fstream fin("src/event_setup.txt", ios::in);
00049
00050
00051 char paramLine[256];
00052 std::string dummy;
00053 bool done = false;
00054 std::string fastneutron;
00055 std::string trigsim;
00056
00057 while(!done){
00058
00059
00060 if(fin.peek() == '#'){
00061
00062 fin.getline(paramLine, 256);
00063 continue;
00064 }
00065
00066
00067 if(fin.peek() == 'e' || fin.peek() == 'E'){
00068
00069 done = true;
00070 continue;
00071 }
00072
00073
00074 if(fin.peek() == 'n' || fin.peek() == 'N'){
00075
00076 fin >> dummy;
00077 while(fin.peek() == ' ' || fin.peek() == '=')
00078 fin.ignore(1);
00079
00080 fin >> nevents;
00081
00082 if(nevents <= 0){
00083 cout << " SETUP ERROR: number of events <= 0!" << endl;
00084 nevents = 1;
00085 }
00086
00087 continue;
00088 }
00089
00090
00091 if(fin.peek() == 'o' || fin.peek() == 'O'){
00092
00093 fin >> dummy;
00094 while(fin.peek() == ' ' || fin.peek() == '=')
00095 fin.ignore(1);
00096
00097 char newFile[256];
00098 fin >> newFile;
00099 fileName = newFile;
00100
00101 continue;
00102 }
00103
00104
00105 if(fin.peek() == 'r' || fin.peek() == 'R'){
00106
00107 fin >> dummy;
00108 while(fin.peek() == ' ' || fin.peek() == '=')
00109 fin.ignore(1);
00110
00111 fin >> r0 >> r1 >> r2 >> rmaxgen;
00112
00113 continue;
00114 }
00115
00116
00117 if(fin.peek() == 'f' || fin.peek() == 'F'){
00118
00119 fin >> dummy;
00120 while(fin.peek() == ' ' || fin.peek() == '=')
00121 fin.ignore(1);
00122
00123 fin >> fastneutron;
00124
00125 if(fastneutron == "off"){
00126 fast = 0;
00127 }
00128 else{
00129 fast = 1;
00130 }
00131
00132 continue;
00133 }
00134
00135
00136 if(fin.peek() == 't' || fin.peek() == 'T'){
00137
00138 fin >> dummy;
00139 while(fin.peek() == ' ' || fin.peek() == '=')
00140 fin.ignore(1);
00141
00142 fin >> trigsim;
00143
00144 if(trigsim == "off"){
00145 trigger = 0;
00146 }
00147 else{
00148 trigger = 1;
00149 }
00150
00151 continue;
00152 }
00153
00154
00155 if(fin.peek() == 'K' || fin.peek() == 'k' ||
00156 fin.peek() == 'P' || fin.peek() == 'p' ||
00157 fin.peek() == 'D' || fin.peek() == 'd'){
00158
00159 fin.getline(paramLine, 256);
00160 continue;
00161 }
00162
00163
00164 if(fin.peek() == ' ' || fin.peek() == '\n'){
00165
00166 fin.ignore(1);
00167 continue;
00168 }
00169 }
00170
00171 cout << "ReactorFsim v0_0_4 by T. Bolton" << endl;
00172 cout << "===============================" << endl;
00173 cout << " " << nevents << " events" << endl;
00174 cout << " Output file: " << fileName << endl;
00175 cout << " R0 = " << r0 << ", R1 = " << r1 << ", R2 = " << r2
00176 << " (Rmaxgen = " << rmaxgen << ") cm" << endl;
00177 cout << " Fast neutrons are " << fastneutron << endl;
00178 cout << " Trigger simulation is " << trigsim << endl;
00179 cout << "===============================" << endl;
00180
00181
00182
00183 int luxury = 2;
00184 int seed = 8796896;
00185 int k1=0;
00186 int k2=0;
00187 rluxgo_(luxury,seed,k1,k2);
00188
00189
00190
00191 hbook_init__();
00192
00193 int iostat;
00194 int ntLun=99;
00195 int ntBuf=1024;
00196 char* dirName = "TEST";
00197 char* chopt = "N";
00198 hropen_(ntLun,dirName,fileName,chopt,ntBuf,iostat,
00199 strlen(dirName),strlen(fileName),strlen(chopt));
00200
00201
00202
00203
00204
00205 ReactorDetector DetectorA;
00206 DetectorA.SetParam_FastNeutronOption(fast);
00207 DetectorA.SetParam_R0(r0);
00208 DetectorA.SetParam_R1(r1);
00209 DetectorA.SetParam_R2(r2);
00210 ReactorNtuple ModelA(1,"model A",dirName);
00211 ReactorTrigger TriggerA;
00212
00213
00214
00215 for (int events = 0;events<nevents;events++){
00216
00217
00218
00219 ReactorEvent Event(nevents,rmaxgen);
00220
00221
00222
00223 DetectorA.LightsOut(Event);
00224 ModelA.Fill(Event,DetectorA);
00225 if(trigger) TriggerA.Trigger(Event,DetectorA);
00226 }
00227
00228
00229
00230 hldir_(" "," ",strlen(" "),strlen(" "));
00231 char dirName2[strlen(dirName)+2];
00232 *dirName2='\0';
00233 strcat(dirName2,"//");
00234 strcat(dirName2,dirName);
00235 hcdir_(dirName2," ",strlen(dirName2),strlen(" "));
00236 int outopt = 0;
00237 int icycle;
00238 hrout_(outopt,icycle," ",strlen(" "));
00239 hrend_(dirName,strlen(dirName));
00240 }