178 const float* flux = 0;
182 if (pdgcode== 12)
return flux[0];
183 if (pdgcode==-12)
return flux[1];
184 if (pdgcode== 14)
return flux[2];
185 if (pdgcode==-14)
return flux[3];
186 if (pdgcode== 16)
return flux[4];
187 if (pdgcode==-16)
return flux[5];
194 double numu,
double numubar,
195 double nutau,
double nutaubar)
205 double numu,
double numubar,
206 double nutau,
double nutaubar)
215 double numu,
double numubar,
216 double nutau,
double nutaubar)
234 const double pimass=.13957;
235 const double mumass=0.105658389;
236 const double kmass=0.49368;
237 const double k0mass=0.49767;
254 std::cerr<<
"dont know ptype "<<
fptype<<
"can not compute new decay weight"<<std::endl;
261 double gamma = Eplab/mass;
262 double beta = sqrt((gamma*gamma-1)/(gamma*gamma));
266 double rnx=1.*(x-
fvx);
267 double rny=1.*(y-
fvy);
268 double rnz=1.*(z-
fvz);
269 double rn=sqrt(rnx*rnx+rny*rny+rnz*rnz);
274 double costhetan = rndotp/(rn*p);
290 MN=1./(gamma*(1-beta*costhetan));
300 const double kRDET = 100.;
301 double san = (1.0-cos(atan( kRDET / rn )))/2.0;
310 double betav[3]={0.};
313 double P_nun[3]={0.};
314 double P_dcm_nun[4]={0.};
315 double P_pcm_mp[4]={0.};
317 betav[0] =
fpdpx/Eplab;
318 betav[1] =
fpdpy/Eplab;
319 betav[2] = fpdpz/Eplab;
321 P_nun[0] = rnx*newE/rn;
322 P_nun[1] = rny*newE/rn;
323 P_nun[2] = rnz*newE/rn;
325 partialn =gamma*(betav[0]*P_nun[0]+betav[1]*P_nun[1]+betav[2]*P_nun[2]);
326 partialn = newE - partialn /(gamma+1.);
328 P_dcm_nun[0] = P_nun[0] - betav[0]*gamma*partialn;
329 P_dcm_nun[1] = P_nun[1] - betav[1]*gamma*partialn;
330 P_dcm_nun[2] = P_nun[2] - betav[2]*gamma*partialn;
331 P_dcm_nun[3] = sqrt(pow(P_dcm_nun[0],2)
333 +pow(P_dcm_nun[2],2));
341 partial = gamma*(betav[0]*
fmuparpx +
343 partial =
fmupare-partial/(gamma+1.);
344 P_pcm_mp[0] =
fmuparpx - betav[0]*gamma*partial;
345 P_pcm_mp[1] =
fmuparpy - betav[1]*gamma*partial;
346 P_pcm_mp[2] =
fmuparpz - betav[2]*gamma*partial;
347 P_pcm_mp[3] = sqrt(pow(P_pcm_mp[0],2)+
353 if(P_dcm_nun[3]!=0&&P_pcm_mp[3]!=0){
354 costhn = ( P_dcm_nun[0]*P_pcm_mp[0]+
355 P_dcm_nun[1]*P_pcm_mp[1]+
356 P_dcm_nun[2]*P_pcm_mp[2])/(P_dcm_nun[3]*P_pcm_mp[3]);
370 double xnu = 2.*fnecm/mass;
371 wt_ration = ( (3.-2.*xnu) - (1.-2.*xnu)*costhn ) / (3.-2.*xnu);
387 <<
"MCFlux:" << std::endl
389 <<
" flux job " << std::setw(11) << mcflux.
frun <<
" " 390 <<
" pot # " << std::setw(11) << mcflux.
fevtno << std::endl
391 <<
" ntype " << std::setw(11) << mcflux.
fntype <<
" " 392 <<
" ptype " << std::setw(11) << mcflux.
fptype <<
" " 393 <<
" impwt " << std::setw(11) << mcflux.
fnimpwt << std::endl
394 <<
" ndecay " << std::setw(11) << mcflux.
fndecay <<
" " 395 <<
" ppmedium " << std::setw(11) << mcflux.
fppmedium <<
" " 396 <<
" tptype " << std::setw(11) << mcflux.
ftptype << std::endl
397 <<
" vxyz " << std::setw(11) << mcflux.
fvx <<
" " 398 << std::setw(11) << mcflux.
fvy <<
" " 399 << std::setw(11) << mcflux.
fvz << std::endl
400 <<
" pdpxyz " << std::setw(11) << mcflux.
fpdpx <<
" " 401 << std::setw(11) << mcflux.
fpdpy <<
" " 402 << std::setw(11) << mcflux.
fpdpz << std::endl
403 <<
" tpxyz " << std::setw(11) << mcflux.
ftpx <<
" " 404 << std::setw(11) << mcflux.
ftpy <<
" " 405 << std::setw(11) << mcflux.
ftpz << std::endl
406 <<
" ppd[xy]dz " << std::setw(11) << mcflux.
fppdxdz <<
" " 407 << std::setw(11) << mcflux.
fppdydz << std::endl
408 <<
" pppxyz " << std::setw(11) << mcflux.
fppdxdz*mcflux.
fpppz <<
" " 410 << std::setw(11) << mcflux.
fpppz << std::endl
411 <<
" muparpxyz " << std::setw(11) << mcflux.
fmuparpx <<
" " 412 << std::setw(11) << mcflux.
fmuparpy <<
" " 413 << std::setw(11) << mcflux.
fmuparpz <<
" " 414 <<
" mupare " << std::setw(11) << mcflux.
fmupare << std::endl
415 <<
" necm " << std::setw(11) << mcflux.
fnecm <<
" " 416 <<
" dk2gen " << std::setw(11) << mcflux.
fdk2gen << std::endl
417 <<
" near E " << std::setw(11) << mcflux.
fnenergyn <<
" " 418 <<
" wgt " << std::setw(11) << mcflux.
fnwtnear << std::endl
419 <<
" far E " << std::setw(11) << mcflux.
fnenergyf <<
" " 420 <<
" wgt " << std::setw(11) << mcflux.
fnwtfar << std::endl;
double Flux(int pdgcode, int which=0) const
double fgen2vtx
distance from ray origin to event vtx
Flux for negative horn focus.
double fgenx
origin of ray from flux generator
constexpr auto abs(T v)
Returns the absolute value of the argument.
float fFluxNeg[6]
Fluxes as aboce, for negative horn focus.
A bogus flux assumed by the generator.
double fdk2gen
distance from decay to ray origin
object containing MC flux information
Flux for positive horn focus.
void SetFluxNeg(double nue, double nuebar, double numu, double numubar, double nutau, double nutaubar)
void SetFluxGen(double nue, double nuebar, double numu, double numubar, double nutau, double nutaubar)
float fFluxGen[6]
Fluxes as above, assumed by generator.
void SetFluxPos(double nue, double nuebar, double numu, double numubar, double nutau, double nutaubar)
friend std::ostream & operator<<(std::ostream &output, const simb::MCFlux &mcflux)
float fFluxPos[6]
e,ebar,mu,mubar,tau,taubar flux, +horn focus
void ReDecay(double &newE, double &newW, double x, double y, double z)