]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | //-------------------------------------------------------------------------- |
2 | // | |
3 | // Environment: | |
4 | // This software is part of the EvtGen package developed jointly | |
5 | // for the BaBar and CLEO collaborations. If you use all or part | |
6 | // of it, please give an appropriate acknowledgement. | |
7 | // | |
8 | // Copyright Information: See EvtGen/COPYRIGHT | |
9 | // Copyright (C) 2000 Caltech, UCSB | |
10 | // | |
11 | // Module: EvtHelAmp.cc | |
12 | // | |
13 | // Description: Decay model for implementation of generic 2 body | |
14 | // decay specified by the partial wave amplitudes | |
15 | // | |
16 | // | |
17 | // Modification history: | |
18 | // | |
19 | // fkw February 2, 2001 changes to satisfy KCC | |
20 | // RYD September 7, 2000 Module created | |
21 | // | |
22 | //------------------------------------------------------------------------ | |
23 | // | |
24 | #include "EvtGenBase/EvtPatches.hh" | |
25 | #include <stdlib.h> | |
26 | #include "EvtGenBase/EvtParticle.hh" | |
27 | #include "EvtGenBase/EvtPDL.hh" | |
28 | #include "EvtGenBase/EvtVector4C.hh" | |
29 | #include "EvtGenBase/EvtReport.hh" | |
30 | #include "EvtGenBase/EvtEvalHelAmp.hh" | |
31 | #include "EvtGenBase/EvtId.hh" | |
32 | #include "EvtGenBase/EvtConst.hh" | |
33 | #include "EvtGenBase/EvtdFunction.hh" | |
34 | #include "EvtGenBase/EvtAmp.hh" | |
35 | using std::endl; | |
36 | ||
37 | ||
38 | EvtEvalHelAmp::~EvtEvalHelAmp() { | |
39 | ||
40 | //deallocate memory | |
41 | delete [] _lambdaA2; | |
42 | delete [] _lambdaB2; | |
43 | delete [] _lambdaC2; | |
44 | ||
45 | int ia,ib,ic; | |
46 | for(ib=0;ib<_nB;ib++){ | |
47 | delete [] _HBC[ib]; | |
48 | } | |
49 | ||
50 | delete [] _HBC; | |
51 | ||
52 | ||
53 | for(ia=0;ia<_nA;ia++){ | |
54 | delete [] _RA[ia]; | |
55 | } | |
56 | delete [] _RA; | |
57 | ||
58 | for(ib=0;ib<_nB;ib++){ | |
59 | delete [] _RB[ib]; | |
60 | } | |
61 | delete [] _RB; | |
62 | ||
63 | for(ic=0;ic<_nC;ic++){ | |
64 | delete [] _RC[ic]; | |
65 | } | |
66 | delete [] _RC; | |
67 | ||
68 | ||
69 | for(ia=0;ia<_nA;ia++){ | |
70 | for(ib=0;ib<_nB;ib++){ | |
71 | delete [] _amp[ia][ib]; | |
72 | delete [] _amp1[ia][ib]; | |
73 | delete [] _amp3[ia][ib]; | |
74 | } | |
75 | delete [] _amp[ia]; | |
76 | delete [] _amp1[ia]; | |
77 | delete [] _amp3[ia]; | |
78 | } | |
79 | ||
80 | delete [] _amp; | |
81 | delete [] _amp1; | |
82 | delete [] _amp3; | |
83 | ||
84 | } | |
85 | ||
86 | ||
87 | EvtEvalHelAmp::EvtEvalHelAmp(EvtId idA, | |
88 | EvtId idB, | |
89 | EvtId idC, | |
90 | EvtComplexPtrPtr HBC){ | |
91 | ||
92 | ||
93 | EvtSpinType::spintype typeA=EvtPDL::getSpinType(idA); | |
94 | EvtSpinType::spintype typeB=EvtPDL::getSpinType(idB); | |
95 | EvtSpinType::spintype typeC=EvtPDL::getSpinType(idC); | |
96 | ||
97 | //find out how many states each particle have | |
98 | _nA=EvtSpinType::getSpinStates(typeA); | |
99 | _nB=EvtSpinType::getSpinStates(typeB); | |
100 | _nC=EvtSpinType::getSpinStates(typeC); | |
101 | ||
102 | //find out what 2 times the spin is | |
103 | _JA2=EvtSpinType::getSpin2(typeA); | |
104 | _JB2=EvtSpinType::getSpin2(typeB); | |
105 | _JC2=EvtSpinType::getSpin2(typeC); | |
106 | ||
107 | ||
108 | //allocate memory | |
109 | _lambdaA2=new int[_nA]; | |
110 | _lambdaB2=new int[_nB]; | |
111 | _lambdaC2=new int[_nC]; | |
112 | ||
113 | _HBC=new EvtComplexPtr[_nB]; | |
114 | int ia,ib,ic; | |
115 | for(ib=0;ib<_nB;ib++){ | |
116 | _HBC[ib]=new EvtComplex[_nC]; | |
117 | } | |
118 | ||
119 | ||
120 | _RA=new EvtComplexPtr[_nA]; | |
121 | for(ia=0;ia<_nA;ia++){ | |
122 | _RA[ia]=new EvtComplex[_nA]; | |
123 | } | |
124 | _RB=new EvtComplexPtr[_nB]; | |
125 | for(ib=0;ib<_nB;ib++){ | |
126 | _RB[ib]=new EvtComplex[_nB]; | |
127 | } | |
128 | _RC=new EvtComplexPtr[_nC]; | |
129 | for(ic=0;ic<_nC;ic++){ | |
130 | _RC[ic]=new EvtComplex[_nC]; | |
131 | } | |
132 | ||
133 | _amp=new EvtComplexPtrPtr[_nA]; | |
134 | _amp1=new EvtComplexPtrPtr[_nA]; | |
135 | _amp3=new EvtComplexPtrPtr[_nA]; | |
136 | for(ia=0;ia<_nA;ia++){ | |
137 | _amp[ia]=new EvtComplexPtr[_nB]; | |
138 | _amp1[ia]=new EvtComplexPtr[_nB]; | |
139 | _amp3[ia]=new EvtComplexPtr[_nB]; | |
140 | for(ib=0;ib<_nB;ib++){ | |
141 | _amp[ia][ib]=new EvtComplex[_nC]; | |
142 | _amp1[ia][ib]=new EvtComplex[_nC]; | |
143 | _amp3[ia][ib]=new EvtComplex[_nC]; | |
144 | } | |
145 | } | |
146 | ||
147 | //find the allowed helicities (actually 2*times the helicity!) | |
148 | ||
149 | fillHelicity(_lambdaA2,_nA,_JA2,idA); | |
150 | fillHelicity(_lambdaB2,_nB,_JB2,idB); | |
151 | fillHelicity(_lambdaC2,_nC,_JC2,idC); | |
152 | ||
153 | for(ib=0;ib<_nB;ib++){ | |
154 | for(ic=0;ic<_nC;ic++){ | |
155 | _HBC[ib][ic]=HBC[ib][ic]; | |
156 | } | |
157 | } | |
158 | } | |
159 | ||
160 | ||
161 | ||
162 | ||
163 | ||
164 | ||
165 | double EvtEvalHelAmp::probMax(){ | |
166 | ||
167 | double c=1.0/sqrt(4*EvtConst::pi/(_JA2+1)); | |
168 | ||
169 | int ia,ib,ic; | |
170 | ||
171 | ||
172 | double theta; | |
173 | int itheta; | |
174 | ||
175 | double maxprob=0.0; | |
176 | ||
177 | for(itheta=-10;itheta<=10;itheta++){ | |
178 | theta=acos(0.099999*itheta); | |
179 | for(ia=0;ia<_nA;ia++){ | |
180 | double prob=0.0; | |
181 | for(ib=0;ib<_nB;ib++){ | |
182 | for(ic=0;ic<_nC;ic++){ | |
183 | _amp[ia][ib][ic]=0.0; | |
184 | if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) { | |
185 | _amp[ia][ib][ic]=c*_HBC[ib][ic]* | |
186 | EvtdFunction::d(_JA2,_lambdaA2[ia], | |
187 | _lambdaB2[ib]-_lambdaC2[ic],theta); | |
188 | prob+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic])); | |
189 | } | |
190 | } | |
191 | } | |
192 | ||
193 | prob*=sqrt(1.0*_nA); | |
194 | ||
195 | if (prob>maxprob) maxprob=prob; | |
196 | ||
197 | } | |
198 | } | |
199 | ||
200 | return maxprob; | |
201 | ||
202 | } | |
203 | ||
204 | ||
205 | void EvtEvalHelAmp::evalAmp( EvtParticle *p, EvtAmp& amp){ | |
206 | ||
207 | //find theta and phi of the first daughter | |
208 | ||
209 | EvtVector4R pB=p->getDaug(0)->getP4(); | |
210 | ||
211 | double theta=acos(pB.get(3)/pB.d3mag()); | |
212 | double phi=atan2(pB.get(2),pB.get(1)); | |
213 | ||
214 | double c=sqrt((_JA2+1)/(4*EvtConst::pi)); | |
215 | ||
216 | int ia,ib,ic; | |
217 | ||
218 | double prob1=0.0; | |
219 | ||
220 | for(ia=0;ia<_nA;ia++){ | |
221 | for(ib=0;ib<_nB;ib++){ | |
222 | for(ic=0;ic<_nC;ic++){ | |
223 | _amp[ia][ib][ic]=0.0; | |
224 | if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) { | |
225 | double dfun=EvtdFunction::d(_JA2,_lambdaA2[ia], | |
226 | _lambdaB2[ib]-_lambdaC2[ic],theta); | |
227 | ||
228 | _amp[ia][ib][ic]=c*_HBC[ib][ic]* | |
229 | exp(EvtComplex(0.0,phi*0.5*(_lambdaA2[ia]-_lambdaB2[ib]+ | |
230 | _lambdaC2[ic])))*dfun; | |
231 | } | |
232 | prob1+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic])); | |
233 | } | |
234 | } | |
235 | } | |
236 | ||
237 | setUpRotationMatrices(p,theta,phi); | |
238 | ||
239 | applyRotationMatrices(); | |
240 | ||
241 | double prob2=0.0; | |
242 | ||
243 | for(ia=0;ia<_nA;ia++){ | |
244 | for(ib=0;ib<_nB;ib++){ | |
245 | for(ic=0;ic<_nC;ic++){ | |
246 | prob2+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic])); | |
247 | if (_nA==1){ | |
248 | if (_nB==1){ | |
249 | if (_nC==1){ | |
250 | amp.vertex(_amp[ia][ib][ic]); | |
251 | } | |
252 | else{ | |
253 | amp.vertex(ic,_amp[ia][ib][ic]); | |
254 | } | |
255 | } | |
256 | else{ | |
257 | if (_nC==1){ | |
258 | amp.vertex(ib,_amp[ia][ib][ic]); | |
259 | } | |
260 | else{ | |
261 | amp.vertex(ib,ic,_amp[ia][ib][ic]); | |
262 | } | |
263 | } | |
264 | }else{ | |
265 | if (_nB==1){ | |
266 | if (_nC==1){ | |
267 | amp.vertex(ia,_amp[ia][ib][ic]); | |
268 | } | |
269 | else{ | |
270 | amp.vertex(ia,ic,_amp[ia][ib][ic]); | |
271 | } | |
272 | } | |
273 | else{ | |
274 | if (_nC==1){ | |
275 | amp.vertex(ia,ib,_amp[ia][ib][ic]); | |
276 | } | |
277 | else{ | |
278 | amp.vertex(ia,ib,ic,_amp[ia][ib][ic]); | |
279 | } | |
280 | } | |
281 | } | |
282 | } | |
283 | } | |
284 | } | |
285 | ||
286 | if (fabs(prob1-prob2)>0.000001*prob1){ | |
287 | report(INFO,"EvtGen") << "prob1,prob2:"<<prob1<<" "<<prob2<<endl; | |
288 | ::abort(); | |
289 | } | |
290 | ||
291 | return ; | |
292 | ||
293 | } | |
294 | ||
295 | ||
296 | void EvtEvalHelAmp::fillHelicity(int* lambda2,int n,int J2, EvtId id){ | |
297 | ||
298 | int i; | |
299 | ||
300 | //photon is special case! | |
301 | if (n==2&&J2==2) { | |
302 | lambda2[0]=2; | |
303 | lambda2[1]=-2; | |
304 | return; | |
305 | } | |
306 | ||
307 | //and so is the neutrino! | |
308 | if (n==1&&J2==1) { | |
309 | if (EvtPDL::getStdHep(id)>0){ | |
310 | //particle i.e. lefthanded | |
311 | lambda2[0]=-1; | |
312 | }else{ | |
313 | //anti particle i.e. righthanded | |
314 | lambda2[0]=1; | |
315 | } | |
316 | return; | |
317 | } | |
318 | ||
319 | ||
320 | assert(n==J2+1); | |
321 | ||
322 | for(i=0;i<n;i++){ | |
323 | lambda2[i]=n-i*2-1; | |
324 | } | |
325 | ||
326 | return; | |
327 | ||
328 | } | |
329 | ||
330 | ||
331 | void EvtEvalHelAmp::setUpRotationMatrices(EvtParticle* p,double theta, double phi){ | |
332 | ||
333 | switch(_JA2){ | |
334 | ||
335 | case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8: | |
336 | ||
337 | { | |
338 | ||
339 | EvtSpinDensity R=p->rotateToHelicityBasis(); | |
340 | ||
341 | ||
342 | int i,j,n; | |
343 | ||
344 | n=R.getDim(); | |
345 | ||
346 | assert(n==_nA); | |
347 | ||
348 | ||
349 | for(i=0;i<n;i++){ | |
350 | for(j=0;j<n;j++){ | |
351 | _RA[i][j]=R.get(i,j); | |
352 | } | |
353 | } | |
354 | ||
355 | } | |
356 | ||
357 | break; | |
358 | ||
359 | default: | |
360 | report(ERROR,"EvtGen") << "Spin2(_JA2)="<<_JA2<<" not supported!"<<endl; | |
361 | ::abort(); | |
362 | } | |
363 | ||
364 | ||
365 | switch(_JB2){ | |
366 | ||
367 | ||
368 | case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8: | |
369 | ||
370 | { | |
371 | ||
372 | int i,j,n; | |
373 | ||
374 | EvtSpinDensity R=p->getDaug(0)->rotateToHelicityBasis(phi,theta,-phi); | |
375 | ||
376 | n=R.getDim(); | |
377 | ||
378 | assert(n==_nB); | |
379 | ||
380 | for(i=0;i<n;i++){ | |
381 | for(j=0;j<n;j++){ | |
382 | _RB[i][j]=conj(R.get(i,j)); | |
383 | } | |
384 | } | |
385 | ||
386 | } | |
387 | ||
388 | break; | |
389 | ||
390 | default: | |
391 | report(ERROR,"EvtGen") << "Spin2(_JB2)="<<_JB2<<" not supported!"<<endl; | |
392 | ::abort(); | |
393 | } | |
394 | ||
395 | switch(_JC2){ | |
396 | ||
397 | case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8: | |
398 | ||
399 | { | |
400 | ||
401 | int i,j,n; | |
402 | ||
403 | EvtSpinDensity R=p->getDaug(1)->rotateToHelicityBasis(phi,EvtConst::pi+theta,phi-EvtConst::pi); | |
404 | ||
405 | n=R.getDim(); | |
406 | ||
407 | assert(n==_nC); | |
408 | ||
409 | for(i=0;i<n;i++){ | |
410 | for(j=0;j<n;j++){ | |
411 | _RC[i][j]=conj(R.get(i,j)); | |
412 | } | |
413 | } | |
414 | ||
415 | } | |
416 | ||
417 | break; | |
418 | ||
419 | default: | |
420 | report(ERROR,"EvtGen") << "Spin2(_JC2)="<<_JC2<<" not supported!"<<endl; | |
421 | ::abort(); | |
422 | } | |
423 | ||
424 | ||
425 | ||
426 | } | |
427 | ||
428 | ||
429 | void EvtEvalHelAmp::applyRotationMatrices(){ | |
430 | ||
431 | int ia,ib,ic,i; | |
432 | ||
433 | EvtComplex temp; | |
434 | ||
435 | ||
436 | ||
437 | for(ia=0;ia<_nA;ia++){ | |
438 | for(ib=0;ib<_nB;ib++){ | |
439 | for(ic=0;ic<_nC;ic++){ | |
440 | temp=0; | |
441 | for(i=0;i<_nC;i++){ | |
442 | temp+=_RC[i][ic]*_amp[ia][ib][i]; | |
443 | } | |
444 | _amp1[ia][ib][ic]=temp; | |
445 | } | |
446 | } | |
447 | } | |
448 | ||
449 | ||
450 | ||
451 | for(ia=0;ia<_nA;ia++){ | |
452 | for(ic=0;ic<_nC;ic++){ | |
453 | for(ib=0;ib<_nB;ib++){ | |
454 | temp=0; | |
455 | for(i=0;i<_nB;i++){ | |
456 | temp+=_RB[i][ib]*_amp1[ia][i][ic]; | |
457 | } | |
458 | _amp3[ia][ib][ic]=temp; | |
459 | } | |
460 | } | |
461 | } | |
462 | ||
463 | ||
464 | ||
465 | for(ib=0;ib<_nB;ib++){ | |
466 | for(ic=0;ic<_nC;ic++){ | |
467 | for(ia=0;ia<_nA;ia++){ | |
468 | temp=0; | |
469 | for(i=0;i<_nA;i++){ | |
470 | temp+=_RA[i][ia]*_amp3[i][ib][ic]; | |
471 | } | |
472 | _amp[ia][ib][ic]=temp; | |
473 | } | |
474 | } | |
475 | } | |
476 | ||
477 | ||
478 | } | |
479 | ||
480 | ||
481 | ||
482 | ||
483 | ||
484 | ||
485 | ||
486 | ||
487 | ||
488 | ||
489 | ||
490 |