Public Member Functions | Public Attributes | Private Attributes | List of all members
NeutrinoFluxAuxiliar::NuWeight Class Reference

Get the weight to get the neutrino probability flux in one point. More...

#include <NuWeight.h>

Public Member Functions

 NuWeight (std::vector< double > &posdet)
 The constructor has to have the position of the detector. More...
 
 ~NuWeight ()
 Destructor. More...
 
void calculate_weight (bsim::Dk2Nu *nu)
 Calculate weight based on dk2nu entry. After this, enu and wgt have meaningful values. More...
 

Public Attributes

double enu
 
double wgt
 

Private Attributes

TDatabasePDG * particle
 
double xdet
 
double ydet
 
double zdet
 

Detailed Description

Get the weight to get the neutrino probability flux in one point.

Definition at line 19 of file NuWeight.h.

Constructor & Destructor Documentation

NeutrinoFluxAuxiliar::NuWeight::NuWeight ( std::vector< double > &  posdet)

The constructor has to have the position of the detector.

Definition at line 7 of file NuWeight.cpp.

7  {
8 
9  particle = TDatabasePDG::Instance();
10  NuWeight::enu = -1.0;
11  NuWeight::wgt = -1.0;
12  xdet = posdet[0];
13  ydet = posdet[1];
14  zdet = posdet[2];
15  }
NeutrinoFluxAuxiliar::NuWeight::~NuWeight ( )

Destructor.

Definition at line 118 of file NuWeight.cpp.

118  {
119  }

Member Function Documentation

void NeutrinoFluxAuxiliar::NuWeight::calculate_weight ( bsim::Dk2Nu *  nu)

Calculate weight based on dk2nu entry. After this, enu and wgt have meaningful values.

Definition at line 16 of file NuWeight.cpp.

16  {
17  //assumes units are GeV and cm
18 
19  int likely_ion = 100000000;
20  double mpar=-1.;
21  int size = (nu->ancestor).size();
22  int pdg = nu->ancestor[size-2].pdg;
23 
24  if(pdg>likely_ion || pdg<-1*likely_ion){
25  std::cout<< "Not handling ions yet, pdg: "<<pdg<<std::endl;
26  exit (1);
27  }
28  else{
29  mpar = particle->GetParticle(pdg)->Mass();
30  }
31 
32  double ppar = sqrt( pow(nu->decay.pdpx,2) + pow(nu->decay.pdpy,2) + pow(nu->decay.pdpz,2) );
33  double epar = sqrt(pow(ppar,2)+pow(mpar,2));
34  double gamma = epar / mpar;
35  double beta = sqrt(pow(gamma,2)-1)/gamma;
36 
37  double rr = sqrt(pow(xdet-nu->decay.vx,2)+pow(ydet-nu->decay.vy,2)+pow(zdet-nu->decay.vz,2));
38  double cos_theta = ( (nu->decay.pdpx)*(xdet-nu->decay.vx) + (nu->decay.pdpy)*(ydet-nu->decay.vy) + (nu->decay.pdpz)*(zdet-nu->decay.vz) )/ (ppar*rr);
39 
40  if(cos_theta > 1 || cos_theta < -1){
41  std::cout<< "Cosine of neutrino not allowed: "<<cos_theta<<std::endl;
42  exit (1);
43  }
44  double MM = 1.0/(gamma*(1.0-beta*cos_theta));
45  double angdet = (pow(rdet,2) /pow(rr,2)/ 4.);
46  NuWeight::enu = MM*(nu->decay.necm);
47  NuWeight::wgt = angdet * pow(MM,2);
48 
49  //done for all except polarized muon
50  // in which case need to modify weight
51  if (pdg == 13 || pdg == -13){
52 
53  //boost new neutrino to mu decay cm
54  double vbeta[3];
55  vbeta[0] = nu->decay.pdpx / epar;
56  vbeta[1] = nu->decay.pdpy / epar;
57  vbeta[2] = nu->decay.pdpz / epar;
58 
59  double p_nu[3]; //nu momentum
60  p_nu[0] = (xdet- nu->decay.vx) * (NuWeight::enu) / rr;
61  p_nu[1] = (ydet- nu->decay.vy) * (NuWeight::enu) / rr;
62  p_nu[2] = (zdet- nu->decay.vz) * (NuWeight::enu) / rr;
63 
64  double partial = gamma*(vbeta[0]*p_nu[0]+vbeta[1]*p_nu[1]+vbeta[2]*p_nu[2]);
65  partial = (NuWeight::enu)-partial / (gamma+1.);
66 
67  double p_dcm_nu[4];
68  for (int i=0;i<3;i++) p_dcm_nu[i]=p_nu[i]-vbeta[i]*gamma*partial;
69  p_dcm_nu[3]=0.;
70  for (int i=0;i<3;i++) p_dcm_nu[3]+=p_dcm_nu[i]*p_dcm_nu[i];
71  p_dcm_nu[3]=sqrt(p_dcm_nu[3]);
72 
73  //boost parent of mu to mu production cm
74  gamma = nu->decay.ppenergy / mpar;
75  vbeta[0] = nu->decay.ppdxdz * nu->decay.pppz / nu->decay.ppenergy;
76  vbeta[1] = nu->decay.ppdydz * nu->decay.pppz / nu->decay.ppenergy;
77  vbeta[2] = nu->decay.pppz / nu->decay.ppenergy;
78  partial = gamma*(vbeta[0]*nu->decay.muparpx+vbeta[1]*nu->decay.muparpy+vbeta[2]*nu->decay.muparpz);
79  partial = nu->decay.mupare - partial / (gamma+1.);
80 
81  double p_pcm_mp[4];
82  p_pcm_mp[0]=nu->decay.muparpx-vbeta[0]*gamma*partial;
83  p_pcm_mp[1]=nu->decay.muparpy-vbeta[1]*gamma*partial;
84  p_pcm_mp[2]=nu->decay.muparpz-vbeta[2]*gamma*partial;
85  p_pcm_mp[3]=0.;
86  for (int i=0;i<3;i++) p_pcm_mp[3]+=p_pcm_mp[i]*p_pcm_mp[i];
87  p_pcm_mp[3]=sqrt(p_pcm_mp[3]);
88 
89  double wt_ratio = 1.;
90  //have to check p_pcm_mp
91  //it can be 0 if mupar..=0. (I guess muons created in target??)
92  if (p_pcm_mp[3] != 0. ) {
93  //calc new decay angle w.r.t. (anti)spin direction
94  double costh = (p_dcm_nu[0]*p_pcm_mp[0]+
95  p_dcm_nu[1]*p_pcm_mp[1]+
96  p_dcm_nu[2]*p_pcm_mp[2])/(p_dcm_nu[3]*p_pcm_mp[3]);
97 
98  if (costh>1.) costh = 1.;
99  else if (costh<-1.) costh = -1.;
100 
101  //calc relative weight due to angle difference
102  if(nu->decay.ntype == 12 || nu->decay.ntype == -12){
103  wt_ratio = 1.-costh;
104  }
105  else if(nu->decay.ntype == 14 || nu->decay.ntype == -14){
106  double mumass = particle->GetParticle(13)->Mass();
107  double xnu = 2.* nu->decay.necm / mumass;
108  wt_ratio = ( (3.-2.*xnu) - (1.-2.*xnu)*costh ) / (3.-2.*xnu);
109  } else {
110  std::cout << "NuWeight:: Bad neutrino type"<<std::endl;
111  }
112  }
113  NuWeight::wgt *= wt_ratio;
114  }
115 
116  }
pionana::beta
double beta(double KE, const simb::MCParticle *part)
Definition: TruthAnalyzer_module.cc:59
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
rdet
const double rdet
Definition: NuWeight.cpp:3
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
pionana::gamma
double gamma(double KE, const simb::MCParticle *part)
Definition: TruthAnalyzer_module.cc:55
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030

Member Data Documentation

double NeutrinoFluxAuxiliar::NuWeight::enu

Definition at line 30 of file NuWeight.h.

TDatabasePDG* NeutrinoFluxAuxiliar::NuWeight::particle
private

Definition at line 34 of file NuWeight.h.

double NeutrinoFluxAuxiliar::NuWeight::wgt

Definition at line 31 of file NuWeight.h.

double NeutrinoFluxAuxiliar::NuWeight::xdet
private

Definition at line 35 of file NuWeight.h.

double NeutrinoFluxAuxiliar::NuWeight::ydet
private

Definition at line 36 of file NuWeight.h.

double NeutrinoFluxAuxiliar::NuWeight::zdet
private

Definition at line 37 of file NuWeight.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11