BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtEvalHelAmp Class Reference

#include <EvtEvalHelAmp.hh>

Public Member Functions

 EvtEvalHelAmp (EvtSpinType::spintype A, EvtSpinType::spintype B, EvtSpinType::spintype C, EvtComplexPtrPtr HBC)
virtual ~EvtEvalHelAmp ()
double probMax ()
void evalAmp (EvtParticle *p, EvtAmp &amp)

Detailed Description

Definition at line 33 of file EvtEvalHelAmp.hh.

Constructor & Destructor Documentation

◆ EvtEvalHelAmp()

EvtEvalHelAmp::EvtEvalHelAmp ( EvtSpinType::spintype A,
EvtSpinType::spintype B,
EvtSpinType::spintype C,
EvtComplexPtrPtr HBC )

Definition at line 78 of file EvtEvalHelAmp.cc.

79 {
80
81 // find out how many states each particle have
82 _nA = EvtSpinType::getSpinStates( typeA );
83 _nB = EvtSpinType::getSpinStates( typeB );
84 _nC = EvtSpinType::getSpinStates( typeC );
85
86 // find out what 2 times the spin is
87 _JA2 = EvtSpinType::getSpin2( typeA );
88 _JB2 = EvtSpinType::getSpin2( typeB );
89 _JC2 = EvtSpinType::getSpin2( typeC );
90
91 // allocate memory
92 _lambdaA2 = new int[_nA];
93 _lambdaB2 = new int[_nB];
94 _lambdaC2 = new int[_nC];
95
96 _HBC = new EvtComplexPtr[_nB];
97 int ia, ib, ic;
98 for ( ib = 0; ib < _nB; ib++ ) { _HBC[ib] = new EvtComplex[_nC]; }
99
100 _RA = new EvtComplexPtr[_nA];
101 for ( ia = 0; ia < _nA; ia++ ) { _RA[ia] = new EvtComplex[_nA]; }
102 _RB = new EvtComplexPtr[_nB];
103 for ( ib = 0; ib < _nB; ib++ ) { _RB[ib] = new EvtComplex[_nB]; }
104 _RC = new EvtComplexPtr[_nC];
105 for ( ic = 0; ic < _nC; ic++ ) { _RC[ic] = new EvtComplex[_nC]; }
106
107 _amp = new EvtComplexPtrPtr[_nA];
108 _amp1 = new EvtComplexPtrPtr[_nA];
109 _amp3 = new EvtComplexPtrPtr[_nA];
110 for ( ia = 0; ia < _nA; ia++ )
111 {
112 _amp[ia] = new EvtComplexPtr[_nB];
113 _amp1[ia] = new EvtComplexPtr[_nB];
114 _amp3[ia] = new EvtComplexPtr[_nB];
115 for ( ib = 0; ib < _nB; ib++ )
116 {
117 _amp[ia][ib] = new EvtComplex[_nC];
118 _amp1[ia][ib] = new EvtComplex[_nC];
119 _amp3[ia][ib] = new EvtComplex[_nC];
120 }
121 }
122
123 // find the allowed helicities (actually 2*times the helicity!)
124
125 fillHelicity( _lambdaA2, _nA, _JA2 );
126 fillHelicity( _lambdaB2, _nB, _JB2 );
127 fillHelicity( _lambdaC2, _nC, _JC2 );
128
129 for ( ib = 0; ib < _nB; ib++ )
130 {
131 for ( ic = 0; ic < _nC; ic++ ) { _HBC[ib][ic] = HBC[ib][ic]; }
132 }
133}
EvtComplexPtr
EvtComplex * EvtComplexPtr
Definition EvtComplex.hh:68
EvtComplexPtrPtr
EvtComplexPtr * EvtComplexPtrPtr
Definition EvtComplex.hh:69
EvtSpinType::getSpin2
static int getSpin2(spintype stype)
Definition EvtSpinType.hh:45
EvtSpinType::getSpinStates
static int getSpinStates(spintype stype)
Definition EvtSpinType.hh:67

◆ ~EvtEvalHelAmp()

EvtEvalHelAmp::~EvtEvalHelAmp ( )
virtual

Definition at line 39 of file EvtEvalHelAmp.cc.

39 {
40
41 // allocate memory
42 delete[] _lambdaA2;
43 delete[] _lambdaB2;
44 delete[] _lambdaC2;
45
46 int ia, ib, ic;
47 for ( ib = 0; ib < _nB; ib++ ) { delete[] _HBC[ib]; }
48
49 delete[] _HBC;
50
51 for ( ia = 0; ia < _nA; ia++ ) { delete[] _RA[ia]; }
52 delete[] _RA;
53
54 for ( ib = 0; ib < _nB; ib++ ) { delete[] _RB[ib]; }
55 delete[] _RB;
56
57 for ( ic = 0; ic < _nC; ic++ ) { delete[] _RC[ic]; }
58 delete[] _RC;
59
60 for ( ia = 0; ia < _nA; ia++ )
61 {
62 for ( ib = 0; ib < _nB; ib++ )
63 {
64 delete[] _amp[ia][ib];
65 delete[] _amp1[ia][ib];
66 delete[] _amp3[ia][ib];
67 }
68 delete[] _amp[ia];
69 delete[] _amp1[ia];
70 delete[] _amp3[ia];
71 }
72
73 delete[] _amp;
74 delete[] _amp1;
75 delete[] _amp3;
76}

Member Function Documentation

◆ evalAmp()

void EvtEvalHelAmp::evalAmp ( EvtParticle * p,
EvtAmp & amp )

Definition at line 176 of file EvtEvalHelAmp.cc.

176 {
177
178 // find theta and phi of the first daughter
179
180 EvtVector4R pB = p->getDaug( 0 )->getP4();
181 EvtVector4R pC = p->getDaug( 1 )->getP4();
182 EvtVector4R pP = pB + pC;
183
184 EvtHelSys angles( pP, pB );
185 double theta = angles.getHelAng( 1 );
186 double phi = angles.getHelAng( 2 );
187 // double theta=acos(pB.get(3)/pB.d3mag()); //pingrg ,2007,4,15
188 // double phi=atan2(pB.get(2),pB.get(1));
189
190 double c = sqrt( ( _JA2 + 1 ) / ( 4 * EvtConst::pi ) );
191
192 int ia, ib, ic;
193
194 double prob1 = 0.0;
195
196 for ( ia = 0; ia < _nA; ia++ )
197 {
198 for ( ib = 0; ib < _nB; ib++ )
199 {
200 for ( ic = 0; ic < _nC; ic++ )
201 {
202 _amp[ia][ib][ic] = 0.0;
203 if ( abs( _lambdaB2[ib] - _lambdaC2[ic] ) <= _JA2 )
204 {
205 double dfun =
206 EvtdFunction::d( _JA2, _lambdaA2[ia], _lambdaB2[ib] - _lambdaC2[ic], theta );
207
208 _amp[ia][ib][ic] =
209 c * _HBC[ib][ic] *
210 exp( EvtComplex( 0.0, phi * 0.5 *
211 ( _lambdaA2[ia] - _lambdaB2[ib] + _lambdaC2[ic] ) ) ) *
212 dfun;
213 }
214 prob1 += real( _amp[ia][ib][ic] * conj( _amp[ia][ib][ic] ) );
215 }
216 }
217 }
218
219 setUpRotationMatrices( p, theta, phi );
220
221 applyRotationMatrices();
222
223 double prob2 = 0.0;
224
225 for ( ia = 0; ia < _nA; ia++ )
226 {
227 for ( ib = 0; ib < _nB; ib++ )
228 {
229 for ( ic = 0; ic < _nC; ic++ )
230 {
231 prob2 += real( _amp[ia][ib][ic] * conj( _amp[ia][ib][ic] ) );
232 if ( _nA == 1 )
233 {
234 if ( _nB == 1 )
235 {
236 if ( _nC == 1 ) { amp.vertex( _amp[ia][ib][ic] ); }
237 else { amp.vertex( ic, _amp[ia][ib][ic] ); }
238 }
239 else
240 {
241 if ( _nC == 1 ) { amp.vertex( ib, _amp[ia][ib][ic] ); }
242 else
243 {
244 amp.vertex( ib, ic, _amp[ia][ib][ic] );
245 // report(INFO,"EvtGen") << "ib,ic:"<<ib<<" "<<ic<<" "<<_amp[ia][ib][ic]<<endl;
246 }
247 }
248 }
249 else
250 {
251 if ( _nB == 1 )
252 {
253 if ( _nC == 1 ) { amp.vertex( ia, _amp[ia][ib][ic] ); }
254 else { amp.vertex( ia, ic, _amp[ia][ib][ic] ); }
255 }
256 else
257 {
258 if ( _nC == 1 ) { amp.vertex( ia, ib, _amp[ia][ib][ic] ); }
259 else { amp.vertex( ia, ib, ic, _amp[ia][ib][ic] ); }
260 }
261 }
262 }
263 }
264 }
265
266 if ( fabs( prob1 - prob2 ) > 0.000001 * prob1 )
267 {
268 report( INFO, "EvtGen" ) << "prob1,prob2:" << prob1 << " " << prob2 << endl;
269 ::abort();
270 }
271
272 return;
273}
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition Evt3Rank3C.cc:166
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
report
ostream & report(Severity severity, const char *facility)
Definition EvtReport.cc:34
INFO
@ INFO
Definition EvtReport.hh:52
EvtAmp::vertex
void vertex(const EvtComplex &amp)
Definition EvtAmp.cc:441
EvtConst::pi
static const double pi
Definition EvtConst.hh:27
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition EvtParticle.cc:108
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:80
EvtdFunction::d
static double d(int j, int m1, int m2, double theta)
Definition EvtdFunction.cc:29

Referenced by EvtKstarstargamma::decay().

◆ probMax()

double EvtEvalHelAmp::probMax ( )

Definition at line 135 of file EvtEvalHelAmp.cc.

135 {
136
137 double c = 1.0 / sqrt( 4 * EvtConst::pi / ( _JA2 + 1 ) );
138
139 int ia, ib, ic;
140
141 double theta;
142 int itheta;
143
144 double maxprob = 0.0;
145
146 for ( itheta = -10; itheta <= 10; itheta++ )
147 {
148 theta = acos( 0.099999 * itheta );
149 for ( ia = 0; ia < _nA; ia++ )
150 {
151 double prob = 0.0;
152 for ( ib = 0; ib < _nB; ib++ )
153 {
154 for ( ic = 0; ic < _nC; ic++ )
155 {
156 _amp[ia][ib][ic] = 0.0;
157 if ( abs( _lambdaB2[ib] - _lambdaC2[ic] ) <= _JA2 )
158 {
159 _amp[ia][ib][ic] =
160 c * _HBC[ib][ic] *
161 EvtdFunction::d( _JA2, _lambdaA2[ia], _lambdaB2[ib] - _lambdaC2[ic], theta );
162 prob += real( _amp[ia][ib][ic] * conj( _amp[ia][ib][ic] ) );
163 }
164 }
165 }
166
167 prob *= sqrt( 1.0 * _nA );
168
169 if ( prob > maxprob ) maxprob = prob;
170 }
171 }
172
173 return maxprob;
174}
The documentation for this class was generated from the following files: