reweight.h File Reference

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)

Function Documentation

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.

References NueConvention::nue, NueConvention::numu, and Munits::rad.

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 }


Generated on 2 Nov 2017 for loon by  doxygen 1.6.1