Line data Source code
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 0 : EvtEvalHelAmp::~EvtEvalHelAmp() {
39 :
40 : //deallocate memory
41 0 : delete [] _lambdaA2;
42 0 : delete [] _lambdaB2;
43 0 : delete [] _lambdaC2;
44 :
45 : int ia,ib,ic;
46 0 : for(ib=0;ib<_nB;ib++){
47 0 : delete [] _HBC[ib];
48 : }
49 :
50 0 : delete [] _HBC;
51 :
52 :
53 0 : for(ia=0;ia<_nA;ia++){
54 0 : delete [] _RA[ia];
55 : }
56 0 : delete [] _RA;
57 :
58 0 : for(ib=0;ib<_nB;ib++){
59 0 : delete [] _RB[ib];
60 : }
61 0 : delete [] _RB;
62 :
63 0 : for(ic=0;ic<_nC;ic++){
64 0 : delete [] _RC[ic];
65 : }
66 0 : delete [] _RC;
67 :
68 :
69 0 : for(ia=0;ia<_nA;ia++){
70 0 : for(ib=0;ib<_nB;ib++){
71 0 : delete [] _amp[ia][ib];
72 0 : delete [] _amp1[ia][ib];
73 0 : delete [] _amp3[ia][ib];
74 : }
75 0 : delete [] _amp[ia];
76 0 : delete [] _amp1[ia];
77 0 : delete [] _amp3[ia];
78 : }
79 :
80 0 : delete [] _amp;
81 0 : delete [] _amp1;
82 0 : delete [] _amp3;
83 :
84 0 : }
85 :
86 :
87 : EvtEvalHelAmp::EvtEvalHelAmp(EvtId idA,
88 : EvtId idB,
89 : EvtId idC,
90 0 : EvtComplexPtrPtr HBC){
91 :
92 :
93 0 : EvtSpinType::spintype typeA=EvtPDL::getSpinType(idA);
94 0 : EvtSpinType::spintype typeB=EvtPDL::getSpinType(idB);
95 0 : EvtSpinType::spintype typeC=EvtPDL::getSpinType(idC);
96 :
97 : //find out how many states each particle have
98 0 : _nA=EvtSpinType::getSpinStates(typeA);
99 0 : _nB=EvtSpinType::getSpinStates(typeB);
100 0 : _nC=EvtSpinType::getSpinStates(typeC);
101 :
102 : //find out what 2 times the spin is
103 0 : _JA2=EvtSpinType::getSpin2(typeA);
104 0 : _JB2=EvtSpinType::getSpin2(typeB);
105 0 : _JC2=EvtSpinType::getSpin2(typeC);
106 :
107 :
108 : //allocate memory
109 0 : _lambdaA2=new int[_nA];
110 0 : _lambdaB2=new int[_nB];
111 0 : _lambdaC2=new int[_nC];
112 :
113 0 : _HBC=new EvtComplexPtr[_nB];
114 : int ia,ib,ic;
115 0 : for(ib=0;ib<_nB;ib++){
116 0 : _HBC[ib]=new EvtComplex[_nC];
117 : }
118 :
119 :
120 0 : _RA=new EvtComplexPtr[_nA];
121 0 : for(ia=0;ia<_nA;ia++){
122 0 : _RA[ia]=new EvtComplex[_nA];
123 : }
124 0 : _RB=new EvtComplexPtr[_nB];
125 0 : for(ib=0;ib<_nB;ib++){
126 0 : _RB[ib]=new EvtComplex[_nB];
127 : }
128 0 : _RC=new EvtComplexPtr[_nC];
129 0 : for(ic=0;ic<_nC;ic++){
130 0 : _RC[ic]=new EvtComplex[_nC];
131 : }
132 :
133 0 : _amp=new EvtComplexPtrPtr[_nA];
134 0 : _amp1=new EvtComplexPtrPtr[_nA];
135 0 : _amp3=new EvtComplexPtrPtr[_nA];
136 0 : for(ia=0;ia<_nA;ia++){
137 0 : _amp[ia]=new EvtComplexPtr[_nB];
138 0 : _amp1[ia]=new EvtComplexPtr[_nB];
139 0 : _amp3[ia]=new EvtComplexPtr[_nB];
140 0 : for(ib=0;ib<_nB;ib++){
141 0 : _amp[ia][ib]=new EvtComplex[_nC];
142 0 : _amp1[ia][ib]=new EvtComplex[_nC];
143 0 : _amp3[ia][ib]=new EvtComplex[_nC];
144 : }
145 : }
146 :
147 : //find the allowed helicities (actually 2*times the helicity!)
148 :
149 0 : fillHelicity(_lambdaA2,_nA,_JA2,idA);
150 0 : fillHelicity(_lambdaB2,_nB,_JB2,idB);
151 0 : fillHelicity(_lambdaC2,_nC,_JC2,idC);
152 :
153 0 : for(ib=0;ib<_nB;ib++){
154 0 : for(ic=0;ic<_nC;ic++){
155 0 : _HBC[ib][ic]=HBC[ib][ic];
156 : }
157 : }
158 0 : }
159 :
160 :
161 :
162 :
163 :
164 :
165 : double EvtEvalHelAmp::probMax(){
166 :
167 0 : 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 0 : for(itheta=-10;itheta<=10;itheta++){
178 0 : theta=acos(0.099999*itheta);
179 0 : for(ia=0;ia<_nA;ia++){
180 : double prob=0.0;
181 0 : for(ib=0;ib<_nB;ib++){
182 0 : for(ic=0;ic<_nC;ic++){
183 0 : _amp[ia][ib][ic]=0.0;
184 0 : if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
185 0 : _amp[ia][ib][ic]=c*_HBC[ib][ic]*
186 0 : EvtdFunction::d(_JA2,_lambdaA2[ia],
187 0 : _lambdaB2[ib]-_lambdaC2[ic],theta);
188 0 : prob+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
189 0 : }
190 : }
191 : }
192 :
193 0 : prob*=sqrt(1.0*_nA);
194 :
195 0 : if (prob>maxprob) maxprob=prob;
196 :
197 : }
198 : }
199 :
200 0 : 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 0 : EvtVector4R pB=p->getDaug(0)->getP4();
210 :
211 0 : double theta=acos(pB.get(3)/pB.d3mag());
212 0 : double phi=atan2(pB.get(2),pB.get(1));
213 :
214 0 : double c=sqrt((_JA2+1)/(4*EvtConst::pi));
215 :
216 : int ia,ib,ic;
217 :
218 : double prob1=0.0;
219 :
220 0 : for(ia=0;ia<_nA;ia++){
221 0 : for(ib=0;ib<_nB;ib++){
222 0 : for(ic=0;ic<_nC;ic++){
223 0 : _amp[ia][ib][ic]=0.0;
224 0 : if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
225 0 : double dfun=EvtdFunction::d(_JA2,_lambdaA2[ia],
226 : _lambdaB2[ib]-_lambdaC2[ic],theta);
227 :
228 0 : _amp[ia][ib][ic]=c*_HBC[ib][ic]*
229 0 : exp(EvtComplex(0.0,phi*0.5*(_lambdaA2[ia]-_lambdaB2[ib]+
230 0 : _lambdaC2[ic])))*dfun;
231 0 : }
232 0 : prob1+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
233 : }
234 : }
235 : }
236 :
237 0 : setUpRotationMatrices(p,theta,phi);
238 :
239 0 : applyRotationMatrices();
240 :
241 : double prob2=0.0;
242 :
243 0 : for(ia=0;ia<_nA;ia++){
244 0 : for(ib=0;ib<_nB;ib++){
245 0 : for(ic=0;ic<_nC;ic++){
246 0 : prob2+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
247 0 : if (_nA==1){
248 0 : if (_nB==1){
249 0 : if (_nC==1){
250 0 : amp.vertex(_amp[ia][ib][ic]);
251 0 : }
252 : else{
253 0 : amp.vertex(ic,_amp[ia][ib][ic]);
254 : }
255 : }
256 : else{
257 0 : if (_nC==1){
258 0 : amp.vertex(ib,_amp[ia][ib][ic]);
259 0 : }
260 : else{
261 0 : amp.vertex(ib,ic,_amp[ia][ib][ic]);
262 : }
263 : }
264 : }else{
265 0 : if (_nB==1){
266 0 : if (_nC==1){
267 0 : amp.vertex(ia,_amp[ia][ib][ic]);
268 0 : }
269 : else{
270 0 : amp.vertex(ia,ic,_amp[ia][ib][ic]);
271 : }
272 : }
273 : else{
274 0 : if (_nC==1){
275 0 : amp.vertex(ia,ib,_amp[ia][ib][ic]);
276 0 : }
277 : else{
278 0 : amp.vertex(ia,ib,ic,_amp[ia][ib][ic]);
279 : }
280 : }
281 : }
282 : }
283 : }
284 : }
285 :
286 0 : if (fabs(prob1-prob2)>0.000001*prob1){
287 0 : report(Severity::Info,"EvtGen") << "prob1,prob2:"<<prob1<<" "<<prob2<<endl;
288 0 : ::abort();
289 : }
290 :
291 : return ;
292 :
293 0 : }
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 0 : if (n==2&&J2==2) {
302 0 : lambda2[0]=2;
303 0 : lambda2[1]=-2;
304 0 : return;
305 : }
306 :
307 : //and so is the neutrino!
308 0 : if (n==1&&J2==1) {
309 0 : if (EvtPDL::getStdHep(id)>0){
310 : //particle i.e. lefthanded
311 0 : lambda2[0]=-1;
312 0 : }else{
313 : //anti particle i.e. righthanded
314 0 : lambda2[0]=1;
315 : }
316 0 : return;
317 : }
318 :
319 :
320 0 : assert(n==J2+1);
321 :
322 0 : for(i=0;i<n;i++){
323 0 : lambda2[i]=n-i*2-1;
324 : }
325 :
326 0 : return;
327 :
328 0 : }
329 :
330 :
331 : void EvtEvalHelAmp::setUpRotationMatrices(EvtParticle* p,double theta, double phi){
332 :
333 0 : 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 0 : EvtSpinDensity R=p->rotateToHelicityBasis();
340 :
341 :
342 : int i,j,n;
343 :
344 0 : n=R.getDim();
345 :
346 0 : assert(n==_nA);
347 :
348 :
349 0 : for(i=0;i<n;i++){
350 0 : for(j=0;j<n;j++){
351 0 : _RA[i][j]=R.get(i,j);
352 : }
353 : }
354 :
355 0 : }
356 :
357 : break;
358 :
359 : default:
360 0 : report(Severity::Error,"EvtGen") << "Spin2(_JA2)="<<_JA2<<" not supported!"<<endl;
361 0 : ::abort();
362 : }
363 :
364 :
365 0 : 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 0 : EvtSpinDensity R=p->getDaug(0)->rotateToHelicityBasis(phi,theta,-phi);
375 :
376 0 : n=R.getDim();
377 :
378 0 : assert(n==_nB);
379 :
380 0 : for(i=0;i<n;i++){
381 0 : for(j=0;j<n;j++){
382 0 : _RB[i][j]=conj(R.get(i,j));
383 : }
384 : }
385 :
386 0 : }
387 :
388 : break;
389 :
390 : default:
391 0 : report(Severity::Error,"EvtGen") << "Spin2(_JB2)="<<_JB2<<" not supported!"<<endl;
392 0 : ::abort();
393 : }
394 :
395 0 : 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 0 : EvtSpinDensity R=p->getDaug(1)->rotateToHelicityBasis(phi,EvtConst::pi+theta,phi-EvtConst::pi);
404 :
405 0 : n=R.getDim();
406 :
407 0 : assert(n==_nC);
408 :
409 0 : for(i=0;i<n;i++){
410 0 : for(j=0;j<n;j++){
411 0 : _RC[i][j]=conj(R.get(i,j));
412 : }
413 : }
414 :
415 0 : }
416 :
417 : break;
418 :
419 : default:
420 0 : report(Severity::Error,"EvtGen") << "Spin2(_JC2)="<<_JC2<<" not supported!"<<endl;
421 0 : ::abort();
422 : }
423 :
424 :
425 :
426 0 : }
427 :
428 :
429 : void EvtEvalHelAmp::applyRotationMatrices(){
430 :
431 : int ia,ib,ic,i;
432 :
433 0 : EvtComplex temp;
434 :
435 :
436 :
437 0 : for(ia=0;ia<_nA;ia++){
438 0 : for(ib=0;ib<_nB;ib++){
439 0 : for(ic=0;ic<_nC;ic++){
440 0 : temp=0;
441 0 : for(i=0;i<_nC;i++){
442 0 : temp+=_RC[i][ic]*_amp[ia][ib][i];
443 : }
444 0 : _amp1[ia][ib][ic]=temp;
445 : }
446 : }
447 : }
448 :
449 :
450 :
451 0 : for(ia=0;ia<_nA;ia++){
452 0 : for(ic=0;ic<_nC;ic++){
453 0 : for(ib=0;ib<_nB;ib++){
454 0 : temp=0;
455 0 : for(i=0;i<_nB;i++){
456 0 : temp+=_RB[i][ib]*_amp1[ia][i][ic];
457 : }
458 0 : _amp3[ia][ib][ic]=temp;
459 : }
460 : }
461 : }
462 :
463 :
464 :
465 0 : for(ib=0;ib<_nB;ib++){
466 0 : for(ic=0;ic<_nC;ic++){
467 0 : for(ia=0;ia<_nA;ia++){
468 0 : temp=0;
469 0 : for(i=0;i<_nA;i++){
470 0 : temp+=_RA[i][ia]*_amp3[i][ib][ic];
471 : }
472 0 : _amp[ia][ib][ic]=temp;
473 : }
474 : }
475 : }
476 :
477 :
478 0 : }
479 :
480 :
481 :
482 :
483 :
484 :
485 :
486 :
487 :
488 :
489 :
490 :
|
