00001 #include <iostream.h>
00002 #include <math.h>
00003 #include <string.h>
00004 #include "ReactorEvent.hh"
00005 #include "ReactorDetector.hh"
00006 #include "ReactorFortran.hh"
00007 #include "ReactorNtuple.hh"
00008
00009
00010 ReactorNtuple::ReactorNtuple(const ReactorNtuple& Ntuple){
00011
00012 Id = Ntuple.Id;
00013 Nvar = Ntuple.Nvar;
00014 Offset = Ntuple.Offset;
00015 }
00016 ReactorNtuple& ReactorNtuple::operator=(const ReactorNtuple& rhs){
00017
00018 if(this == &rhs)return *this;
00019 Id = rhs.Id;
00020 Nvar = rhs.Nvar;
00021 Offset = rhs.Offset;
00022 return *this;
00023 }
00024
00025 ReactorNtuple::ReactorNtuple(int ntId,
00026 const char* ntName,
00027 char* dirName){
00028
00029
00030
00031
00032
00033
00034 int wordsize = 8;
00035 int ntVarsMax = 100;
00036 int ntVars=0;
00037 char ntTags[ntVarsMax*wordsize];
00038 strcpy(ntTags+ntVars++*wordsize,"wtc ");
00039 strcpy(ntTags+ntVars++*wordsize,"wtf ");
00040 strcpy(ntTags+ntVars++*wordsize,"enu ");
00041 strcpy(ntTags+ntVars++*wordsize,"KEpos ");
00042 strcpy(ntTags+ntVars++*wordsize,"poscos ");
00043 strcpy(ntTags+ntVars++*wordsize,"KEneu ");
00044 strcpy(ntTags+ntVars++*wordsize,"neucos ");
00045 strcpy(ntTags+ntVars++*wordsize,"phposvis ");
00046 strcpy(ntTags+ntVars++*wordsize,"phneuvis ");
00047 strcpy(ntTags+ntVars++*wordsize,"phposrec ");
00048 strcpy(ntTags+ntVars++*wordsize,"phneurec ");
00049 strcpy(ntTags+ntVars++*wordsize,"vnur ");
00050 strcpy(ntTags+ntVars++*wordsize,"vnucos ");
00051 strcpy(ntTags+ntVars++*wordsize,"vnuphi ");
00052 strcpy(ntTags+ntVars++*wordsize,"vnux ");
00053 strcpy(ntTags+ntVars++*wordsize,"vnuy ");
00054 strcpy(ntTags+ntVars++*wordsize,"vnuz ");
00055 strcpy(ntTags+ntVars++*wordsize,"vposr ");
00056 strcpy(ntTags+ntVars++*wordsize,"vposcos ");
00057 strcpy(ntTags+ntVars++*wordsize,"vposphi ");
00058 strcpy(ntTags+ntVars++*wordsize,"vposx ");
00059 strcpy(ntTags+ntVars++*wordsize,"vposy ");
00060 strcpy(ntTags+ntVars++*wordsize,"vposz ");
00061 strcpy(ntTags+ntVars++*wordsize,"vneur ");
00062 strcpy(ntTags+ntVars++*wordsize,"vneucos ");
00063 strcpy(ntTags+ntVars++*wordsize,"vneuphi ");
00064 strcpy(ntTags+ntVars++*wordsize,"vneux ");
00065 strcpy(ntTags+ntVars++*wordsize,"vneuy ");
00066 strcpy(ntTags+ntVars++*wordsize,"vneuz ");
00067 strcpy(ntTags+ntVars++*wordsize,"dnupos ");
00068 strcpy(ntTags+ntVars++*wordsize,"dnuneu ");
00069 strcpy(ntTags+ntVars++*wordsize,"tpos ");
00070 strcpy(ntTags+ntVars++*wordsize,"tneu ");
00071 strcpy(ntTags+ntVars++*wordsize,"npepos ");
00072 strcpy(ntTags+ntVars++*wordsize,"attenpos ");
00073 strcpy(ntTags+ntVars++*wordsize,"sfactpos ");
00074 strcpy(ntTags+ntVars++*wordsize,"nsteps ");
00075 strcpy(ntTags+ntVars++*wordsize,"absorb ");
00076 strcpy(ntTags+ntVars++*wordsize,"npeneu ");
00077 strcpy(ntTags+ntVars++*wordsize,"attenneu ");
00078 strcpy(ntTags+ntVars++*wordsize,"sfactneu ");
00079 strcpy(ntTags+ntVars++*wordsize,"ngampos ");
00080 strcpy(ntTags+ntVars++*wordsize,"ngamneu ");
00081 strcpy(ntTags+ntVars++*wordsize,"fpmtpos ");
00082 strcpy(ntTags+ntVars++*wordsize,"fpmtneu ");
00083 strcpy(ntTags+ntVars++*wordsize,"dposx ");
00084 strcpy(ntTags+ntVars++*wordsize,"dposy ");
00085 strcpy(ntTags+ntVars++*wordsize,"dposz ");
00086 strcpy(ntTags+ntVars++*wordsize,"dneux ");
00087 strcpy(ntTags+ntVars++*wordsize,"dneuy ");
00088 strcpy(ntTags+ntVars++*wordsize,"dneuz ");
00089
00090 hbookn_(ntId,ntName,ntVars,dirName,1000,ntTags,
00091 strlen(ntName),strlen(dirName),wordsize,ntVars);
00092
00093
00094
00095
00096 Nvar = ntVars;
00097 Id = ntId;
00098 Offset = 10000*(ntId-1);
00099
00100 hbook1_(Offset+101,"phpos",200,0,0.1,0,strlen("phpos"));
00101 hbook1_(Offset+102,"phneu",200,0,0.1,0,strlen("phneu"));
00102
00103 hbook1_(Offset+201,"tpos",250,0,25,0,strlen("tpos"));
00104 hbook1_(Offset+202,"tneu",250,0,500000,0,strlen("tneu"));
00105
00106 hbook1_(Offset+301,"ngpos",10,0,10,0,strlen("ngpos"));
00107 hbook1_(Offset+302,"ngneu",20,0,20,0,strlen("ngneu"));
00108
00109 hbook1_(Offset+401,"egpos",200,0,10,0,strlen("egpos"));
00110 hbook1_(Offset+402,"egneu",200,0,20,0,strlen("egneu"));
00111
00112 hbook1_(Offset+501,"rgp-rpos",200,0,100,0,strlen("rgp-rpos"));
00113 hbook1_(Offset+502,"rgn-rneu",200,0,100,0,strlen("rgn-rneu"));
00114
00115 hbook1_(Offset+601,"rgp-rnu",200,0,100,0,strlen("rgp-rnu"));
00116 hbook1_(Offset+602,"rgn-rnu",200,0,100,0,strlen("rgn-rnu"));
00117
00118 hbook1_(Offset+701,"tgp-tpos",250,0,5,0,strlen("tgp-tpos"));
00119 hbook1_(Offset+702,"tgn-tneu",250,0,5,0,strlen("tgn-tneu"));
00120
00121 hbook1_(Offset+801,"tgp-tnu",250,0,25,0,strlen("tgp-tnu"));
00122 hbook1_(Offset+802,"tgn-tnu",250,0,1000000,0,strlen("tgn-tnu"));
00123 }
00124
00125 ReactorNtuple::~ReactorNtuple(){
00126 ;
00127 }
00128
00129 void ReactorNtuple::Fill(ReactorEvent& Event,ReactorDetector& Detector){
00130
00131 float Tuple[Nvar];
00132 float* tuple = Tuple;
00133
00134 double rnu = Event.GetNubarVertexR();
00135 double cnu = Event.GetNubarVertexCosTheta();
00136 double fnu = Event.GetNubarVertexPhi();
00137 double xnu = rnu*cos(fnu)*sqrt(1-cnu*cnu);
00138 double ynu = rnu*sin(fnu)*sqrt(1-cnu*cnu);
00139 double znu = rnu*cnu;
00140
00141 double rpos = Event.GetPositronVertexR();
00142 double cpos = Event.GetPositronVertexCosTheta();
00143 double fpos = Event.GetPositronVertexPhi();
00144 double xpos = rpos*cos(fpos)*sqrt(1-cpos*cpos);
00145 double ypos = rpos*sin(fpos)*sqrt(1-cpos*cpos);
00146 double zpos = rpos*cpos;
00147
00148 double rneu = Event.GetNeutronVertexR();
00149 double cneu = Event.GetNeutronVertexCosTheta();
00150 double fneu = Event.GetNeutronVertexPhi();
00151 double xneu = rneu*cos(fneu)*sqrt(1-cneu*cneu);
00152 double yneu = rneu*sin(fneu)*sqrt(1-cneu*cneu);
00153 double zneu = rneu*cneu;
00154
00155 double dippos[3];
00156 double dipneu[3];
00157 for(int i=0;i<3;i++)dippos[i]=*(Detector.GetPositronDipole()+i);
00158 for(int i=0;i<3;i++)dipneu[i]=*(Detector.GetNeutronDipole()+i);
00159
00160 *tuple++ = Event.GetXCWeight();
00161 *tuple++ = Event.GetFluxWeight();
00162 *tuple++ = Event.GetEnu();
00163 *tuple++ = Event.GetPositronKE();
00164 *tuple++ = Event.GetPositronCosTheta();
00165 *tuple++ = Event.GetNeutronKE();
00166 *tuple++ = Event.GetNeutronCosTheta();
00167 *tuple++ = Detector.GetPHpositron();
00168 *tuple++ = Detector.GetPHneutron();
00169 *tuple++ = Detector.GetPHposReco();
00170 *tuple++ = Detector.GetPHneuReco();
00171 *tuple++ = rnu;
00172 *tuple++ = cnu;
00173 *tuple++ = fnu;
00174 *tuple++ = xnu;
00175 *tuple++ = ynu;
00176 *tuple++ = znu;
00177 *tuple++ = rpos;
00178 *tuple++ = cpos;
00179 *tuple++ = fpos;
00180 *tuple++ = xpos;
00181 *tuple++ = ypos;
00182 *tuple++ = zpos;
00183 *tuple++ = rneu;
00184 *tuple++ = cneu;
00185 *tuple++ = fneu;
00186 *tuple++ = xneu;
00187 *tuple++ = yneu;
00188 *tuple++ = zneu;
00189 *tuple++ = sqrt((xnu-xpos)*(xnu-xpos)+
00190 (ynu-ypos)*(ynu-ypos)+
00191 (znu-zpos)*(znu-zpos));
00192 *tuple++ = sqrt((xnu-xneu)*(xnu-xneu)+
00193 (ynu-yneu)*(ynu-yneu)+
00194 (znu-zneu)*(znu-zneu));
00195 *tuple++ = Event.GetPositronTime();
00196 *tuple++ = Event.GetNeutronTime();
00197 *tuple++ = Detector.GetMeanNpePos();
00198 *tuple++ = Detector.GetMeanAttenPos();
00199 *tuple++ = Detector.GetMeanSfactorPos();
00200 *tuple++ = Detector.GetNeutronSteps();
00201 *tuple++ = Detector.GetAbsorber();
00202 *tuple++ = Detector.GetMeanNpeNeu();
00203 *tuple++ = Detector.GetMeanAttenNeu();
00204 *tuple++ = Detector.GetMeanSfactorNeu();
00205 *tuple++ = Detector.GetNgamma('P');
00206 *tuple++ = Detector.GetNgamma('N');
00207 *tuple++ = Detector.GetHitFracPos();
00208 *tuple++ = Detector.GetHitFracNeu();
00209 *tuple++ = dippos[0];
00210 *tuple++ = dippos[1];
00211 *tuple++ = dippos[2];
00212 *tuple++ = dipneu[0];
00213 *tuple++ = dipneu[1];
00214 *tuple++ = dipneu[2];
00215
00216 hfn_(Id,Tuple);
00217
00218
00219
00220 int npmt=Detector.GetNpmt();
00221 for(int n=0;n<npmt;n++){
00222 int npos;
00223 int nneu;
00224 double xyzpmt[3];
00225 double phpos[100];
00226 double phneu[100];
00227 double tpos[100];
00228 double tneu[100];
00229
00230
00231
00232 Detector.GetPMTdata(n,xyzpmt,npos,phpos,tpos,nneu,phneu,tneu);
00233
00234 for(int nn=0;nn<npos;nn++){
00235 hfill_(Offset+101,*(phpos+nn));
00236 hfill_(Offset+201,*(tpos+nn));
00237 }
00238 for(int nn=0;nn<nneu;nn++){
00239 hfill_(Offset+102,*(phneu+nn));
00240 hfill_(Offset+202,*(tneu+nn));
00241 }
00242 }
00243
00244
00245
00246
00247
00248 int nposgamma = Detector.GetNgamma('P');
00249 hfill_(Offset+301,1.0*nposgamma);
00250 double egp[nposgamma];
00251 double tgp[nposgamma];
00252 double rgp[nposgamma];
00253 double cgp[nposgamma];
00254 double fgp[nposgamma];
00255 Detector.GetGamma('P',egp,tgp,rgp,cgp,fgp);
00256 for (int n=0;n<nposgamma;n++){
00257 hfill_(Offset+401,egp[n]);
00258 hfill_(Offset+701,tgp[n]-Event.GetPositronTime());
00259 hfill_(Offset+801,tgp[n]-Event.GetNubarTime());
00260 }
00261 double xyzp[3*nposgamma];
00262 Detector.GetGamma('P',xyzp);
00263 double* pxyz=xyzp;
00264 for (int n=0;n<nposgamma;n++){
00265 double* qxyz=Event.GetPositronVertexXYZ();
00266 double* rxyz=Event.GetNubarVertexXYZ();
00267 double dr1[3];
00268 double dr2[3];
00269 for(int i=0;i<3;i++){
00270 dr1[i] = *pxyz - *qxyz++;
00271 dr2[i] = *pxyz++ - *rxyz++;
00272 }
00273 hfill_(Offset+501,sqrt(dr1[0]*dr1[0]+dr1[1]*dr1[1]+dr1[2]*dr1[2]));
00274 hfill_(Offset+601,sqrt(dr2[0]*dr2[0]+dr2[1]*dr2[1]+dr2[2]*dr2[2]));
00275 }
00276
00277
00278
00279 int nneugamma = Detector.GetNgamma('N');
00280 hfill_(Offset+302,1.0*nneugamma);
00281 double egn[nneugamma];
00282 double tgn[nneugamma];
00283 double rgn[nneugamma];
00284 double cgn[nneugamma];
00285 double fgn[nneugamma];
00286 Detector.GetGamma('N',egn,tgn,rgn,cgn,fgn);
00287 for (int n=0;n<nneugamma;n++){
00288 hfill_(Offset+402,egn[n]);
00289 hfill_(Offset+702,tgn[n]-Event.GetNeutronTime());
00290 hfill_(Offset+802,tgn[n]-Event.GetNubarTime());
00291 }
00292 double xyzn[3*nneugamma];
00293 Detector.GetGamma('N',xyzn);
00294 pxyz=xyzn;
00295 for (int n=0;n<nneugamma;n++){
00296 double* qxyz=Event.GetNeutronVertexXYZ();
00297 double* rxyz=Event.GetNubarVertexXYZ();
00298 double dr1[3];
00299 double dr2[3];
00300 for(int i=0;i<3;i++){
00301 dr1[i] = *pxyz - *qxyz++;
00302 dr2[i] = *pxyz++ - *rxyz++;
00303 }
00304 hfill_(Offset+502,sqrt(dr1[0]*dr1[0]+dr1[1]*dr1[1]+dr1[2]*dr1[2]));
00305 hfill_(Offset+602,sqrt(dr2[0]*dr2[0]+dr2[1]*dr2[1]+dr2[2]*dr2[2]));
00306 }
00307 }
00308
