KineUtils.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Authors: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  Afroditi Papadopoulou <apapadop \at mit.edu>
11  Massachusetts Institute of Technology
12 
13  Changes required to implement the GENIE Boosted Dark Matter module
14  were installed by Josh Berger (Univ. of Wisconsin)
15 */
16 //____________________________________________________________________________
17 
18 #include <cstdlib>
19 #include <limits>
20 #include <TMath.h>
21 
29 
30 using namespace genie;
31 using namespace genie::constants;
32 
33 //____________________________________________________________________________
35  const Interaction * const in, KinePhaseSpace_t ps)
36 {
37  double vol = 0;
38 
39  const KPhaseSpace & phase_space = in->PhaseSpace();
40 
41  switch(ps) {
42 
43  case(kPSQ2fE):
44  {
45  Range1D_t Q2 = phase_space.Limits(kKVQ2);
46  vol = Q2.max - Q2.min;
47  return vol;
48  break;
49  }
50  case(kPSq2fE):
51  {
52  Range1D_t q2 = phase_space.Limits(kKVq2);
53  vol = q2.min - q2.max;
54  return vol;
55  break;
56  }
57  case(kPSWfE):
58  {
59  Range1D_t W = phase_space.Limits(kKVW);
60  vol = W.max - W.min;
61  return vol;
62  break;
63  }
64  case(kPSWQ2fE):
65  {
66  Range1D_t W = phase_space.Limits(kKVW);
67  if(W.max<0) return 0;
68  const int kNW = 100;
69  double dW = (W.max-W.min)/(kNW-1);
70  Interaction interaction(*in);
71  for(int iw=0; iw<kNW; iw++) {
72  interaction.KinePtr()->SetW(W.min + iw*dW);
73  Range1D_t Q2 = interaction.PhaseSpace().Q2Lim_W();
74  double dQ2 = (Q2.max-Q2.min);
75  vol += (dW*dQ2);
76  }
77  return vol;
78  break;
79  }
80  case(kPSxyfE):
81  {
82  Range1D_t W = phase_space.Limits(kKVW);
83  if(W.max<0) return 0;
84 
85  const InitialState & init_state = in->InitState();
86  double Ev = init_state.ProbeE(kRfHitNucRest);
87  double M = init_state.Tgt().HitNucP4Ptr()->M();
88 
89  const int kNx = 100;
90  const int kNy = 100;
91  const double kminx = controls::kMinX;
92  const double kminy = controls::kMinY;
93  const double kdx = (controls::kMaxX - kminx) / (kNx-1);
94  const double kdy = (controls::kMaxY - kminy) / (kNy-1);
95  const double kdV = kdx*kdy;
96 
97  double cW=-1, cQ2 = -1;
98 
99  Interaction interaction(*in);
100 
101  for(int ix=0; ix<kNx; ix++) {
102  double x = kminx+ix*kdx;
103  for(int iy=0; iy<kNy; iy++) {
104  double y = kminy+iy*kdy;
105 
106  XYtoWQ2(Ev, M, cW, cQ2, x, y);
107  if(!math::IsWithinLimits(cW, W)) continue;
108 
109  interaction.KinePtr()->SetW(cW);
110  Range1D_t Q2 = interaction.PhaseSpace().Q2Lim_W();
111  if(!math::IsWithinLimits(cQ2, Q2)) continue;
112 
113  vol += kdV;
114  }
115  }
116  return vol;
117  break;
118  }
119  default:
120  SLOG("KineLimits", pERROR)
121  << "Couldn't compute phase space volume for "
122  << KinePhaseSpace::AsString(ps) << "using \n" << *in;
123  return 0;
124  }
125  return 0;
126 }
127 //____________________________________________________________________________
129  const Interaction * const i, KinePhaseSpace_t fromps, KinePhaseSpace_t tops)
130 {
131 // Returns the Jacobian for a kinematical transformation:
132 // from_ps (x,y,z,...) -> to_ps (u,v,w,...).
133 //
134 // Note on the convention used here:
135 // The differential cross-section d^n sigma / dxdydz..., is transformed to
136 // d^n sigma / dudvdw... using the Jacobian, computed in this function, as:
137 // as d^n sigma / dudvdw... = J * d^n sigma / dxdydz...
138 // where:
139 //
140 // dxdxdz... = J * dudvdw...
141 //
142 // | dx/du dx/dv dx/dw ... |
143 // | |
144 // d {x,y,z,...} | dy/du dy/dv dy/dw ... |
145 // J = ------------ = | |
146 // d {u,v,w,...} | dz/du dz/dv dz/dw ... |
147 // | |
148 // | ... ... ... ... |
149 //
150  SLOG("KineLimits", pDEBUG)
151  << "Computing Jacobian for transformation: "
152  << KinePhaseSpace::AsString(fromps) << " --> "
153  << KinePhaseSpace::AsString(tops);
154 
155  double J=0;
156  bool forward;
157  const Kinematics & kine = i->Kine();
158 
159  //
160  // transformation: {Q2}|E -> {lnQ2}|E
161  //
162  if ( TransformMatched(fromps,tops,kPSQ2fE,kPSlogQ2fE,forward) )
163  {
164  J = 1. / kine.Q2();
165  }
166 
167  //
168  // transformation: {QD2}|E -> {Q2}|E
169  //
170  else
171  if ( TransformMatched(fromps,tops,kPSQD2fE,kPSQ2fE,forward) )
172  {
173  J = TMath::Power(1+kine.Q2()/controls::kMQD2,-2)/controls::kMQD2;
174  }
175 
176  //
177  // transformation: {x,y}|E -> {lnx,lny}|E
178  //
179  else
180  if ( TransformMatched(fromps,tops,kPSxyfE,kPSlogxlogyfE,forward) )
181  {
182  J = 1. / (kine.x() * kine.y());
183  }
184 
185  //
186  // transformation: {Q2,y}|E -> {lnQ2,lny}|E
187  //
188  else
189  if ( TransformMatched(fromps,tops,kPSQ2yfE,kPSlogQ2logyfE,forward) )
190  {
191  J = 1. / (kine.Q2() * kine.y());
192  }
193 
194  //
195  // transformation: {Q2,y}|E -> {x,y}|E
196  //
197  else
198  if ( TransformMatched(fromps,tops,kPSQ2yfE,kPSxyfE,forward) )
199  {
200  const InitialState & init_state = i->InitState();
201  double Ev = init_state.ProbeE(kRfHitNucRest);
202  double M = init_state.Tgt().HitNucP4Ptr()->M();
203  double y = kine.y();
204  J = 2*y*Ev*M;
205  }
206 
207  //
208  // transformation: {W,Q2}|E -> {W,lnQ2}|E
209  //
210  else if ( TransformMatched(fromps,tops,kPSWQ2fE,kPSWlogQ2fE,forward) )
211  {
212  J = 1. / kine.Q2();
213  }
214 
215  //
216  // transformation: {W,QD2}|E -> {W,Q2}|E
217  //
218  else
219  if ( TransformMatched(fromps,tops,kPSWQD2fE,kPSWQ2fE,forward) )
220  {
221  J = TMath::Power(1+kine.Q2()/controls::kMQD2,-2)/controls::kMQD2;
222  }
223 
224  //
225  // transformation: {W2,Q2}|E --> {x,y}|E
226  //
227  else
228  if ( TransformMatched(fromps,tops,kPSW2Q2fE,kPSxyfE,forward) )
229  {
230  const InitialState & init_state = i->InitState();
231  double Ev = init_state.ProbeE(kRfHitNucRest);
232  double M = init_state.Tgt().HitNucP4Ptr()->M();
233  double y = kine.y();
234  J = TMath::Power(2*M*Ev,2) * y;
235  }
236 
237  //
238  // transformation: {W,Q2}|E -> {x,y}|E
239  //
240  else
241  if ( TransformMatched(fromps,tops,kPSWQ2fE,kPSxyfE,forward) )
242  {
243  const InitialState & init_state = i->InitState();
244  double Ev = init_state.ProbeE(kRfHitNucRest);
245  double M = init_state.Tgt().HitNucP4Ptr()->M();
246  double y = kine.y();
247  double W = kine.W();
248  J = 2*TMath::Power(M*Ev,2) * y/W;
249  }
250 
251  // Transformation: {Omegalep,Omegapi}|E -> {Omegalep,Thetapi}|E
252  else if ( TransformMatched(fromps,tops,kPSElOlOpifE,kPSElOlTpifE,forward) ) {
253  // Use symmetry to turn 4d integral into 3d * 2pi
254  J = 2*constants::kPi;
255  }
256 
257  else {
258  SLOG("KineLimits", pFATAL)
259  << "*** Can not compute Jacobian for transforming: "
260  << KinePhaseSpace::AsString(fromps) << " --> "
261  << KinePhaseSpace::AsString(tops);
262  exit(1);
263  }
264 
265  // if any of the above transforms was reverse, invert the jacobian
266  if(!forward) J = 1./J;
267 
268  return J;
269 }
270 //____________________________________________________________________________
273  KinePhaseSpace_t a, KinePhaseSpace_t b, bool & fwd)
274 {
275 
276  // match a->b or b->a
277  bool matched = ( (a==inpa&&b==inpb) || (a==inpb&&b==inpa) );
278 
279  if(matched) {
280  if(a==inpa&&b==inpb) fwd = true; // matched forward transform: a->b
281  else fwd = false; // matched reverse transform: b->a
282  }
283  return matched;
284 }
285 //____________________________________________________________________________
286 // Kinematical Limits:
287 //____________________________________________________________________________
288 Range1D_t genie::utils::kinematics::InelWLim(double Ev, double M, double ml)
289 {
290 // Computes W limits for inelastic v interactions
291 //
292  double M2 = TMath::Power(M,2);
293  double s = M2 + 2*M*Ev;
294  assert (s>0);
295 
296  Range1D_t W;
298  W.max = TMath::Sqrt(s) - ml;
299  if(W.max<=W.min) {
300  W.min = -1;
301  W.max = -1;
302  return W;
303  }
306  return W;
307 }
308 //____________________________________________________________________________
310  double Ev, double M, double ml, double W, double Q2min_cut)
311 {
312 // Computes Q2 limits (>0) @ the input W for inelastic v interactions
313 
314  Range1D_t Q2;
315  Q2.min = -1;
316  Q2.max = -1;
317 
318  double M2 = TMath::Power(M, 2.);
319  double ml2 = TMath::Power(ml, 2.);
320  double W2 = TMath::Power(W, 2.);
321  double s = M2 + 2*M*Ev;
322 
323  SLOG("KineLimits", pDEBUG) << "s = " << s;
324  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
325  assert (s>0);
326 
327  double auxC = 0.5*(s-M2)/s;
328  double aux1 = s + ml2 - W2;
329  double aux2 = aux1*aux1 - 4*s*ml2;
330 
331  (aux2 < 0) ? ( aux2 = 0 ) : ( aux2 = TMath::Sqrt(aux2) );
332 
333  Q2.max = -ml2 + auxC * (aux1 + aux2); // => 0
334  Q2.min = -ml2 + auxC * (aux1 - aux2); // => 0
335 
336  // guard against overflows
337  Q2.max = TMath::Max(0., Q2.max);
338  Q2.min = TMath::Max(0., Q2.min);
339 
340  // limit the minimum Q2
341  if(Q2.min < Q2min_cut) {Q2.min = Q2min_cut; }
342  if(Q2.max < Q2.min ) {Q2.min = -1; Q2.max = -1;}
343 
344  return Q2;
345 }
346 //____________________________________________________________________________
348  double Ev, double M, double ml, double W, double q2min_cut)
349 {
350 // Computes q2 (<0) limits @ the input W for inelastic v interactions
351 
352  Range1D_t Q2 = utils::kinematics::InelQ2Lim_W(Ev,M,ml,W,-1.*q2min_cut);
353  Range1D_t q2;
354  q2.min = - Q2.max;
355  q2.max = - Q2.min;
356  return q2;
357 }
358 //____________________________________________________________________________
360  double Ev, double M, double ml, double Q2min_cut)
361 {
362 // Computes Q2 (>0) limits irrespective of W for inelastic v interactions
363 
364  Range1D_t Q2;
365  Q2.min = -1;
366  Q2.max = -1;
367 
369  if(W.min<0) return Q2;
370 
371  Q2 = utils::kinematics::InelQ2Lim_W(Ev,M,ml,W.min,Q2min_cut);
372  return Q2;
373 }
374 //____________________________________________________________________________
376  double Ev, double M, double ml, double q2min_cut)
377 {
378 // Computes Q2 (>0) limits irrespective of W for inelastic v interactions
379 
380  Range1D_t Q2 = utils::kinematics::InelQ2Lim(Ev,M,ml,-1.*q2min_cut);
381  Range1D_t q2;
382  q2.min = - Q2.max;
383  q2.max = - Q2.min;
384  return q2;
385 }
386 //____________________________________________________________________________
387 Range1D_t genie::utils::kinematics::InelXLim(double Ev, double M, double ml)
388 
389 {
390 // Computes Bjorken x limits for inelastic v interactions
391 
392  double M2 = TMath::Power(M, 2.);
393  double ml2 = TMath::Power(ml,2.);
394  double s = M2 + 2*M*Ev;
395 
396  SLOG("KineLimits", pDEBUG) << "s = " << s;
397  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
398  assert (s>M2);
399 
400  Range1D_t x;
401  x.min = ml2/(s-M2) + controls::kASmallNum;
402  x.max = 1. - controls::kASmallNum;
403 
404  return x;
405 }
406 //____________________________________________________________________________
407 Range1D_t genie::utils::kinematics::InelYLim(double Ev, double M, double ml)
408 {
409 // Computes y limits for inelastic v interactions
410 
411  Range1D_t y;
412  y.min = 999;
413  y.max = -999;
414 
415  Range1D_t xl = kinematics::InelXLim(Ev,M,ml);
416  assert(xl.min>0 && xl.max>0);
417 
418  const unsigned int N=100;
419  const double logxmin = TMath::Log10(xl.min);
420  const double logxmax = TMath::Log10(xl.max);
421  const double dlogx = (logxmax-logxmin) / (double)(N-1);
422 
423  for(unsigned int i=0; i<N; i++) {
424  double x = TMath::Power(10, logxmin + i*dlogx);
425 
426  Range1D_t y_x = kinematics::InelYLim_X(Ev,M,ml,x);
427  if(y_x.min>=0 && y_x.min<=1) y.min = TMath::Min(y.min, y_x.min);
428  if(y_x.max>=0 && y_x.max<=1) y.max = TMath::Max(y.max, y_x.max);
429  }
430 
431  if(y.max >= 0 && y.max <= 1 && y.min >= 0 && y.min <= 1) {
432  y.min = TMath::Max(y.min, controls::kASmallNum);
433  y.max = TMath::Min(y.max, 1 - controls::kASmallNum);
434  } else {
435  y.min = -1;
436  y.max = -1;
437  }
438  SLOG("KineLimits", pDEBUG) << "y = [" << y.min << ", " << y.max << "]";
439  return y;
440 }
441 //____________________________________________________________________________
443  double Ev, double M, double ml, double x)
444 
445 {
446 // Computes y limits @ the input x for inelastic v interactions
447 
448  Range1D_t y;
449  y.min = -1;
450  y.max = -1;
451 
452  double Ev2 = TMath::Power(Ev,2);
453  double ml2 = TMath::Power(ml,2);
454 
455  SLOG("KineLimits", pDEBUG) << "x = " << x;
456  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
457 
458  assert (Ev>0);
459  assert (x>0&&x<1);
460 
461  double a = 0.5 * ml2/(M*Ev*x);
462  double b = ml2/Ev2;
463  double c = 1 + 0.5*x*M/Ev;
464  double d = TMath::Max(0., TMath::Power(1-a,2.) - b);
465 
466  double A = 0.5 * (1-a-0.5*b)/c;
467  double B = 0.5 * TMath::Sqrt(d)/c;
468 
469  y.min = TMath::Max(0., A-B) + controls::kASmallNum;
470  y.max = TMath::Min(1., A+B) - controls::kASmallNum;
471 
472  return y;
473 }
474 //____________________________________________________________________________
475 // Kinematical Limits for em interactions:
476 //____________________________________________________________________________
478 {
479 // Computes W limits for inelastic em interactions
480 //
481  double M2 = TMath::Power(M,2);
482  double ml2 = TMath::Power(ml,2); // added lepton mass squared to be used in s calculation
483  double s = M2 + 2*M*El + ml2; // non-negligible mass for em interactions
484  assert (s>0);
485 
486  Range1D_t W;
488  W.max = TMath::Sqrt(s) - ml;
489  if(W.max<=W.min) {
490  W.min = -1;
491  W.max = -1;
492  return W;
493  }
496  return W;
497 }
498 //____________________________________________________________________________
500  double El, double ml, double M, double W)
501 {
502 // Computes Q2 limits (>0) @ the input W for inelastic em interactions
503 
504  Range1D_t Q2;
505  Q2.min = -1;
506  Q2.max = -1;
507 
508  double M2 = TMath::Power(M, 2.);
509  double ml2 = TMath::Power(ml, 2.);
510  double W2 = TMath::Power(W, 2.);
511  double s = M2 + 2*M*El + ml2; // added lepton mass squared to be used in s calculation (non-negligible mass for em interactions)
512 
513  SLOG("KineLimits", pDEBUG) << "s = " << s;
514  SLOG("KineLimits", pDEBUG) << "El = " << El;
515  assert (s>0);
516 
517  double auxC = 0.5*(s - M2 - ml2)/s; // subtract ml2 to account for the non-negligible mass of the incoming lepton
518  double aux1 = s + ml2 - W2;
519  double aux2 = aux1*aux1 - 4*s*ml2;
520 
521  (aux2 < 0) ? ( aux2 = 0 ) : ( aux2 = TMath::Sqrt(aux2) );
522 
523  Q2.max = -ml2 + auxC * (aux1 + aux2); // => 0
524  Q2.min = -ml2 + auxC * (aux1 - aux2); // => 0
525 
526  // guard against overflows
527  Q2.max = TMath::Max(0., Q2.max);
528  Q2.min = TMath::Max(0., Q2.min);
529 
530  // limit the minimum Q2
531  if(Q2.min < utils::kinematics::electromagnetic::kMinQ2Limit) {Q2.min = utils::kinematics::electromagnetic::kMinQ2Limit; } // use the relevant threshold for em scattering
532  if(Q2.max < Q2.min ) {Q2.min = -1; Q2.max = -1;}
533 
534  return Q2;
535 }
536 //____________________________________________________________________________
538  double El, double ml, double M, double W)
539 {
540 // Computes q2 (<0) limits @ the input W for inelastic em interactions
541 
543  Range1D_t q2;
544  q2.min = - Q2.max;
545  q2.max = - Q2.min;
546  return q2;
547 }
548 //____________________________________________________________________________
550  double El, double ml, double M)
551 {
552 // Computes Q2 (>0) limits irrespective of W for inelastic em interactions
553 
554  Range1D_t Q2;
555  Q2.min = -1;
556  Q2.max = -1;
557 
559  if(W.min<0) return Q2;
560 
562  return Q2;
563 }
564 //____________________________________________________________________________
566  double El, double ml, double M)
567 {
568 // Computes Q2 (>0) limits irrespective of W for inelastic em interactions
569 
571  Range1D_t q2;
572  q2.min = - Q2.max;
573  q2.max = - Q2.min;
574  return q2;
575 }
576 //____________________________________________________________________________
578 
579 {
580 // Computes Bjorken x limits for inelastic em interactions
581 
582  double M2 = TMath::Power(M, 2.);
583  double ml2 = TMath::Power(ml,2.);
584  double s = M2 + 2*M*El + ml2; // added lepton mass squared to be used in s calculation (non-negligible mass for em interactions)
585 
586  SLOG("KineLimits", pDEBUG) << "s = " << s;
587  SLOG("KineLimits", pDEBUG) << "El = " << El;
588  assert (s > M2 + ml2); // added lepton mass squared to be used in s calculation (non-negligible mass for em interactions)
589 
590  Range1D_t x;
591  x.min = ml2/(s - M2 - ml2) + controls::kASmallNum; // subtracted lepton mass squared to be used in s calculation (non-negligible mass for em interactions)
592  x.max = 1. - controls::kASmallNum;
593 
594  return x;
595 }
596 //____________________________________________________________________________
598 {
599 // Computes y limits for inelastic v interactions
600 
601  Range1D_t y;
602  y.min = 999;
603  y.max = -999;
604 
606  assert(xl.min>0 && xl.max>0);
607 
608  const unsigned int N=100;
609  const double logxmin = TMath::Log10(xl.min);
610  const double logxmax = TMath::Log10(xl.max);
611  const double dlogx = (logxmax-logxmin) / (double)(N-1);
612 
613  for(unsigned int i=0; i<N; i++) {
614  double x = TMath::Power(10, logxmin + i*dlogx);
615 
617  if(y_x.min>=0 && y_x.min<=1) y.min = TMath::Min(y.min, y_x.min);
618  if(y_x.max>=0 && y_x.max<=1) y.max = TMath::Max(y.max, y_x.max);
619  }
620 
621  if(y.max >= 0 && y.max <= 1 && y.min >= 0 && y.min <= 1) {
622  y.min = TMath::Max(y.min, controls::kASmallNum);
623  y.max = TMath::Min(y.max, 1 - controls::kASmallNum);
624  } else {
625  y.min = -1;
626  y.max = -1;
627  }
628  SLOG("KineLimits", pDEBUG) << "y = [" << y.min << ", " << y.max << "]";
629  return y;
630 }
631 //____________________________________________________________________________
633  double El, double ml, double M, double x)
634 
635 {
636 // Computes y limits @ the input x for inelastic em interactions
637 
638  Range1D_t y;
639  y.min = -1;
640  y.max = -1;
641 
642  double El2 = TMath::Power(El,2);
643  double ml2 = TMath::Power(ml,2);
644 
645  SLOG("KineLimits", pDEBUG) << "x = " << x;
646  SLOG("KineLimits", pDEBUG) << "El = " << El;
647 
648  assert (El>0);
649  assert (x>0&&x<1);
650 
651  double a = 0.5 * ml2/(M*El*x);
652  double b = ml2/El2;
653  double c = 1 + 0.5*x*M/El;
654  double d = TMath::Max(0., TMath::Power(1-a,2.) - b);
655 
656  double A = 0.5 * (1-a-0.5*b)/c;
657  double B = 0.5 * TMath::Sqrt(d)/c;
658 
659  y.min = TMath::Max(0., A-B) + controls::kASmallNum;
660  y.max = TMath::Min(1., A+B) - controls::kASmallNum;
661 
662  return y;
663 }
664 //____________________________________________________________________________
666 {
667 // Computes x limits for coherent v interactions
668 
670  return x;
671 }
672 //____________________________________________________________________________
673 Range1D_t genie::utils::kinematics::CohQ2Lim(double Mn, double mpi, double mlep, double Ev)
674 {
675  // The expressions for Q^2 min appears in PRD 74, 054007 (2006) by
676  // Kartavtsev, Paschos, and Gounaris
677 
678  Range1D_t Q2;
679  Q2.min = 0.0;
680  Q2.max = std::numeric_limits<double>::max(); // Value must be overriden in user options
681 
682  double Mn2 = Mn * Mn;
683  double mlep2 = mlep * mlep;
684  double s = Mn2 + 2.0 * Mn * Ev;
685  double W2min = CohW2Min(Mn, mpi);
686 
687  // Looks like Q2min = A * B - C, where A, B, and C are complicated
688  double a = 1.0;
689  double b = mlep2 / s;
690  double c = W2min / s;
691  double lambda = a * a + b * b + c * c - 2.0 * a * b - 2.0 * a * c - 2.0 * b * c;
692  if (lambda > 0) {
693  double A = (s - Mn * Mn) / 2.0;
694  double B = 1 - TMath::Sqrt(lambda);
695  double C = 0.5 * (W2min + mlep2 - Mn2 * (W2min - mlep2) / s );
696  if (A * B - C < 0) {
697  SLOG("KineLimits", pERROR)
698  << "Q2 kinematic limits calculation failed for CohQ2Lim. "
699  << "Assuming Q2min = 0.0";
700  }
701  Q2.min = TMath::Max(0., A * B - C);
702  } else {
703  SLOG("KineLimits", pERROR)
704  << "Q2 kinematic limits calculation failed for CohQ2Lim. "
705  << "Assuming Q2min = 0.0";
706  }
707 
708  return Q2;
709 }
710 //____________________________________________________________________________
711 Range1D_t genie::utils::kinematics::Cohq2Lim(double Mn, double mpi, double mlep, double Ev)
712 {
713  Range1D_t Q2 = utils::kinematics::CohQ2Lim(Mn, mpi, mlep, Ev);
714  Range1D_t q2;
715  q2.min = - Q2.max;
716  q2.max = - Q2.min;
717  return q2;
718 }
719 //____________________________________________________________________________
720 Range1D_t genie::utils::kinematics::CohW2Lim(double Mn, double mpi, double mlep,
721  double Ev, double Q2)
722 {
723  // These expressions for W^2 min and max appear in PRD 74, 054007 (2006) by
724  // Kartavtsev, Paschos, and Gounaris
725 
726  Range1D_t W2l;
727  W2l.min = -1;
728  W2l.max = -1;
729 
730  double s = Mn * Mn + 2.0 * Mn * Ev;
731  double Mnterm = 1 - Mn * Mn / s;
732  double Mlterm = 1 - mlep * mlep / s;
733  // Here T1, T2 are generically "term 1" and "term 2" in a long expression
734  double T1 = 0.25 * s * s * Mnterm * Mnterm * Mlterm;
735  double T2 = Q2 - (0.5 * s * Mnterm) + (0.5 * mlep * mlep * Mnterm);
736 
737  W2l.min = CohW2Min(Mn, mpi);
738  W2l.max = (T1 - T2 * T2 ) *
739  (1.0 / Mnterm) *
740  (1.0 / (Q2 + mlep * mlep));
741 
742  return W2l;
743 }
744 //____________________________________________________________________________
746  double Q2, double Mn, double xsi)
747 {
748  Range1D_t nul;
749  nul.min = -1;
750  nul.max = -1;
751 
752  double nu_min = (W2min + Q2 - Mn * Mn) / (2.0 * Mn);
753  double nu_max = (W2max + Q2 - Mn * Mn) / (2.0 * Mn);
754  double xsiQ = xsi * TMath::Sqrt(Q2);
755 
756  nul.min = (xsiQ > nu_min) ? xsiQ : nu_min;
757  nul.max = nu_max;
758 
759  return nul;
760 }
761 //____________________________________________________________________________
762 Range1D_t genie::utils::kinematics::CohYLim(double Mn, double mpi, double mlep,
763  double Ev, double Q2, double xsi)
764 {
765  Range1D_t ylim;
766  ylim.min = -1;
767  ylim.max = -1;
768 
769  Range1D_t W2lim = genie::utils::kinematics::CohW2Lim(Mn, mpi, mlep, Ev, Q2);
770  if (W2lim.min > W2lim.max) {
771  LOG("KineLimits", pDEBUG)
772  << "Kinematically forbidden region in CohYLim. W2min = " << W2lim.min
773  << "; W2max =" << W2lim.max;
774  LOG("KineLimits", pDEBUG)
775  << " Mn = " << Mn << "; mpi = " << mpi << "; mlep = "
776  << mlep << "; Ev = " << Ev << "; Q2 = " << Q2;
777  return ylim;
778  }
780  Q2, Mn, xsi);
781  ylim.min = nulim.min / Ev;
782  ylim.max = nulim.max / Ev;
783 
784  return ylim;
785 }
786 //____________________________________________________________________________
788 {
789 // Computes y limits for coherent v interactions
790 
792  1.-ml/EvL - controls::kASmallNum);
793  return y;
794 }
795 //____________________________________________________________________________
796 // Helpers for kinematic limits
797 //____________________________________________________________________________
798 double genie::utils::kinematics::CohW2Min(double Mn, double mpi)
799 {
800  // These expressions for W^2 min and max appear in PRD 74, 054007 (2006) by
801  // Kartavtsev, Paschos, and Gounaris
802 
803  return (Mn + mpi) * (Mn + mpi);
804 }
805 //____________________________________________________________________________
806 Range1D_t genie::utils::kinematics::DarkWLim(double Ev, double M, double ml)
807 {
808 // Computes W limits for inelastic v interactions
809 //
810  double M2 = TMath::Power(M,2);
811  double ml2 = TMath::Power(ml,2);
812  double s = M2 + 2*M*Ev + ml2;
813  assert (s>0);
814 
815  Range1D_t W;
817 // W.min = kNeutronMass + kPionMass;
818  W.max = TMath::Sqrt(s) - ml;
819  if(W.max<=W.min) {
820  W.min = -1;
821  W.max = -1;
822  return W;
823  }
826  return W;
827 }
828 //____________________________________________________________________________
830  double Ev, double M, double ml, double W, double Q2min_cut)
831 {
832 // Computes Q2 limits (>0) @ the input W for inelastic v interactions
833 
834  Range1D_t Q2;
835  Q2.min = -1;
836  Q2.max = -1;
837 
838  double M2 = TMath::Power(M, 2.);
839  double ml2 = TMath::Power(ml, 2.);
840  double W2 = TMath::Power(W, 2.);
841  double s = M2 + 2*M*Ev + ml2;
842  assert(s > 0);
843  double sqs = TMath::Sqrt(s);
844  double E1CM = (s + ml2 - M2) / (2.*sqs);
845  double p1CM = TMath::Max(0., E1CM*E1CM - ml2);
846  p1CM = TMath::Sqrt(p1CM);
847  double E3CM = (s + ml2 - W2) / (2.*sqs);
848  double p3CM = TMath::Max(0., E3CM*E3CM - ml2);
849  p3CM = TMath::Sqrt(p3CM);
850 
851  SLOG("KineLimits", pDEBUG) << "s = " << s;
852  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
853  SLOG("KineLimits", pDEBUG) << "M = " << M;
854  SLOG("KineLimits", pDEBUG) << "W = " << W;
855  SLOG("KineLimits", pDEBUG) << "E1_CM = " << E1CM;
856  SLOG("KineLimits", pDEBUG) << "p1_CM = " << p1CM;
857  SLOG("KineLimits", pDEBUG) << "E3_CM = " << E3CM;
858  SLOG("KineLimits", pDEBUG) << "p3_CM = " << p3CM;
859 
860  Q2.min = TMath::Power(p3CM - p1CM,2) - TMath::Power((W2 - M2) / (2.*sqs),2);
861  Q2.max = TMath::Power(p3CM + p1CM,2) - TMath::Power((W2 - M2) / (2.*sqs),2);
862 
863  SLOG("KineLimits", pDEBUG) << "Nominal Q^2 limits: " << Q2.min << " , " << Q2.max;
864  // guard against overflows
865  Q2.max = TMath::Max(0., Q2.max);
866  Q2.min = TMath::Max(0., Q2.min);
867 
868  // limit the minimum Q2
869  if(Q2.min < Q2min_cut) {Q2.min = Q2min_cut; }
870  if(Q2.max < Q2.min ) {Q2.min = -1; Q2.max = -1;}
871 
872  return Q2;
873 }
874 //____________________________________________________________________________
876  double Ev, double M, double ml, double W, double q2min_cut)
877 {
878 // Computes q2 (<0) limits @ the input W for inelastic v interactions
879 
880  Range1D_t Q2 = utils::kinematics::DarkQ2Lim_W(Ev,M,ml,W,-1.*q2min_cut);
881  Range1D_t q2;
882  q2.min = - Q2.max;
883  q2.max = - Q2.min;
884  return q2;
885 }
886 //____________________________________________________________________________
888  double Ev, double M, double ml, double Q2min_cut)
889 {
890 // Computes Q2 (>0) limits irrespective of W for inelastic v interactions
891 
892  Range1D_t Q2;
893  Q2.min = -1;
894  Q2.max = -1;
895 
897  if(W.min<0) return Q2;
898 
899  Q2 = utils::kinematics::DarkQ2Lim_W(Ev,M,ml,W.min,Q2min_cut);
900  return Q2;
901 }
902 //____________________________________________________________________________
904  double Ev, double M, double ml, double q2min_cut)
905 {
906 // Computes Q2 (>0) limits irrespective of W for inelastic v interactions
907 
908  Range1D_t Q2 = utils::kinematics::DarkQ2Lim(Ev,M,ml,-1.*q2min_cut);
909  Range1D_t q2;
910  q2.min = - Q2.max;
911  q2.max = - Q2.min;
912  return q2;
913 }
914 //____________________________________________________________________________
915 Range1D_t genie::utils::kinematics::DarkXLim(double Ev, double M, double ml)
916 
917 {
918 // Computes Bjorken x limits for inelastic interactions
919 // For the dark matter case, it is relatively straightforward
920  Range1D_t Wl = utils::kinematics::DarkWLim(Ev, M, ml);
921  double Wmin = Wl.min;
922  double W2min = Wmin*Wmin;
923  SLOG("KineLimits", pDEBUG) << "W^2_min = " << W2min;
924  Range1D_t Q2l = utils::kinematics::DarkQ2Lim_W(Ev, M, ml, Wmin);
925  SLOG("KineLimits", pDEBUG) << "Q^2 range : " << Q2l.min << " , " << Q2l.max;
926  double M2 = M*M;
927  Range1D_t x;
928  x.min = Q2l.min / (Q2l.min + W2min - M2);
929  x.max = Q2l.max / (Q2l.max + W2min - M2);
930 
931  SLOG("KineLimits", pDEBUG) << "x = [" << x.min << ", " << x.max << "]";
932  return x;
933 }
934 //____________________________________________________________________________
935 Range1D_t genie::utils::kinematics::DarkYLim(double Ev, double M, double ml)
936 {
937  // For dark inelastic scattering, can compute exactly and is much simpler
938  Range1D_t Wl = utils::kinematics::DarkWLim(Ev, M, ml);
939  double Wmin = Wl.min;
940  double W2min = Wmin*Wmin;
941  Range1D_t Q2l = utils::kinematics::DarkQ2Lim_W(Ev, M, ml, Wmin);
942  double M2 = M*M;
943  Range1D_t y;
944  y.min = (Q2l.min + W2min - M2) / (2*Ev*M);
945  y.max = (Q2l.max + W2min - M2) / (2*Ev*M);
946 
947  SLOG("KineLimits", pDEBUG) << "y = [" << y.min << ", " << y.max << "]";
948  return y;
949 }
950 //____________________________________________________________________________
952  double Ev, double M, double ml, double x)
953 
954 {
955  // Computes y limits @ the input x for inelastic interactions
956  // We hit y_min when W = W_min and y_max when Q^2 = Q_min^2 or Q^2 = Q_max^2
957 
958  Range1D_t y;
959  y.min = -1;
960  y.max = -1;
961 
962  Range1D_t Wl = utils::kinematics::DarkWLim(Ev, M, ml);
963  double Wmin = Wl.min;
964  double W2min = Wmin*Wmin;
965  double M2 = M*M;
966  double ml2 = ml*ml;
967  y.min = (W2min - M2) / (1.-x) / (2*Ev*M);
968  y.max = 2.* M * x *(Ev*Ev - ml2) / Ev / (2. * M * Ev * x + M2 * x * x + ml2);
969 
970  return y;
971 }
972 //____________________________________________________________________________
973 // Kinematical Transformations:
974 //____________________________________________________________________________
976 {
977 // Q2 -> QD2 transformation. QD2 takes out the dipole form of the form factors
978 // making the differential cross section to be flatter and speeding up the
979 // kinematical selection.
980 
981  assert(Q2>0);
982  return TMath::Power(1+Q2/controls::kMQD2, -1);
983 }
984 //____________________________________________________________________________
986 {
987  assert(QD2>0);
988  return controls::kMQD2*(1/QD2-1);
989 }
990 //____________________________________________________________________________
991 double genie::utils::kinematics::Q2(const Interaction * const interaction)
992 {
993 // Get Q^2 from kinematics object
994 // If Q^2 is not set but x,y are, then compute Q2 from x,y
995 
996  const Kinematics & kinematics = interaction->Kine();
997 
998  if (kinematics.KVSet(kKVQ2) || kinematics.KVSet(kKVq2)) {
999  double Q2 = kinematics.Q2();
1000  return Q2;
1001  }
1002  if (kinematics.KVSet(kKVy)) {
1003  const InitialState & init_state = interaction->InitState();
1004  double Mn = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
1005  double x = kinematics.x();
1006  double y = kinematics.y();
1007  double Ev = init_state.ProbeE(kRfHitNucRest);
1008  double Q2 = 2*Mn*Ev*x*y;
1009  return Q2;
1010  }
1011  SLOG("KineLimits", pERROR) << "Couldn't compute Q^2 for \n"<< *interaction;
1012  return 0;
1013 }
1014 //____________________________________________________________________________
1015 double genie::utils::kinematics::W(const Interaction * const interaction)
1016 {
1017  const ProcessInfo & process_info = interaction->ProcInfo();
1018 
1019  if(process_info.IsQuasiElastic()) {
1020  // hadronic inv. mass is equal to the recoil nucleon on-shell mass
1021  int rpdgc = interaction->RecoilNucleonPdg();
1022  double M = PDGLibrary::Instance()->Find(rpdgc)->Mass();
1023  return M;
1024  }
1025 
1026  const Kinematics & kinematics = interaction->Kine();
1027  if(kinematics.KVSet(kKVW)) {
1028  double W = kinematics.W();
1029  return W;
1030  }
1031  if(kinematics.KVSet(kKVx) && kinematics.KVSet(kKVy)) {
1032  const InitialState & init_state = interaction->InitState();
1033  double Ev = init_state.ProbeE(kRfHitNucRest);
1034  double M = init_state.Tgt().HitNucP4Ptr()->M();
1035  double M2 = M*M;
1036  double x = kinematics.x();
1037  double y = kinematics.y();
1038  double W2 = TMath::Max(0., M2 + 2*Ev*M*y*(1-x));
1039  double W = TMath::Sqrt(W2);
1040  return W;
1041  }
1042  SLOG("KineLimits", pERROR) << "Couldn't compute W for \n"<< *interaction;
1043  return 0;
1044 }
1045 //___________________________________________________________________________
1047  double Ev, double M, double W, double Q2, double & x, double & y)
1048 {
1049 // Converts (W,Q2) => (x,y)
1050 // Uses the system: a) W^2 - M^2 = 2*Ev*M*y*(1-x) and b) Q^2 = 2*x*y*M*Ev
1051 // Ev is the neutrino energy at the struck nucleon rest frame
1052 // M is the nucleon mass - it does not need to be on the mass shell
1053 
1054  double M2 = TMath::Power(M,2);
1055  double W2 = TMath::Power(W,2);
1056 
1057  x = Q2 / (W2-M2+Q2);
1058  y = (W2-M2+Q2) / (2*M*Ev);
1059 
1060  x = TMath::Min(1.,x);
1061  y = TMath::Min(1.,y);
1062  x = TMath::Max(0.,x);
1063  y = TMath::Max(0.,y);
1064 
1065  LOG("KineLimits", pDEBUG)
1066  << "(W=" << W << ",Q2=" << Q2 << ") => (x="<< x << ", y=" << y<< ")";
1067 }
1068 //___________________________________________________________________________
1070  double Ev, double M, double & W, double & Q2, double x, double y)
1071 {
1072 // Converts (x,y) => (W,Q2)
1073 // Uses the system: a) W^2 - M^2 = 2*Ev*M*y*(1-x) and b) Q^2 = 2*x*y*M*Ev
1074 // Ev is the neutrino energy at the struck nucleon rest frame
1075 // M is the nucleon mass - it does not need to be on the mass shell
1076 
1077  double M2 = TMath::Power(M,2);
1078  double W2 = M2 + 2*Ev*M*y*(1-x);
1079 
1080  W = TMath::Sqrt(TMath::Max(0., W2));
1081  Q2 = 2*x*y*M*Ev;
1082 
1083  LOG("KineLimits", pDEBUG)
1084  << "(x=" << x << ",y=" << y << " => (W=" << W << ",Q2=" << Q2 << ")";
1085 }
1086 //___________________________________________________________________________
1088  double Ev, double M, double x, double y)
1089 {
1090 // Converts (x,y) => W
1091 // Ev is the neutrino energy at the struck nucleon rest frame
1092 // M is the nucleon mass - it does not need to be on the mass shell
1093 
1094  double M2 = TMath::Power(M,2);
1095  double W2 = M2 + 2*Ev*M*y*(1-x);
1096  double W = TMath::Sqrt(TMath::Max(0., W2));
1097 
1098  LOG("KineLimits", pDEBUG) << "(x=" << x << ",y=" << y << ") => W=" << W;
1099 
1100  return W;
1101 }
1102 //___________________________________________________________________________
1104  double Ev, double M, double x, double y)
1105 {
1106 // Converts (x,y) => Q2
1107 // Ev is the neutrino energy at the struck nucleon rest frame
1108 // M is the nucleon mass - it does not need to be on the mass shell
1109 
1110  double Q2 = 2*x*y*M*Ev;
1111 
1112  LOG("KineLimits", pDEBUG) << "(x=" << x << ",y=" << y << ") => Q2=" << Q2;
1113 
1114  return Q2;
1115 }
1116 //___________________________________________________________________________
1118  double Ev, double M, double Q2, double y)
1119 {
1120 // Converts (Q^2,y) => x
1121 // Ev is the neutrino energy at the struck nucleon rest frame
1122 // M is the nucleon mass - it does not need to be on the mass shell
1123  assert(Ev > 0. && M > 0. && Q2 > 0. && y > 0.);
1124 
1125  double x = Q2 / (2. * y * M * Ev);
1126 
1127  LOG("KineLimits", pDEBUG) << "(Ev=" << Ev << ",Q2=" << Q2
1128  << ",y=" << y << ",M=" << M << ") => x=" << x;
1129 
1130  return x;
1131 }
1132 //___________________________________________________________________________
1134  double x, double Q2, double M, double mc)
1135 {
1136 // x : scaling variable (plain or modified)
1137 // Q2: momentum transfer
1138 // M : hit nucleon "mass" (nucleon can be off the mass shell)
1139 // mc: charm mass
1140 //
1141  double M2 = TMath::Power(M,2);
1142  double v = 0.5*Q2/(M*x);
1143  double W2 = TMath::Max(0., M2+2*M*v-Q2);
1144  double W = TMath::Sqrt(W2);
1145  double Wmin = M + kLightestChmHad;
1146  double xc = utils::kinematics::SlowRescalingVar(x,Q2,M,mc);
1147 
1148  if(xc>=1 || W<=Wmin) return false;
1149  else return true;
1150 }
1151 //____________________________________________________________________________
1153  double x, double Q2, double M, double mc)
1154 {
1155 // x : scaling variable (plain or modified)
1156 // Q2: momentum transfer
1157 // M : hit nucleon "mass" (nucleon can be off the mass shell)
1158 // mc: charm mass
1159 
1160  double mc2 = TMath::Power(mc,2);
1161  double v = 0.5*Q2/(M*x);
1162  double xc = x + 0.5*mc2/(M*v);
1163  return xc;
1164 }
1165 //____________________________________________________________________________
1167  Range1D_t & range, double min_cut, double max_cut)
1168 {
1169  // if the min,max are within the existing limits, the cut can be applied
1170  // by narrowing down the xisting limits
1171  if ( utils::math::IsWithinLimits(min_cut, range ) ) range.min = min_cut;
1172  if ( utils::math::IsWithinLimits(max_cut, range ) ) range.max = max_cut;
1173 
1174  // if the min-cut is above the existing max-limit or
1175  // if the max-cut is below the existing min-limit then
1176  // the range should be invalidated
1177 
1178  if (min_cut > range.max || max_cut < range.min) {
1179 
1180  range.min = 0;
1181  range.max = 0;
1182  }
1183 }
1184 //___________________________________________________________________________
1186 {
1187  Kinematics * kine = in->KinePtr();
1188 
1189  if(kine->KVSet(kKVx) && kine->KVSet(kKVy)) {
1190  const InitialState & init_state = in->InitState();
1191  double Ev = init_state.ProbeE(kRfHitNucRest);
1192  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off mass shell
1193  double x = kine->x();
1194  double y = kine->y();
1195 
1196  double Q2=-1,W=-1;
1197  kinematics::XYtoWQ2(Ev,M,W,Q2,x,y);
1198  kine->SetQ2(Q2);
1199  kine->SetW(W);
1200  }
1201 }
1202 //___________________________________________________________________________
1204 {
1205  Kinematics * kine = in->KinePtr();
1206 
1207  if(kine->KVSet(kKVW) && kine->KVSet(kKVQ2)) {
1208  const InitialState & init_state = in->InitState();
1209  double Ev = init_state.ProbeE(kRfHitNucRest);
1210  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off mass shell
1211  double W = kine->W();
1212  double Q2 = kine->Q2();
1213 
1214  double x=-1,y=-1;
1215  kinematics::WQ2toXY(Ev,M,W,Q2,x,y);
1216  kine->Setx(x);
1217  kine->Sety(y);
1218  }
1219 }
1220 //___________________________________________________________________________
1222 {
1223  Kinematics * kine = in->KinePtr();
1224 
1225  if(kine->KVSet(kKVy) && kine->KVSet(kKVQ2)) {
1226 
1227  const ProcessInfo & pi = in->ProcInfo();
1228  const InitialState & init_state = in->InitState();
1229  double M = 0.0;
1230  double Ev = 0.0;
1231 
1232  if (pi.IsCoherent()) {
1233  M = in->InitState().Tgt().Mass(); // nucleus mass
1234  Ev = init_state.ProbeE(kRfLab);
1235  }
1236  else {
1237  M = in->InitState().Tgt().HitNucP4Ptr()->M(); //can be off m/shell
1238  Ev = init_state.ProbeE(kRfHitNucRest);
1239  }
1240 
1241  double y = kine->y();
1242  double Q2 = kine->Q2();
1243 
1244  double x = kinematics::Q2YtoX(Ev,M,Q2,y);
1245  kine->Setx(x);
1246  }
1247 }
1248 //____________________________________________________________________________
1250  double * x, double * par)
1251 {
1252  //-- inputs
1253  double QD2 = x[0]; // QD2 (Q2 transformed to take out the dipole form)
1254  double W = x[1]; // invariant mass
1255 
1256  //-- parameters
1257  double mD = par[0]; // resonance mass
1258  double gD = par[1]; // resonance width
1259  double xsmax = par[2]; // safety factor * max cross section in (W,Q2)
1260  double Wmax = par[3]; // kinematically allowed Wmax
1261  //double E = par[4]; // neutrino energy
1262 
1263 // return 3*xsmax; ////////////////NOTE
1264 
1265  xsmax*=5;
1266 
1267  double func = 0;
1268 
1269  if(Wmax > mD) {
1270  // -- if the resonance mass is within the kinematical range then
1271  // -- the envelope is defined as a plateau above the resonance peak,
1272  // -- a steeply falling leading edge (low-W side) and a more slowly
1273  // -- falling trailing edge (high-W side)
1274 
1275  double hwfe = mD+2*gD; // high W falling edge
1276  double lwfe = mD-gD/2; // low W falling edge
1277 
1278  if(W < lwfe) {
1279  //low-W falling edge
1280  func = xsmax / (1 + 5* TMath::Power((W-lwfe)/gD,2));
1281  } else if (W > hwfe) {
1282  //high-W falling edge
1283  func = xsmax / (1 + TMath::Power((W-hwfe)/gD,2));
1284  } else {
1285  // plateau
1286  func = xsmax;
1287  }
1288  } else {
1289  // -- if the resonance mass is above the kinematical range then the
1290  // -- envelope is a small plateau just bellow Wmax and a falling edge
1291  // -- at lower W
1292 
1293  double plateau_edge = Wmax-0.1;
1294 
1295  if (W > plateau_edge) {
1296  // plateau
1297  func = xsmax;
1298  } else {
1299  //low-W falling edge
1300  func = xsmax / (1 + TMath::Power((W-plateau_edge)/gD,2));
1301  }
1302  }
1303 
1304  if(QD2<0.5) {
1305  func *= (1 - (0.5-QD2));
1306  }
1307 
1308  return func;
1309 }
1310 //___________________________________________________________________________
1312  double * x, double * par)
1313 {
1314  //-- inputs
1315  double xb = x[0]; // scaling variable x brorken
1316  //double y = x[1]; // inelasticity y
1317 
1318  //-- parameters
1319  double xpeak = par[0]; // x peak
1320  double xmin = par[1]; // min x
1321  double xmax = 1.; // max x
1322  double xsmax = par[2]; // safety factor * max cross section in (x,y)
1323 
1324  double func = 0;
1325 
1326  if(xb < xpeak/20.) {
1327  //low-x falling edge
1328  double plateau_edge = xpeak/20.;
1329  double slope = xsmax/(xmin - plateau_edge);
1330  func = xsmax - slope * (xb - plateau_edge);
1331  } else if (xb > 2*xpeak) {
1332  //high-x falling edge
1333  double plateau_edge = 2*xpeak;
1334  double slope = xsmax/(xmax - plateau_edge);
1335  func = xsmax - slope * (xb - plateau_edge);
1336  } else {
1337  // plateau
1338  func = xsmax;
1339  }
1340  return func;
1341 }
1342 //___________________________________________________________________________
1344  double * x, double * par)
1345 {
1346  //-- inputs
1347  double xb = x[0]; // x
1348  double yb = x[1]; // y
1349 
1350  //-- parameters
1351  double xsmax = 3*par[0]; // safety factor * max cross section in (x,y)
1352  double Ev = par[1]; // neutrino energy;
1353 
1354  if(yb<0.|| yb>1.) return 0.;
1355  if(xb<0.|| xb>1.) return 0.;
1356 
1357  if(Ev<1) return xsmax;
1358  if(xb/Ev<1E-4 && yb>0.95) return 5*xsmax;
1359 
1360  double func = 0;
1361  double xp = 0.1;
1362  double yp = (Ev>2.5) ? 2.5/Ev : 1;
1363 
1364  if(xb>xp) {
1365  double xs0=0;
1366  if(yb<yp) {
1367  xs0 = xsmax;
1368  } else {
1369  xs0 = xsmax - (yb-yp)*xsmax;
1370  }
1371  double d = TMath::Power( (xb-xp)/0.075, 2);
1372  func = xs0/(1 + d);
1373  } else {
1374  if(yb>yp) {
1375  func = xsmax - (yb-yp)*xsmax;
1376  } else {
1377  func = xsmax;
1378  }
1379  }
1380  return func;
1381 }
1382 //___________________________________________________________________________
double COHImportanceSamplingEnvelope(double *x, double *par)
Definition: KineUtils.cxx:1343
static const double kMaxX
Definition: Controls.h:44
double W(bool selected=false) const
Definition: Kinematics.cxx:167
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
bool TransformMatched(KinePhaseSpace_t ia, KinePhaseSpace_t ib, KinePhaseSpace_t a, KinePhaseSpace_t b, bool &fwd)
Definition: KineUtils.cxx:271
std::map< std::string, double > xmax
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Range1D_t InelXLim(double El, double ml, double M)
Definition: KineUtils.cxx:577
Range1D_t InelWLim(double El, double ml, double M)
Definition: KineUtils.cxx:477
#define pERROR
Definition: Messenger.h:60
void UpdateXYFromWQ2(const Interaction *in)
Definition: KineUtils.cxx:1203
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
Range1D_t Inelq2Lim_W(double El, double ml, double M, double W)
Definition: KineUtils.cxx:537
double DISImportanceSamplingEnvelope(double *x, double *par)
Definition: KineUtils.cxx:1311
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:265
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Range1D_t InelWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:288
int RecoilNucleonPdg(void) const
recoil nucleon pdg
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Range1D_t Darkq2Lim(double Ev, double M, double ml, double q2min_cut=-1 *controls::kMinQ2Limit)
Definition: KineUtils.cxx:903
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
static const double kMaxY
Definition: Controls.h:46
#define pFATAL
Definition: Messenger.h:57
Range1D_t DarkQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:829
bool IsWithinLimits(double x, Range1D_t range)
Definition: MathUtils.cxx:265
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
static const double kPhotontest
Definition: Constants.h:80
Int_t par
Definition: SimpleIterate.C:24
Range1D_t CohYLim(double Mn, double mpi, double mlep, double Ev, double Q2, double xsi)
Definition: KineUtils.cxx:762
double Wmax[kNWBins]
double x(bool selected=false) const
Definition: Kinematics.cxx:109
static const double kLightestChmHad
Definition: Constants.h:79
double CohW2Min(double Mn, double mpi)
Definition: KineUtils.cxx:798
static const double kMinY
Definition: Controls.h:45
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
Range1D_t InelQ2Lim(double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:359
enum genie::EKinePhaseSpace KinePhaseSpace_t
double Mass(void) const
Definition: Target.cxx:241
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
double XYtoW(double Ev, double M, double x, double y)
Definition: KineUtils.cxx:1087
double PhaseSpaceVolume(const Interaction *const i, KinePhaseSpace_t ps)
Definition: KineUtils.cxx:34
Range1D_t InelQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:309
double y(bool selected=false) const
Definition: Kinematics.cxx:122
double W(const Interaction *const i)
Definition: KineUtils.cxx:1015
Double_t q2[12][num]
Definition: f2_nu.C:137
const double C
const XML_Char * s
Definition: expat.h:262
Summary information for an interaction.
Definition: Interaction.h:56
double SlowRescalingVar(double x, double Q2, double M, double mc)
Definition: KineUtils.cxx:1152
Range1D_t InelXLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:387
Range1D_t DarkQ2Lim(double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:887
Range1D_t Inelq2Lim_W(double Ev, double M, double ml, double W, double q2min_cut=-1 *controls::kMinQ2Limit)
Definition: KineUtils.cxx:347
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Range1D_t CohNuLim(double W2min, double W2max, double Q2, double Mn, double xsi)
Definition: KineUtils.cxx:745
bool KVSet(KineVar_t kv) const
Definition: Kinematics.cxx:327
double QD2toQ2(double QD2)
Definition: KineUtils.cxx:985
Range1D_t Darkq2Lim_W(double Ev, double M, double ml, double W, double q2min_cut=-1 *controls::kMinQ2Limit)
Definition: KineUtils.cxx:875
static string AsString(KinePhaseSpace_t kps)
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
const double a
double XYtoQ2(double Ev, double M, double x, double y)
Definition: KineUtils.cxx:1103
void XYtoWQ2(double Ev, double M, double &W, double &Q2, double x, double y)
Definition: KineUtils.cxx:1069
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1046
const Kinematics & Kine(void) const
Definition: Interaction.h:71
Kinematical phase space.
Definition: KPhaseSpace.h:34
static const double kNeutronMass
Definition: Constants.h:77
Range1D_t Cohq2Lim(double Mn, double mpi, double mlep, double Ev)
Definition: KineUtils.cxx:711
Range1D_t CohQ2Lim(double Mn, double mpi, double mlep, double Ev)
Definition: KineUtils.cxx:673
Float_t d
Definition: plot.C:236
static const double kMinX
Definition: Controls.h:43
double func(double x, double y)
Range1D_t CohXLim(void)
Definition: KineUtils.cxx:665
static const double kASmallNum
Definition: Controls.h:40
bool IsAboveCharmThreshold(double x, double Q2, double M, double mc)
Definition: KineUtils.cxx:1133
double Q2toQD2(double Q2)
Definition: KineUtils.cxx:975
Range1D_t InelQ2Lim(double El, double ml, double M)
Definition: KineUtils.cxx:549
Range1D_t Inelq2Lim(double El, double ml, double M)
Definition: KineUtils.cxx:565
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:241
static const double kMQD2
Definition: Controls.h:65
void UpdateWQ2FromXY(const Interaction *in)
Definition: KineUtils.cxx:1185
Range1D_t InelYLim(double El, double ml, double M)
Definition: KineUtils.cxx:597
static const double A
Definition: Units.h:82
double max
Definition: Range1.h:54
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
ifstream in
Definition: comparison.C:7
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
Range1D_t InelQ2Lim_W(double El, double ml, double M, double W)
Definition: KineUtils.cxx:499
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
static const double kPionMass
Definition: Constants.h:74
Range1D_t InelYLim_X(double Ev, double M, double ml, double x)
Definition: KineUtils.cxx:442
double Q2YtoX(double Ev, double M, double Q2, double y)
Definition: KineUtils.cxx:1117
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
Range1D_t DarkYLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:935
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:253
void UpdateXFromQ2Y(const Interaction *in)
Definition: KineUtils.cxx:1221
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
Range1D_t InelYLim_X(double El, double ml, double M, double x)
Definition: KineUtils.cxx:632
exit(0)
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double Wmin[kNWBins]
double min
Definition: Range1.h:53
Range1D_t DarkYLim_X(double Ev, double M, double ml, double x)
Definition: KineUtils.cxx:951
Range1D_t Inelq2Lim(double Ev, double M, double ml, double q2min_cut=-1 *controls::kMinQ2Limit)
Definition: KineUtils.cxx:375
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
const Target & Tgt(void) const
Definition: InitialState.h:67
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
void kinematics()
Definition: kinematics.C:10
#define W(x)
Range1D_t InelYLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:407
double ProbeE(RefFrame_t rf) const
static const double kPi
Definition: Constants.h:38
Range1D_t DarkXLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:915
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Range1D_t CohW2Lim(double Mn, double mpi, double mlep, double Ev, double Q2)
Definition: KineUtils.cxx:720
void ApplyCutsToKineLimits(Range1D_t &r, double min, double max)
Definition: KineUtils.cxx:1166
Range1D_t DarkWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:806
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
bool IsCoherent(void) const
Definition: ProcessInfo.cxx:97
double RESImportanceSamplingEnvelope(double *x, double *par)
Definition: KineUtils.cxx:1249