calcLocationWeights.cxx
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 
4 #include "tree/calcLocationWeights.h"
5 
6 #include "tree/dkmeta.h"
7 #include "tree/dk2nu.h"
8 
9 /// user interface
11 {
12  size_t nloc = dkmeta->location.size();
13  for (size_t iloc = 0; iloc < nloc; ++iloc ) {
14  // skip calculation for random location ... should already be filled
15  const std::string rkey = "random decay";
16  if ( dkmeta->location[iloc].name == rkey ) {
17  if ( iloc != 0 ) {
18  std::cerr << "calcLocationWeights \"" << rkey << "\""
19  << " isn't the 0-th entry" << std::endl;
20  return; // assert(0);
21  }
22  if ( dk2nu->nuray.size() != 1 ) {
23  std::cerr << "calcLocationWeights \"" << rkey << "\""
24  << " nuenergy[" << iloc << "] not filled" << std::endl;
25  return; //assert(0);
26  }
27  continue;
28  }
29  TVector3 xyzDet(dkmeta->location[iloc].x,
30  dkmeta->location[iloc].y,
31  dkmeta->location[iloc].z); // position to evaluate
32  double enu_xy = 0; // give a default value
33  double wgt_xy = 0; // give a default value
34  int status = bsim::calcEnuWgt(dk2nu,xyzDet,enu_xy,wgt_xy);
35  if ( status != 0 ) {
36  std::cerr << "bsim::calcEnuWgt returned " << status << " for "
37  << dkmeta->location[iloc].name << std::endl;
38  }
39  // with the recalculated energy compute the momentum components
40  TVector3 xyzDk(dk2nu->decay.vx,dk2nu->decay.vy,dk2nu->decay.vz); // origin of decay
41  TVector3 p3 = enu_xy * (xyzDet - xyzDk).Unit();
42  bsim::NuRay anuray(p3.x(), p3.y(), p3.z(), enu_xy, wgt_xy);
43  dk2nu->nuray.push_back(anuray);
44  }
45 }
46 
47 //___________________________________________________________________________
48 int bsim::calcEnuWgt(const bsim::Decay& decay, const TVector3& xyz,
49  double& enu, double& wgt_xy)
50 {
51  // Neutrino Energy and Weight at arbitrary point
52  // Based on:
53  // NuMI-NOTE-BEAM-0109 (MINOS DocDB # 109)
54  // Title: Neutrino Beam Simulation using PAW with Weighted Monte Carlos
55  // Author: Rick Milburn
56  // Date: 1995-10-01
57 
58  // History:
59  // jzh 3/21/96 grab R.H.Milburn's weighing routine
60  // jzh 5/ 9/96 substantially modify the weighting function use dot product
61  // instead of rotation vecs to get theta get all info except
62  // det from ADAMO banks neutrino parent is in Particle.inc
63  // Add weighting factor for polarized muon decay
64  // jzh 4/17/97 convert more code to double precision because of problems
65  // with Enu>30 GeV
66  // rwh 10/ 9/08 transliterate function from f77 to C++
67 
68  // Original function description:
69  // Real function for use with PAW Ntuple To transform from destination
70  // detector geometry to the unit sphere moving with decaying hadron with
71  // velocity v, BETA=v/c, etc.. For (pseudo)scalar hadrons the decays will
72  // be isotropic in this sphere so the fractional area (out of 4-pi) is the
73  // fraction of decays that hit the target. For a given target point and
74  // area, and given x-y components of decay transverse location and slope,
75  // and given decay distance from target ans given decay GAMMA and
76  // rest-frame neutrino energy, the lab energy at the target and the
77  // fractional solid angle in the rest-frame are determined.
78  // For muon decays, correction for non-isotropic nature of decay is done.
79 
80  // Arguments:
81  // dk2nu :: contains current decay information
82  // xyz :: 3-vector of position to evaluate
83  // in *beam* frame coordinates (cm units)
84  // enu :: resulting energy
85  // wgt_xy :: resulting weight
86  // Return:
87  // (int) :: error code
88  // Assumptions:
89  // Energies given in GeV
90  // Particle codes have been translated from GEANT into PDG codes
91 
92  // for now ... these masses _should_ come from TDatabasePDG
93  // but use these hard-coded values to "exactly" reproduce old code
94  //
95  // old mass values are v01_07_* and before
96  // new mass values (v01_08_* and after) are Geant4 v4.10.3 values
97  //
98 #ifdef HISTORIC_MASS
99  const double kPIMASS = 0.13957;
100  const double kKMASS = 0.49368;
101  const double kK0MASS = 0.49767;
102  const double kMUMASS = 0.105658389;
103  const double kOMEGAMASS = 1.67245;
104 #else
105  const double kPIMASS = 0.1395701; // 0.13957;
106  const double kKMASS = 0.493677; // 0.49368;
107  const double kK0MASS = 0.497614; // 0.49767;
108  const double kMUMASS = 0.1056583715; // 0.105658389;
109  const double kOMEGAMASS = 1.67245; // 1.67245;
110 #endif
111  // from CLHEP/Units/PhysicalConstants.h
112  // used by Geant as CLHEP::neutron_mass_c2
113  const double kNEUTRONMASS = 0.93956536;
114 
115  const int kpdg_nue = 12; // extended Geant 53
116  const int kpdg_nuebar = -12; // extended Geant 52
117  const int kpdg_numu = 14; // extended Geant 56
118  const int kpdg_numubar = -14; // extended Geant 55
119 
120  const int kpdg_muplus = -13; // Geant 5
121  const int kpdg_muminus = 13; // Geant 6
122  const int kpdg_pionplus = 211; // Geant 8
123  const int kpdg_pionminus = -211; // Geant 9
124  const int kpdg_k0long = 130; // Geant 10 ( K0=311, K0S=310 )
125  const int kpdg_k0short = 310; // Geant 16
126  const int kpdg_k0mix = 311;
127  const int kpdg_kaonplus = 321; // Geant 11
128  const int kpdg_kaonminus = -321; // Geant 12
129  const int kpdg_omegaminus = 3334; // Geant 24
130  const int kpdg_omegaplus = -3334; // Geant 32
131  const int kpdg_neutron = 2112;
132  const int kpdg_antineutron = -2112;
133 
134  const double kRDET = 100.0; // set to flux per 100 cm radius
135 
136  double xpos = xyz.X();
137  double ypos = xyz.Y();
138  double zpos = xyz.Z();
139 
140  enu = 0.0; // don't know what the final value is
141  wgt_xy = 0.0; // but set these in case we return early due to error
142 
143 
144  // in principle we should get these from the particle DB
145  // but for consistency testing use the hardcoded values
146  double parent_mass = kPIMASS;
147  switch ( decay.ptype ) {
148  case kpdg_pionplus:
149  case kpdg_pionminus:
150  parent_mass = kPIMASS;
151  break;
152  case kpdg_kaonplus:
153  case kpdg_kaonminus:
154  parent_mass = kKMASS;
155  break;
156  case kpdg_k0long:
157  case kpdg_k0short:
158  case kpdg_k0mix:
159  parent_mass = kK0MASS;
160  break;
161  case kpdg_muplus:
162  case kpdg_muminus:
163  parent_mass = kMUMASS;
164  break;
165  case kpdg_omegaminus:
166  case kpdg_omegaplus:
167  parent_mass = kOMEGAMASS;
168  break;
169  case kpdg_neutron:
170  case kpdg_antineutron:
171  parent_mass = kNEUTRONMASS;
172  break;
173  default:
174  std::cerr << "bsim::calcEnuWgt unknown particle type " << decay.ptype
175  << std::endl << std::flush;
176  enu = 0.0;
177  wgt_xy = 0.0;
178  return 1;
179  }
180 
181  double parentp2 = ( decay.pdpx*decay.pdpx +
182  decay.pdpy*decay.pdpy +
183  decay.pdpz*decay.pdpz );
184  double parent_energy = TMath::Sqrt( parentp2 +
185  parent_mass*parent_mass);
186  double parentp = TMath::Sqrt( parentp2 );
187 
188  double gamma = parent_energy / parent_mass;
189  double gamma_sqr = gamma * gamma;
190  double beta_mag = TMath::Sqrt( ( gamma_sqr - 1.0 )/gamma_sqr );
191 
192  // Get the neutrino energy in the parent decay CM
193  double enuzr = decay.necm;
194  // Get angle from parent line of flight to chosen point in beam frame
195  double rad = TMath::Sqrt( (xpos-decay.vx)*(xpos-decay.vx) +
196  (ypos-decay.vy)*(ypos-decay.vy) +
197  (zpos-decay.vz)*(zpos-decay.vz) );
198 
199  double emrat = 1.0;
200  double costh_pardet = -999., theta_pardet = -999.;
201 
202  // boost correction, but only if parent hasn't stopped
203  if ( parentp > 0. ) {
204  costh_pardet = ( decay.pdpx*(xpos-decay.vx) +
205  decay.pdpy*(ypos-decay.vy) +
206  decay.pdpz*(zpos-decay.vz) )
207  / ( parentp * rad);
208  if ( costh_pardet > 1.0 ) costh_pardet = 1.0;
209  if ( costh_pardet < -1.0 ) costh_pardet = -1.0;
210  theta_pardet = TMath::ACos(costh_pardet);
211 
212  // Weighted neutrino energy in beam, approx, good for small theta
213  emrat = 1.0 / ( gamma * ( 1.0 - beta_mag * costh_pardet ));
214  }
215 
216  enu = emrat * enuzr; // the energy ... normally
217 
218  // Get solid angle/4pi for detector element
219  // small angle approximation, fixed by Alex Radovic
220  //SAA// double sangdet = ( kRDET*kRDET /
221  //SAA// ( (zpos-decay.vz)*(zpos-decay.vz) ) ) / 4.0;
222  double sanddetcomp = TMath::Sqrt( ( (xpos-decay.vx)*(xpos-decay.vx) ) +
223  ( (ypos-decay.vy)*(ypos-decay.vy) ) +
224  ( (zpos-decay.vz)*(zpos-decay.vz) ) );
225  double sangdet = (1.0-TMath::Cos(TMath::ATan( kRDET / sanddetcomp )))/2.0;
226 
227  // Weight for solid angle and lorentz boost
228  wgt_xy = sangdet * ( emrat * emrat ); // ! the weight ... normally
229 
230  // Done for all except polarized muon decay
231  // in which case need to modify weight
232  // (must be done in double precision)
233  // BUT do this only for case of muon decay, not muon capture
234  // until beamline simulation code gets updated these generally show up as
235  // decay.ndecay == 0, but certainly not dkp_mup_nusep or dkp_mum_nusep
236  // so was:
237  // if ( decay.ptype == kpdg_muplus ||
238  // decay.ptype == kpdg_muminus ) {
239  // now:
240  if ( decay.ndecay == bsim::dkp_mup_nusep ||
241  decay.ndecay == bsim::dkp_mum_nusep ) {
242  double beta[3], p_dcm_nu[4], p_nu[3], p_pcm_mp[3], partial;
243 
244  // Boost neu neutrino to mu decay CM
245  beta[0] = decay.pdpx / parent_energy;
246  beta[1] = decay.pdpy / parent_energy;
247  beta[2] = decay.pdpz / parent_energy;
248  p_nu[0] = (xpos-decay.vx)*enu/rad;
249  p_nu[1] = (ypos-decay.vy)*enu/rad;
250  p_nu[2] = (zpos-decay.vz)*enu/rad;
251  partial = gamma *
252  (beta[0]*p_nu[0] + beta[1]*p_nu[1] + beta[2]*p_nu[2] );
253  partial = enu - partial/(gamma+1.0);
254  // the following calculation is numerically imprecise
255  // especially p_dcm_nu[2] leads to taking the difference of numbers
256  // of order ~10's and getting results of order ~0.02's
257  // for g3numi we're starting with floats (ie. good to ~1 part in 10^7)
258  p_dcm_nu[0] = p_nu[0] - beta[0]*gamma*partial;
259  p_dcm_nu[1] = p_nu[1] - beta[1]*gamma*partial;
260  p_dcm_nu[2] = p_nu[2] - beta[2]*gamma*partial;
261  p_dcm_nu[3] = TMath::Sqrt( p_dcm_nu[0]*p_dcm_nu[0] +
262  p_dcm_nu[1]*p_dcm_nu[1] +
263  p_dcm_nu[2]*p_dcm_nu[2] );
264 
265  // Boost parent of mu to mu production CM
266  double particle_energy = decay.ppenergy;
267  gamma = particle_energy/parent_mass;
268  beta[0] = decay.ppdxdz * decay.pppz / particle_energy;
269  beta[1] = decay.ppdydz * decay.pppz / particle_energy;
270  beta[2] = decay.pppz / particle_energy;
271  partial = gamma * ( beta[0]*decay.muparpx +
272  beta[1]*decay.muparpy +
273  beta[2]*decay.muparpz );
274  partial = decay.mupare - partial/(gamma+1.0);
275  p_pcm_mp[0] = decay.muparpx - beta[0]*gamma*partial;
276  p_pcm_mp[1] = decay.muparpy - beta[1]*gamma*partial;
277  p_pcm_mp[2] = decay.muparpz - beta[2]*gamma*partial;
278  double p_pcm = TMath::Sqrt ( p_pcm_mp[0]*p_pcm_mp[0] +
279  p_pcm_mp[1]*p_pcm_mp[1] +
280  p_pcm_mp[2]*p_pcm_mp[2] );
281 
282  const double eps = 1.0e-30; // ? what value to use
283  if ( p_pcm < eps || p_dcm_nu[3] < eps ) {
284  return 3; // mu missing parent info?
285  }
286  // Calc new decay angle w.r.t. (anti)spin direction
287  double costh = ( p_dcm_nu[0]*p_pcm_mp[0] +
288  p_dcm_nu[1]*p_pcm_mp[1] +
289  p_dcm_nu[2]*p_pcm_mp[2] ) /
290  ( p_dcm_nu[3]*p_pcm );
291  // protect against small excursions
292  if ( costh > 1.0 ) costh = 1.0;
293  if ( costh < -1.0 ) costh = -1.0;
294  // Calc relative weight due to angle difference
295  double wgt_ratio = 0.0;
296  switch ( decay.ntype ) {
297  case kpdg_nue:
298  case kpdg_nuebar:
299  wgt_ratio = 1.0 - costh;
300  break;
301  case kpdg_numu:
302  case kpdg_numubar:
303  {
304  double xnu = 2.0 * enuzr / kMUMASS;
305  wgt_ratio = ( (3.0-2.0*xnu ) - (1.0-2.0*xnu)*costh ) / (3.0-2.0*xnu);
306 
307  if ( wgt_ratio < 0.0 ) {
308  std::cerr << "bsim::calcEnuWgt encountered serious problem: "
309  << " wgt_ratio " << wgt_ratio
310  << " enu " << enu << " costh " << costh << " xnu " << xnu
311  << " enuzr=decay.necm " << enuzr << " kMUMASS " << kMUMASS
312  << " norig " << decay.norig
313  << " ndecay " << decay.ndecay
314  << " ntype " << decay.ntype
315  << " ptype " << decay.ptype
316  << std::endl;
317  enu = 0;
318  wgt_xy = 0;
319  return 4; // bad, bad, bad calculation
320  }
321  break;
322  }
323  default:
324  enu = 0.0;
325  wgt_xy = 0.0;
326  return 2; // bad neutrino type for muon decay
327  }
328 
329  wgt_xy = wgt_xy * wgt_ratio;
330 
331  } // ptype is muon
332 
333  return 0;
334 }
335 //___________________________________________________________________________
336 
337 int bsim::calcEnuWgt(const bsim::Dk2Nu* dk2nu, const TVector3& xyz,
338  double& enu, double& wgt_xy)
339 {
340  return bsim::calcEnuWgt(dk2nu->decay,xyz,enu,wgt_xy);
341 }
342 //___________________________________________________________________________
Double_t ppdydz
% direction of nu parent at its production point
Definition: dk2nu.h:139
std::vector< bsim::Location > location
locations
Definition: dkmeta.h:119
Double_t muparpx
%
Definition: dk2nu.h:160
Double_t pppz
% z momentum of nu parent at its production point
Definition: dk2nu.h:140
int status
Definition: fabricate.py:1613
Double_t pdpx
% px momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:133
std::vector< bsim::NuRay > nuray
rays through detector fixed points
Definition: dk2nu.h:328
Double_t necm
% nu energy in center-of-mass frame
Definition: dk2nu.h:165
OStream cerr
Definition: OStream.cxx:7
bsim::Decay decay
basic decay information
Definition: dk2nu.h:327
Double_t beta
bsim::Dk2Nu * dk2nu
Double_t ppdxdz
% direction of nu parent at its production point
Definition: dk2nu.h:138
mu- => nu_mu + nu_e_bar + e-
Definition: dk2nu.h:82
Definition: lz4.cxx:387
Double_t vz
% neutrino production vertex z
Definition: dk2nu.h:132
double enu
Definition: runWimpSim.h:113
Double_t pdpz
% pz momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:135
Double_t muparpy
%
Definition: dk2nu.h:161
Double_t muparpz
%
Definition: dk2nu.h:162
Double_t ppenergy
% energy of nu parent at its production point
Definition: dk2nu.h:141
Int_t ntype
% neutrino flavor (PDG? code)
Definition: dk2nu.h:128
TVector3 Unit() const
static constexpr Double_t rad
Definition: Munits.h:162
int calcEnuWgt(const bsim::Decay &decay, const TVector3 &xyz, double &enu, double &wgt_xy)
workhorse routine
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
bsim::DkMeta * dkmeta
Definition: decay.py:1
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
mu+ => nu_mu_bar + nu_e + e+
Definition: dk2nu.h:81
Double_t vx
% neutrino production vertex x
Definition: dk2nu.h:130
void calcLocationWeights(const bsim::DkMeta *dkmeta, bsim::Dk2Nu *dk2nu)
user interface
Int_t ndecay
decay process (see dkproc_t)
Definition: dk2nu.h:127
Double_t mupare
% energy of nu grandparent
Definition: dk2nu.h:163
Double_t vy
% neutrino production vertex y
Definition: dk2nu.h:131
Int_t norig
not used?
Definition: dk2nu.h:126
Double_t pdpy
% py momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:134
Int_t ptype
% nu parent species (PDG? code)
Definition: dk2nu.h:144