BWFunc.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  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab - November 22, 2004
9 
10  For documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ May 01, 2016 - Libo Jiang
14  Add W dependence to Delta->N gamma
15 
16 */
17 //____________________________________________________________________________
18 
19 #include <cassert>
20 
21 #include <TMath.h>
22 
23 #include "Framework/Utils/BWFunc.h"
25 
26 using namespace genie;
27 using namespace genie::constants;
28 
29 //
31  double W, int L, double mass, double width0, double norm)
32 {
33 //Inputs:
34 // - W: Invariant mass (GeV)
35 // - L: Resonance orbital angular momentum
36 // - mass: Resonance mass (GeV)
37 // - width0: Resonance width
38 // - norm: Breit Wigner norm
39 
40  //-- sanity checks
41  assert(mass > 0);
42  assert(width0 > 0);
43  assert(norm > 0);
44  assert(W > 0);
45  assert(L >= 0);
46 
47  //-- auxiliary parameters
48  double mN = kNucleonMass;
49  double mPi = kPi0Mass;
50  double m_2 = TMath::Power(mass, 2);
51  double mN_2 = TMath::Power(mN, 2);
52  double W_2 = TMath::Power(W, 2);
53 
54  double m=mass;
55  //m_aux1 m_aux2
56  double m_aux1= TMath::Power(mN+mPi, 2);
57  double m_aux2= TMath::Power(mN-mPi, 2);
58 
59  double BRPi0 = 0.994; //Npi Branching Ratio
60  double BRgamma0 = 0.006; //Ngamma Branching Ratio
61 
62  double widPi0 = width0*BRPi0;
63  double widgamma0= width0*BRgamma0;
64 
65  //-- calculate the L-dependent resonance width
66  double EgammaW= (W_2-mN_2)/(2*W);
67  double Egammam= (m_2-mN_2)/(2*m);
68 
69 
70  if(EgammaW<0) {
71 // cout<< "Two small W!!! W is lower than one Nucleon Mass!!!!"<<endl;
72  return 0;
73  }
74  //pPiW pion momentum
75  double pPiW = 0;
76  //
77  if(W_2>m_aux1) pPiW = TMath::Sqrt((W_2-m_aux1)*(W_2-m_aux2))/(2*W);
78 
79  double pPim = TMath::Sqrt((m_2-m_aux1)*(m_2-m_aux2))/(2*m);
80 
81  //double TPiW = pPiW*TMath::Power(pPiW*rDelta, 2*L)/(1+TMath::Power(pPiW*rDelta, 2));
82  //double TPim = pPim*TMath::Power(pPim*rDelta, 2*L)/(1+TMath::Power(pPim*rDelta, 2));
83 
84  // Form factors
85  //double fgammaW= 1/(TMath::Power(1+EgammaW*EgammaW/0.706, 2)*(1+EgammaW*EgammaW/3.519));
86  //double fgammam= 1/(TMath::Power(1+Egammam*Egammam/0.706, 2)*(1+Egammam*Egammam/3.519));
87  double fgammaW= 1/(TMath::Power(1+EgammaW*EgammaW/0.706, 2));
88  double fgammam= 1/(TMath::Power(1+Egammam*Egammam/0.706, 2));
89 
90 
91  double EgammaW_3=TMath::Power(EgammaW, 3);
92  double Egammam_3=TMath::Power(Egammam, 3);
93  double fgammaW_2=TMath::Power(fgammaW, 2);
94  double fgammam_2=TMath::Power(fgammam, 2);
95 
96  //double width = widPi0*(TPiW/TPim)+widgamma0*(EgammaW_3*fgammaW_2/(Egammam_3*fgammam_2));
97  double width = widPi0*TMath::Power((pPiW/pPim),3)+widgamma0*(EgammaW_3*fgammaW_2/(Egammam_3*fgammam_2));
98  //-- calculate the Breit Wigner function for the input W
99  double width_2 = TMath::Power( width, 2);
100  double W_m_2 = TMath::Power( W-mass, 2);
101 
102  double bw = (0.5/kPi) * (width/norm) / (W_m_2 + 0.25*width_2);
103 
104  return bw;
105 }
106 //______________________________________________________________________
108  double W, int L, double mass, double width0, double norm)
109 {
110 //Inputs:
111 // - W: Invariant mass (GeV)
112 // - L: Resonance orbital angular momentum
113 // - mass: Resonance mass (GeV)
114 // - width0: Resonance width
115 // - norm: Breit Wigner norm
116 
117  //-- sanity checks
118  assert(mass > 0);
119  assert(width0 > 0);
120  assert(norm > 0);
121  assert(W > 0);
122  assert(L >= 0);
123 
124  //-- auxiliary parameters
125  double mN = kNucleonMass;
126  double mPi = kPi0Mass;
127  double m_2 = TMath::Power(mass, 2);
128  double mN_2 = TMath::Power(mN, 2);
129  double mPi_2 = TMath::Power(mPi, 2);
130  double W_2 = TMath::Power(W, 2);
131 
132  //-- calculate the L-dependent resonance width
133  double qpW_2 = ( TMath::Power(W_2 - mN_2 - mPi_2, 2) - 4*mN_2*mPi_2 );
134  double qpM_2 = ( TMath::Power(m_2 - mN_2 - mPi_2, 2) - 4*mN_2*mPi_2 );
135  if(qpW_2 < 0) qpW_2 = 0;
136  if(qpM_2 < 0) qpM_2 = 0;
137  double qpW = TMath::Sqrt(qpW_2) / (2*W);
138  double qpM = TMath::Sqrt(qpM_2) / (2*mass);
139  double width = width0 * TMath::Power( qpW/qpM, 2*L+1 );
140 
141  //-- calculate the Breit Wigner function for the input W
142  double width_2 = TMath::Power( width, 2);
143  double W_m_2 = TMath::Power( W-mass, 2);
144 
145  double bw = (0.5/kPi) * (width/norm) / (W_m_2 + 0.25*width_2);
146  return bw;
147 }
148 //______________________________________________________________________
150  double W, double mass, double width, double norm)
151 {
152 //Inputs:
153 // - W: Invariant mass (GeV)
154 // - mass: Resonance mass (GeV)
155 // - width: Resonance width
156 // - norm: Breit Wigner norm
157 
158  //-- sanity checks
159  assert(mass > 0);
160  assert(width > 0);
161  assert(norm > 0);
162  assert(W > 0);
163 
164  //-- auxiliary parameters
165  double width_2 = TMath::Power( width, 2);
166  double W_m_2 = TMath::Power( W-mass, 2);
167 
168  //-- calculate the Breit Wigner function for the input W
169  double bw = (0.5/kPi) * (width/norm) / (W_m_2 + 0.25*width_2);
170  return bw;
171 }
172 //______________________________________________________________________
173 
double BreitWigner(double W, double mass, double width, double norm)
Definition: BWFunc.cxx:149
Basic constants.
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
static const double kNucleonMass
Definition: Constants.h:78
static const double kPi0Mass
Definition: Constants.h:75
double BreitWignerLGamma(double W, int L, double mass, double width0, double norm)
Definition: BWFunc.cxx:30
double BreitWignerL(double W, int L, double mass, double width0, double norm)
Definition: BWFunc.cxx:107
static constexpr double L
Float_t norm
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
assert(nhit_max >=nhit_nbins)
#define W(x)
static const double kPi
Definition: Constants.h:38