Functions
kinematics.C File Reference
#include "code/singlekaon_xsec.cxx"
#include <TMath.h>
#include <TVector3.h>

Go to the source code of this file.

Functions

void kinematics ()
 

Function Documentation

void kinematics ( )

Definition at line 10 of file kinematics.C.

References abs(), std::acos(), std::cos(), om::cout, dE, rsim::DEBUG, singlekaon_xsec::diffxsec(), e, make_cached_def::e1, make_cached_def::e2, Emax, allTimeWatchdog::endl, evt, singlekaon_xsec::init(), python.hepunit::pi, std_candles::pl, printf(), central_limit::rand, canMan::sign(), std::sin(), std::sqrt(), string, dump_to_simple_cpp::threshold, Unit(), DeclareReco::VERBOSE, and xsec.

Referenced by genie::QPMDISStrucFuncBase::Calculate(), genie::COHHadronicSystemGenerator::CalculateHadronicSystem_AlvarezRuso(), genie::SKHadronicSystemGenerator::CalculateHadronicSystem_AtharSingleKaon(), genie::COHPrimaryLeptonGenerator::CalculatePrimaryLepton_AlvarezRuso(), ClassImp(), genie::QELEventGeneratorSM::ComputeMaxXSec(), genie::QELEventGeneratorSM::ComputeMaxXSec2(), genie::SmithMonizQELCCPXSec::d2sQES_dQ2dv_SM(), genie::SmithMonizQELCCPXSec::d3sQES_dQ2dvdkF_SM(), genie::utils::gsl::d3Xsec_dTldTkdCosThetal::DoEval(), genie::utils::gsl::d2Xsec_dTCosth::DoEval(), genie::utils::gsl::d2Xsec_dQ2dv::DoEval(), genie::utils::gsl::d5XSecAR::DoEval(), genie::utils::gsl::d5Xsec_dEldOmegaldOmegapi::DoEval(), genie::utils::gsl::d4Xsec_dEldThetaldOmegapi::DoEval(), genie::utils::gsl::d3Xsec_dOmegaldThetapi::DoEval(), genie::KovalenkoQELCharmPXSec::DR(), genie::SmithMonizQELCCPXSec::dsQES_dQ2_SM(), genie::BergerSehgalFMCOHPiPXSec2015::ExactKinematicTerm(), genie::BergerSehgalCOHPiPXSec2015::ExactKinematicTerm(), genie::LwlynSmithQELCCPXSec::FullDifferentialXSec(), genie::PythiaHadronization::Hadronize(), genie::CharmHadronization::Hadronize(), genie::PythiaHadronization::MultiplicityProb(), genie::BergerSehgalFMCOHPiPXSec2015::PionCOMAbsMomentum(), genie::BergerSehgalCOHPiPXSec2015::PionCOMAbsMomentum(), genie::QELEventGeneratorSM::ProcessEventRecord(), genie::QPMDISStrucFuncBase::Q2(), genie::utils::kinematics::Q2(), genie::MECGenerator::SelectNSVLeptonKinematics(), genie::LwlynSmithFF::tau(), genie::utils::kinematics::W(), genie::COHElasticPXSec::XSec(), genie::EmpiricalMECPXSec2015::XSec(), genie::SlowRsclCharmDISPXSecLO::XSec(), genie::AlamSimoAtharVacasSKPXSec2014::XSec(), genie::AhrensNCELPXSec::XSec(), genie::ReinDFRPXSec::XSec(), genie::RosenbluthPXSec::XSec(), genie::QPMDISPXSec::XSec(), genie::AlvarezRusoCOHPiPXSec::XSec(), genie::AivazisCharmPXSecLO::XSec(), genie::IMDAnnihilationPXSec::XSec(), genie::AhrensDMELPXSec::XSec(), genie::P33PaschosLalakulichPXSec::XSec(), genie::BergerSehgalFMCOHPiPXSec2015::XSec(), genie::QPMDMDISPXSec::XSec(), genie::LwlynSmithQELCCPXSec::XSec(), genie::BergerSehgalCOHPiPXSec2015::XSec(), genie::NuElectronPXSec::XSec(), genie::PaisQELLambdaPXSec::XSec(), genie::ReinSehgalCOHPiPXSec::XSec(), genie::KovalenkoQELCharmPXSec::XSec(), genie::ReinSehgalRESPXSec::XSec(), genie::BSKLNBaseRESPXSec2014::XSec(), and genie::NievesQELCCPXSec::XSec().

10  {
11  // Validate kinematics of single kaon production
12  // Have: T_l, T_K, cos(theta_l), phi_l, phi_Kq (?)
13  // Want: p_l, p_K, p_N (3-vectors)
14 
15  // Specify amount of output
16  const int VERBOSE = 0; // verbosity (0/1)
17  const int DEBUG = 0; // debugging mode
18 
19  // Use random input variables (0: use fixed inputs to verify GENIE results)
20  const int RANDOM = 1;
21 
22  // Initialise random number seed
23  srand (time(NULL));
24  const double pi = TMath::Pi();
25 
26  const int type = 2; // lepton: 1=electron, 2=muon, 3=tau
27  const int reac = 3; // reaction: 1=NN, 2=NP, 3=PP
28 
29  std::string strl;
30  double ml = 0.; // lepton mass [GeV]
31  if (type==1) { ml = 0.510998928e-3; strl = "e"; }
32  else if (type==2) { ml = 0.1056583715; strl = "mu"; }
33  else if (type==3) { ml = 1.77682; strl = "tau"; }
34  else {std::cout<<"ERROR: Invalid lepton type!"<<std::endl; return;}
35 
36  std::string strK;
37  double mK = 0.; // kaon mass [GeV]
38  if (reac==1) { mK = 0.493677; strK = "K+"; }
39  else if (reac==2) { mK = 0.497614; strK = "K0"; }
40  else if (reac==3) { mK = 0.493677; strK = "K+"; }
41  else {std::cout<<"ERROR: Invalid reaction!"<<std::endl; return;}
42 
43  std::string strN0, strN1;
44  double mN = 0.; // nucleon mass [GeV]
45  if (reac==1) { mN = 0.939565379; strN0 = "n"; strN1 = "n"; }
46  else if (reac==2) { mN = 0.939565379; strN0 = "n"; strN1 = "p"; }
47  else if (reac==3) { mN = 0.938272046; strN0 = "p"; strN1 = "p"; }
48  else {std::cout<<"ERROR: Invalid reaction!"<<std::endl; return;}
49 
50  std::string strReac = Form( "nu_%s + %s -> %s + %s + %s", strl.c_str(), strN0.c_str(),
51  strl.c_str(), strN1.c_str(), strK.c_str() );
52 
53  // Reaction threshold
54  const double threshold = ((ml+mN+mK)*(ml+mN+mK)-mN*mN)/(2.0*mN);
55  const double Emax = 3.1;
56 
57  // Neutrino energy [GeV] - chosen at random between threshold and Emax
58  double Enu;
59  if (RANDOM) Enu = threshold + (Emax-threshold)*(rand()/(double)RAND_MAX);
60  else Enu = 1.0;
61 
62  // Initialize reaction
63  singlekaon_xsec *inst = new singlekaon_xsec();
64  inst->init(Enu, type, reac);
65 
66  // INPUT PARAMETERS
67  double TK_max, TK, Tl_max, Tl, costh_l, phi_l, phi_Kq, xsec;
68  int evt;
69  if (RANDOM) {
70  xsec = -999.;
71  while(xsec <= 0) {
72  TK_max = Enu - mK - ml; // maximal allowed kaon kinetic energy
73  TK = TK_max*(rand()/(double)RAND_MAX); // kaon kinetic energy [0, TK_max]
74  Tl_max = Enu - mK - ml - TK; // maximal allowed lepton kinetic energy
75  Tl = Tl_max*(rand()/(double)RAND_MAX); // lepton kinetic energy [0, Tl_max]
76  costh_l = 2.*(rand()/(double)RAND_MAX)-1.; // lepton polar angle [-1, 1]
77  phi_l = pi*(2.*(rand()/(double)RAND_MAX)-1.); // lepton azimuthal angle [-pi, pi]
78  phi_Kq = pi*(2.*(rand()/(double)RAND_MAX)-1.); // kaon azimuthal angle [-pi, pi]
79  xsec = inst->diffxsec(Tl,TK,acos(costh_l),phi_Kq); // 4D-differential cross-section
80  }
81  } else {
82  xsec = -999.;
83  evt = 0; // choose event from input table
84  double input[9][5] = {
85  {0.188312, 0.108214, 0.706911, -0.66872, 5.75773},
86  {0.234579, 0.0893942, 0.74982, 3.035, 6.10204},
87  {0.196701, 0.0362663, 0.663934, 2.64581, 4.82931},
88  {0.0191452, 0.0575003, 0.67804, -2.96401, 4.11426},
89  {0.0529443, 0.0119084, 0.127732, -1.57085, 1.10203},
90  {0.139563, 0.143596, 0.946546, 0.80543, 2.34818},
91  {0.200567, 0.13117, 0.851108, -0.806353, 1.98515},
92  {0.0372796, 0.0501203, 0.148006, 0.81699, 4.35815},
93  {0.0806993, 0.172774, 0.903585, 2.42196, 3.6006}
94  };
95 
96  TK = input[evt][0];
97  Tl = input[evt][1];
98  costh_l = input[evt][2];
99  phi_l = input[evt][3];
100  phi_Kq = input[evt][4];
101  xsec = inst->diffxsec(Tl,TK,acos(costh_l),phi_Kq);
102  if (xsec <= 0) {
103  printf("ERROR: INVALID XSEC! (%lf)\n", xsec);
104  return;
105  }
106  }
107 
108  // Print input parameters to screen
109  printf("\n\n");
110  printf("INITIAL INPUT PARAMETERS\n");
111  printf("------------------------\n");
112  printf("Neutrino energy:\t\tE_nu\t= %6.3lf GeV\n",Enu);
113  printf("Lepton kinetic energy:\t\tT_l\t= %6.3lf GeV\n",Tl);
114  printf("Kaon kinetic energy:\t\tT_K\t= %6.3lf GeV\n",TK);
115  printf("Lepton polar angle:\t\tcosth_l\t= %6.3lf\n",costh_l);
116  printf("Lepton azimuthal angle:\t\tphi_l\t= %6.3lf pi\n",phi_l/pi);
117  printf("Kaon azimuthal angle:\t\tphi_Kq\t= %6.3lf pi\n",phi_Kq/pi);
118  printf("Resulting cross-section:\td4xsec\t= %10.3e [units]\n",xsec);
119  printf("\n\n");
120 
121  // Some directly obtainable parameters
122  const double theta_l = acos(costh_l);
123  const double TN = Enu - Tl - TK - ml - mK;
124  const double El = Tl + ml;
125  const double EK = TK + mK;
126  const double EN = TN + mN;
127 
128  // Momentum 3-vectors
129  TVector3 p_l, p_Q, p_K, p_N, p_Kq, p_Nq;
130 
131  // Unit 3-vectors
132  TVector3 e1(1,0,0);
133  TVector3 e2(0,1,0);
134  TVector3 e3(0,0,1);
135 
136  double pl, pK, pN; // Total momentum
137  double bl, bK, bN; // Lorentz beta
138  double gl, gK, gN; // Lorentz gamma
139 
140  // Lepton momentum
141  gl = Tl/ml + 1;
142  bl = sqrt(1.-1./gl/gl);
143  pl = bl*gl*ml;
144 
145  // Kaon momentum
146  gK = TK/mK + 1;
147  bK = sqrt(1.-1./gK/gK);
148  pK = bK*gK*mK;
149 
150  // Nucleon momentum
151  gN = TN/mN + 1;
152  bN = sqrt(1.-1./gN/gN);
153  pN = bN*gN*mN;
154 
155  // Calculate 3-momentum of lepton
156  p_l[0] = pl*sin(theta_l)*cos(phi_l);
157  p_l[1] = pl*sin(theta_l)*sin(phi_l);
158  p_l[2] = pl*costh_l;
159  if (VERBOSE) {
160  printf("Lepton 3-mom in lab:\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
161  p_l[0], p_l[1], p_l[2], pl);
162  }
163 
164  // Calculate 3-momentum of momentum transfer to hadronic system
165  p_Q = -p_l;
166  p_Q[2] += Enu;
167  double pQ = p_Q.Mag();
168  if (VERBOSE) {
169  printf("Q-squared 3-mom in lab:\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
170  p_Q[0], p_Q[1], p_Q[2], pQ);
171  printf("\nQ-squared 3-mom in Q2:\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
172  0., 0., pQ, pQ);
173  }
174  TVector3 dirQ = p_Q.Unit();
175 
176  // In the momentum transfer plane, new angles... (see eq.17 of notes)
177  double costh_Kq = (pQ*pQ + pK*pK + mN*mN - (Enu-El-EK+mN)*(Enu-El-EK+mN))/(2.*pQ*pK);
178  double theta_Kq = acos(costh_Kq);
179 
180  // Get rotation angle from z-axis to Q^2-direction
181  double rot_Z = acos(e3*(p_Q.Unit())); // rotate counter-clockwise, towards x-axis by [0,pi]
182 
183  // Get rotation angle from x-axis to Q^2-direction (projected to transverse plane)
184  int sign = 0;
185  if (p_Q[1] > 0) sign = +1;
186  else sign = -1;
187  TVector3 p_Qt(p_Q[0],p_Q[1],0);
188  double rot_X = sign*acos(e1*(p_Qt.Unit())); // rotate around z-axis by [-pi,pi]
189 
190  // Make sure that x' is in the nu-q plane
191  TVector3 p_nu = Enu*e3;
192  TVector3 yrot = (p_Q.Cross(p_nu)).Unit(); // must be parallel to y' axis
193  p_nu.RotateZ(-rot_X);
194  p_nu.RotateY(-rot_Z);
195  if (VERBOSE) {
196  printf("Neutrino 3-mom in Q2:\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c",
197  p_nu[0], p_nu[1], p_nu[2], p_nu.Mag());
198  if ( abs(p_nu[1])>1.e-6 )
199  printf("--> WARNING: NOT IN X'Z'-PLANE!\n");
200  else
201  printf("\n");
202  printf("\n");
203  }
204 
205  // Make sure that q points the same way as the y'-axis
206  TVector3 unu = p_nu.Unit();
207  TVector3 uQ = p_Q.Unit();
208  uQ.RotateZ(-rot_X);
209  uQ.RotateY(-rot_Z);
210  if (DEBUG) {
211  printf("Unit vector nu in Q2:\t(%6.3lf | %6.3lf | %6.3lf )\n",unu[0],unu[1],unu[2]);
212  printf("Unit vector q in Q2:\t(%6.3lf | %6.3lf | %6.3lf )",uQ[0],uQ[1],uQ[2]);
213  if ( abs(uQ[0])>1.e-6 || abs(uQ[1])>1.e-6 || abs(uQ[2])<(1-1.e-6) )
214  printf(" --> WARNING: NOT EQUAL TO Z' UNIT VECTOR!\n");
215  else
216  printf("\n");
217  printf("\n");
218  }
219 
220  // Check what happens to the unit vectors when rotated (for debugging purposes)
221  if (DEBUG) {
222  TVector3 u1=e1, u2=e2, u3=e3; // unit vectors
223  u1.RotateZ(-rot_X); u2.RotateZ(-rot_X); u3.RotateZ(-rot_X);
224  u1.RotateY(-rot_Z); u2.RotateY(-rot_Z); u3.RotateY(-rot_Z);
225  printf("Unit vector X in Q2:\t(%6.3lf | %6.3lf | %6.3lf )\n",u1[0],u1[1],u1[2]);
226  printf("Unit vector Y in Q2:\t(%6.3lf | %6.3lf | %6.3lf )\n",u2[0],u2[1],u2[2]);
227  printf("Unit vector Z in Q2:\t(%6.3lf | %6.3lf | %6.3lf )\n",u3[0],u3[1],u3[2]);
228  printf("\n");
229  u1=e1; u2=e2; u3=e3;
230 
231  u1.RotateY(rot_Z); u2.RotateY(rot_Z); u3.RotateY(rot_Z);
232  u1.RotateZ(rot_X); u2.RotateZ(rot_X); u3.RotateZ(rot_X);
233  printf("Unit vector X' in lab:\t(%6.3lf | %6.3lf | %6.3lf )\n",u1[0],u1[1],u1[2]);
234  printf("Unit vector Y' in lab:\t(%6.3lf | %6.3lf | %6.3lf )\n",u2[0],u2[1],u2[2]);
235  printf("Unit vector Z' in lab:\t(%6.3lf | %6.3lf | %6.3lf )\n",u3[0],u3[1],u3[2]);
236  printf("\n");
237  printf("CHECK yrot in lab:\t(%6.3lf | %6.3lf | %6.3lf )\n",yrot[0],yrot[1],yrot[2]);
238  printf("CHECK q1*x2' = q2*x1':\t%6.3lf = %6.3lf\n", p_Q[0]*u1[1], p_Q[1]*u1[0]);
239  printf("\n");
240  }
241 
242  // Calculate 3-momentum of kaon (in Q-frame)
243  p_Kq[0] = pK*sin(theta_Kq)*cos(phi_Kq);
244  p_Kq[1] = pK*sin(theta_Kq)*sin(phi_Kq);
245  p_Kq[2] = pK*costh_Kq;
246 
247  // Calculate 3-momentum of nucleon (in Q-frame)
248  p_Nq = -p_Kq;
249  p_Nq[2] += pQ;
250  if (VERBOSE) {
251  printf("Kaon 3-mom in Q2:\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
252  p_Kq[0], p_Kq[1], p_Kq[2], p_Kq.Mag());
253  printf("Nucleon 3-mom in Q2:\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n\n",
254  p_Nq[0], p_Nq[1], p_Nq[2], p_Nq.Mag());
255  }
256 
257  // Rotate particles into correct frame
258  p_K = p_Kq;
259  p_K.RotateY(rot_Z);
260  p_K.RotateZ(rot_X);
261  p_N = p_Nq;
262  p_N.RotateY(rot_Z);
263  p_N.RotateZ(rot_X);
264  if (VERBOSE) {
265  printf("Kaon 3-mom in lab:\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
266  p_K[0], p_K[1], p_K[2], p_K.Mag());
267  printf("Nucleon 3-mom in lab:\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
268  p_N[0], p_N[1], p_N[2], p_N.Mag());
269  printf("\n\n");
270  }
271 
272  // Compare alternative rotation by Chris
273  if (DEBUG) {
274  TVector3 p_K2 = p_Kq;
275  p_K2.RotateUz( p_Q.Unit() );
276  TVector3 p_N2 = p_Nq;
277  p_N2.RotateUz( p_Q.Unit() );
278  printf("Kaon 3-mom in Q2:\t\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
279  p_Kq[0], p_Kq[1], p_Kq[2], p_Kq.Mag());
280  printf("Kaon 3-mom in lab (Martti):\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
281  p_K[0], p_K[1], p_K[2], p_K.Mag());
282  printf("Kaon 3-mom in lab (Chris):\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
283  p_K2[0], p_K2[1], p_K2[2], p_K2.Mag());
284  printf("\n");
285  printf("Nucleon 3-mom in Q2:\t\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
286  p_Nq[0], p_Nq[1], p_Nq[2], p_Nq.Mag());
287  printf("Nucleon 3-mom in lab (Martti):\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
288  p_N[0], p_N[1], p_N[2], p_N.Mag());
289  printf("Nucleon 3-mom in lab (Chris):\t(%6.3lf | %6.3lf | %6.3lf ) --> %6.3lf GeV/c\n",
290  p_N2[0], p_N2[1], p_N2[2], p_N2.Mag());
291  printf("\n\n");
292  }
293 
294  // Check total momenta
295  if (DEBUG) {
296  printf("Check lepton momentum:\t%5.3lf = %5.3lf\n",pl,p_l.Mag());
297  printf("Check nucleon momentum:\t%5.3lf = %5.3lf\n",pN,p_N.Mag());
298  printf("Check kaon momentum:\t%5.3lf = %5.3lf\n",pK,p_K.Mag());
299  printf("\n\n");
300  }
301 
302  if (!RANDOM) {
303  printf("COMPARE THIS TO EVENT #%d FROM CHRIS! (email 30.10.2014)\n\n",evt+1);
304  return;
305  }
306 
307  // Check that [numu -> mu+Q] conserves energy+momentum
308  // -----------------------------------------------------
309  if (VERBOSE) {
310  double dX1 = p_l[0]+p_Q[0];
311  double dY1 = p_l[1]+p_Q[1];
312  double dZ1 = p_l[2]+p_Q[2] - Enu;
313  double dE1 = (Enu-El + El) - (Enu);
314 
315  printf("CONSERVATION OF REACTION (nu -> l + q)\n");
316  printf("---------------------------------------------------------\n");
317  printf(" Particle | px py pz E m \n");
318  printf("---------------------------------------------------------\n");
319  printf(" neutrino | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",0.,0.,Enu,Enu,0.);
320  printf("---------------------------------------------------------\n");
321  printf(" lepton | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",p_l[0],p_l[1],p_l[2],El,ml);
322  printf(" Q^2 | %6.3lf %6.3lf %6.3lf %6.3lf --- \n",p_Q[0],p_Q[1],p_Q[2],Enu-El);
323  printf("---------------------------------------------------------\n");
324  printf(" FIN-INIT | %6.3lf %6.3lf %6.3lf %6.3lf \n",dX1,dY1,dZ1,dE1);
325  printf("---------------------------------------------------------\n");
326  printf("\n\n");
327  }
328 
329  // Check that [Q+N -> K+N] conserves energy+momentum
330  // -------------------------------------------------
331  if (VERBOSE) {
332  double dX2 = p_Kq[0]+p_Nq[0] - pQ*e3[0];
333  double dY2 = p_Kq[1]+p_Nq[1] - pQ*e3[1];
334  double dZ2 = p_Kq[2]+p_Nq[2] - pQ*e3[2];
335  double dE2 = (EN+EK) - (Enu-El+mN);
336 
337  printf("CONSERVATION OF REACTION (q + N -> K + N) - ROTATED FRAME\n");
338  printf("---------------------------------------------------------\n");
339  printf(" Particle | px py pz E m \n");
340  printf("---------------------------------------------------------\n");
341  printf(" Q^2 | %6.3lf %6.3lf %6.3lf %6.3lf --- \n",pQ*e3[0],pQ*e3[1],pQ*e3[2],Enu-El);
342  printf(" nucleon | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",0.,0.,0.,mN,mN);
343  printf("---------------------------------------------------------\n");
344  printf(" kaon | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",p_Kq[0],p_Kq[1],p_Kq[2],EK,mK);
345  printf(" nucleon | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",p_Nq[0],p_Nq[1],p_Nq[2],EN,mN);
346  printf("---------------------------------------------------------\n");
347  printf(" FIN-INIT | %6.3lf %6.3lf %6.3lf %6.3lf \n",dX2,dY2,dZ2,dE2);
348  printf("---------------------------------------------------------\n");
349  printf("\n\n");
350  }
351 
352  // Check that [Q+N -> K+N] conserves energy+momentum
353  // -------------------------------------------------
354  if (VERBOSE) {
355  double dX3 = p_K[0]+p_N[0] - p_Q[0];
356  double dY3 = p_K[1]+p_N[1] - p_Q[1];
357  double dZ3 = p_K[2]+p_N[2] - p_Q[2];
358  double dE3 = (EN+EK) - (Enu-El+mN);
359 
360  printf("CONSERVATION OF REACTION (q + N -> K + N) - NORMAL FRAME\n");
361  printf("---------------------------------------------------------\n");
362  printf(" Particle | px py pz E m \n");
363  printf("---------------------------------------------------------\n");
364  printf(" Q^2 | %6.3lf %6.3lf %6.3lf %6.3lf --- \n",p_Q[0],p_Q[1],p_Q[2],Enu-El);
365  printf(" nucleon | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",0.,0.,0.,mN,mN);
366  printf("---------------------------------------------------------\n");
367  printf(" kaon | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",p_K[0],p_K[1],p_K[2],EK,mK);
368  printf(" nucleon | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",p_N[0],p_N[1],p_N[2],EN,mN);
369  printf("---------------------------------------------------------\n");
370  printf(" FIN-INIT | %6.3lf %6.3lf %6.3lf %6.3lf \n",dX3,dY3,dZ3,dE3);
371  printf("---------------------------------------------------------\n");
372  printf("\n\n");
373  }
374 
375  // Check that entire reaction conserves energy+momentum
376  // ----------------------------------------------------
377  double dX = p_l[0]+p_N[0]+p_K[0];
378  double dY = p_l[1]+p_N[1]+p_K[1];
379  double dZ = p_l[2]+p_N[2]+p_K[2] - Enu;
380  double dE = (El+EN+EK) - (Enu+mN);
381 
382  printf("KINEMATICS OF PARTICLES IN REACTION (%s)\n",strReac.c_str());
383  printf("---------------------------------------------------------\n");
384  printf(" Particle | px py pz E m \n");
385  printf("---------------------------------------------------------\n");
386  printf(" neutrino | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",0.,0.,Enu,Enu,0.);
387  printf(" nucleon | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",0.,0.,0.,mN,mN);
388  printf("---------------------------------------------------------\n");
389  printf(" lepton | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",p_l[0],p_l[1],p_l[2],El,ml);
390  printf(" nucleon | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",p_N[0],p_N[1],p_N[2],EN,mN);
391  printf(" kaon | %6.3lf %6.3lf %6.3lf %6.3lf %6.3lf \n",p_K[0],p_K[1],p_K[2],EK,mK);
392  printf("---------------------------------------------------------\n");
393  printf(" FIN-INIT | %6.3lf %6.3lf %6.3lf %6.3lf \n",dX,dY,dZ,dE);
394  printf("---------------------------------------------------------\n");
395  printf("\n\n");
396 
397  const double EPS = 1e-6;
398  if (abs(dX)<EPS && abs(dY)<EPS && abs(dZ)<EPS && abs(dE)<EPS)
399  printf("INFO: ALL OK - energy & momentum are conserved.\n\n");
400  else
401  printf("WARNING: something is wrong, check E/p conservation!\n\n");
402  return;
403 }
DeclareReco.VERBOSE
bool VERBOSE
Definition: DeclareReco.py:8
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
dump_to_simple_cpp.threshold
threshold
Definition: dump_to_simple_cpp.py:7
singlekaon_xsec
Definition: singlekaon_xsec.h:6
make_cached_def.e2
e2
Definition: make_cached_def.py:53
abs
void abs(TH1 *hist)
Definition: wrong_sign_uncertainty.C:716
std::acos
T acos(T number)
Definition: d0nt_math.hpp:54
xml_schema::time
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
dE
double dE
Definition: write_out_h20_xsec.C:9
evt
int evt
Definition: FillTrainingTrees.h:22
printf
printf("%d Experimental points found\n", nlines)
singlekaon_xsec::diffxsec
double diffxsec(double Tlep, double Tkaon, double theta, double phikq)
Definition: singlekaon_xsec.cxx:90
make_cached_def.e1
e1
Definition: make_cached_def.py:40
Unit
TVector3 Unit() const
Definition: SRVector3DAddons.h:14
xsec
Double_t xsec[nknots]
Definition: testXsec.C:47
om::cout
OStream cout
Definition: OStream.cxx:6
std::sin
T sin(T number)
Definition: d0nt_math.hpp:132
std::cos
T cos(T number)
Definition: d0nt_math.hpp:78
Emax
double Emax
Definition: write_out_h20_xsec.C:8
singlekaon_xsec::init
void init(double Etot, int type, int reac)
Definition: singlekaon_xsec.cxx:3
e
Float_t e
Definition: plot.C:35
central_limit.rand
rand
Definition: central_limit.py:20
std_candles.pl
pl
Definition: std_candles.py:675
canMan.sign
def sign(x)
Definition: canMan.py:197
string
enum BeamMode string
Generated on Wed Nov 25 2020 05:17:39 for NOvA by   doxygen 1.8.11