24 const double pi = TMath::Pi();
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"; }
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+"; }
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"; }
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() );
54 const double threshold = ((ml+mN+mK)*(ml+mN+mK)-mN*mN)/(2.0*mN);
55 const double Emax = 3.1;
59 if (RANDOM) Enu = threshold + (Emax-
threshold)*(
rand()/(double)RAND_MAX);
64 inst->
init(Enu, type, reac);
67 double TK_max, TK, Tl_max, Tl, costh_l, phi_l, phi_Kq,
xsec;
72 TK_max = Enu - mK - ml;
73 TK = TK_max*(
rand()/(double)RAND_MAX);
74 Tl_max = Enu - mK - ml - TK;
75 Tl = Tl_max*(
rand()/(double)RAND_MAX);
76 costh_l = 2.*(
rand()/(double)RAND_MAX)-1.;
77 phi_l = pi*(2.*(
rand()/(double)RAND_MAX)-1.);
78 phi_Kq = pi*(2.*(
rand()/(double)RAND_MAX)-1.);
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}
98 costh_l = input[
evt][2];
99 phi_l = input[
evt][3];
100 phi_Kq = input[
evt][4];
103 printf(
"ERROR: INVALID XSEC! (%lf)\n", xsec);
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);
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;
129 TVector3 p_l, p_Q, p_K, p_N, p_Kq, p_Nq;
142 bl =
sqrt(1.-1./gl/gl);
147 bK =
sqrt(1.-1./gK/gK);
152 bN =
sqrt(1.-1./gN/gN);
156 p_l[0] = pl*
sin(theta_l)*
cos(phi_l);
157 p_l[1] = pl*
sin(theta_l)*
sin(phi_l);
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);
167 double pQ = p_Q.Mag();
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",
174 TVector3 dirQ = p_Q.Unit();
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);
181 double rot_Z =
acos(e3*(p_Q.Unit()));
185 if (p_Q[1] > 0) sign = +1;
187 TVector3 p_Qt(p_Q[0],p_Q[1],0);
188 double rot_X = sign*
acos(
e1*(p_Qt.Unit()));
191 TVector3 p_nu = Enu*e3;
192 TVector3 yrot = (p_Q.Cross(p_nu)).
Unit();
193 p_nu.RotateZ(-rot_X);
194 p_nu.RotateY(-rot_Z);
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");
206 TVector3 unu = p_nu.Unit();
207 TVector3 uQ = p_Q.Unit();
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");
222 TVector3 u1=
e1, u2=
e2, u3=e3;
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]);
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]);
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]);
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;
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());
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());
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());
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());
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());
303 printf(
"COMPARE THIS TO EVENT #%d FROM CHRIS! (email 30.10.2014)\n\n",evt+1);
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);
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");
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);
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");
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);
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");
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);
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");
397 const double EPS = 1
e-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");
401 printf(
"WARNING: something is wrong, check E/p conservation!\n\n");
::xsd::cxx::tree::time< char, simple_type > time
printf("%d Experimental points found\n", nlines)
double diffxsec(double Tlep, double Tkaon, double theta, double phikq)
void init(double Etot, int type, int reac)