Go to the source code of this file.
Functions | |
| void | reweight_neutrino (double det_pos[], double vertex[], int neutrino_type, double neutrino_momentum[], int parent_particle_type, double parent_momentum[], double muparent_momentum[], double &eneu, double &weight) |
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 5 of file reweight.cxx. Referenced by GnumiTree::paint_flux(). 00014 {
00015 eneu = weight = 0.0;
00016
00017 const double pimass=0.13957, kmass=0.49368, k0mass=0.49767, mumass=0.105658389;
00018 const int nue=53, nuebar=52, numu=56, numubar=55, muplus=5, muminus=6;
00019
00020 double parent_mass;
00021 switch (parent_particle_type) {
00022 case 8: case 9: parent_mass = pimass; break;
00023 case 11: case 12: parent_mass = kmass; break;
00024 case 10: parent_mass = k0mass; break;
00025 case 5: case 6: parent_mass = mumass; break;
00026 default:
00027 cerr << "Unknown particle: " << parent_particle_type << endl;
00028 return;
00029 }
00030
00031 double parent_energy =0;
00032 for(int i=0; i<3; ++i)
00033 parent_energy += parent_momentum[i]*parent_momentum[i];
00034 parent_energy += parent_mass*parent_mass;
00035
00036 parent_energy = sqrt(parent_energy);
00037 double gamma = parent_energy / parent_mass;
00038
00039 double beta[3];
00040 for (int i=0; i<3; ++i) beta[i] = parent_momentum[i]/parent_energy;
00041
00042 double neutrino_energy = 0;
00043 double partial = 0;
00044 for (int i=0; i<3; ++i) {
00045 neutrino_energy += neutrino_momentum[i]*neutrino_momentum[i];
00046 partial += neutrino_momentum[i]*beta[i];
00047 }
00048 neutrino_energy = sqrt(neutrino_energy);
00049 partial *= gamma;
00050 double enuzr = gamma * neutrino_energy - partial;
00051
00052 //C get angle from parent line of flight to detector in lab frame
00053
00054 double rad = 0;
00055 double parentp = 0;
00056 for (int i=0; i<3; ++i) {
00057 double tmp = det_pos[i]-vertex[i];
00058 rad += tmp*tmp;
00059 parentp += parent_momentum[i]*parent_momentum[i];
00060 }
00061 rad = sqrt(rad);
00062 parentp = sqrt(parentp);
00063 double costh_pardet = 0;
00064 for (int i=0; i<3; ++i) costh_pardet += parent_momentum[i]*(det_pos[i]-vertex[i]);
00065 costh_pardet /= parentp * rad;
00066 if (costh_pardet > +1.0) costh_pardet = +1.0;
00067 if (costh_pardet < -1.0) costh_pardet = -1.0;
00068 double theta_pardet = acos(costh_pardet);
00069
00070 //C weighted neutrino energy in lab, approx, good for small theta
00071 double emrat= (2.0*gamma) / (1.0 + gamma*gamma*theta_pardet*theta_pardet);
00072 eneu=emrat*enuzr;
00073
00074 //C Get solid angle/4pi for detector element
00075 const double rdet = 100.0;
00076 double sangdet = (rdet*rdet/((det_pos[2]-vertex[2])*(det_pos[2]-vertex[2])))/4.0;
00077
00078 //C weight for solid angle and lorentz boost
00079 double wgt = sangdet * (emrat*emrat);
00080
00081 if (parent_particle_type==muplus || parent_particle_type == muminus) {
00082
00083 //C boost new neutrino to mu decay cm
00084 double Pnu[3];
00085 partial = 0;
00086 for (int i=0; i<3; ++i) {
00087 Pnu[i] = (det_pos[i]-vertex[i])*eneu/rad;
00088 partial += Pnu[i]*beta[i];
00089 }
00090 partial *= gamma;
00091 partial = eneu - partial/(gamma+1.0);
00092
00093 double P_dcm_nu[4];
00094 double P_dcm_nu2=0;
00095 for (int i=0; i<3; ++i) {
00096 P_dcm_nu[i] = Pnu[i] - beta[i]*gamma*partial;
00097 P_dcm_nu2 += P_dcm_nu[i]*P_dcm_nu[i];
00098 }
00099 P_dcm_nu[3] = sqrt(P_dcm_nu2);
00100
00101 //C boost parent of mu to mu production cm
00102 gamma = parent_energy/parent_mass;
00103 partial = 0;
00104 double muparent_energy = 0;
00105 for (int i=0; i<3; ++i) {
00106 beta[i] = parent_momentum[i]/parent_energy;
00107 partial += beta[i]*muparent_momentum[i];
00108 muparent_energy += muparent_momentum[i]*muparent_momentum[i];
00109 }
00110 muparent_energy = sqrt(muparent_energy + mumass*mumass);
00111 partial *= gamma;
00112 partial = muparent_energy - partial/(gamma+1.0);
00113
00114 double P_pcm_mp[3], P_pcm=0;
00115 for (int i=0; i<3; ++i) {
00116 P_pcm_mp[i] = muparent_momentum[i] - beta[i]*gamma*partial;
00117 P_pcm += P_pcm_mp[i]*P_pcm_mp[i];
00118 }
00119 P_pcm = sqrt(P_pcm);
00120
00121 //C calc new decay angle w.r.t. (anti)spin direction
00122 double costh = 0;
00123 for (int i=0; i<3; ++i) costh += P_dcm_nu[i]*P_pcm_mp[i];
00124 costh /= P_dcm_nu[3]*P_pcm;
00125 if (costh > +1) costh = +0.9999;
00126 if (costh < -1) costh = -0.9999;
00127
00128 //C calc relative weight due to angle difference
00129 double wt_ratio=0;
00130 if (neutrino_type == nue || neutrino_type == nuebar)
00131 wt_ratio = 1.0 - costh;
00132 else if(neutrino_type == numu || neutrino_type == numubar) {
00133 double xnu = 2. * enuzr / mumass;
00134 wt_ratio = ((3.0-2.0*xnu) - (1.0-2.0*xnu)*costh) / (3.0-2.0*xnu);
00135 }
00136 else {
00137 cerr << "Warning, unknown neutrino type: "
00138 << neutrino_type << endl;
00139 weight = 0;
00140 }
00141 wgt *= wt_ratio;
00142 }
00143
00144 weight = wgt;
00145 }
|
1.3.9.1