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/EvtGenKine.hh"
28 : #include "EvtGenBase/EvtPDL.hh"
29 : #include "EvtGenBase/EvtVector4C.hh"
30 : #include "EvtGenBase/EvtTensor4C.hh"
31 : #include "EvtGenBase/EvtReport.hh"
32 : #include "EvtGenModels/EvtPartWave.hh"
33 : #include "EvtGenBase/EvtEvalHelAmp.hh"
34 : #include "EvtGenBase/EvtId.hh"
35 : #include <string>
36 : #include "EvtGenBase/EvtConst.hh"
37 : #include "EvtGenBase/EvtKine.hh"
38 : #include "EvtGenBase/EvtCGCoefSingle.hh"
39 : #include <algorithm>
40 : using std::endl;
41 0 : EvtPartWave::~EvtPartWave() {}
42 :
43 : std::string EvtPartWave::getName(){
44 :
45 0 : return "PARTWAVE";
46 :
47 : }
48 :
49 :
50 : EvtDecayBase* EvtPartWave::clone(){
51 :
52 0 : return new EvtPartWave;
53 :
54 0 : }
55 :
56 : void EvtPartWave::init(){
57 :
58 0 : checkNDaug(2);
59 :
60 : //find out how many states each particle have
61 0 : int _nA=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getParentId()));
62 0 : int _nB=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getDaug(0)));
63 0 : int _nC=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getDaug(1)));
64 :
65 0 : if (verbose()){
66 0 : report(Severity::Info,"EvtGen")<<"_nA,_nB,_nC:"
67 0 : <<_nA<<","<<_nB<<","<<_nC<<endl;
68 0 : }
69 :
70 : //find out what 2 times the spin is
71 0 : int _JA2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getParentId()));
72 0 : int _JB2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(0)));
73 0 : int _JC2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(1)));
74 :
75 0 : if (verbose()){
76 0 : report(Severity::Info,"EvtGen")<<"_JA2,_JB2,_JC2:"
77 0 : <<_JA2<<","<<_JB2<<","<<_JC2<<endl;
78 0 : }
79 :
80 :
81 : //allocate memory
82 0 : int* _lambdaA2=new int[_nA];
83 0 : int* _lambdaB2=new int[_nB];
84 0 : int* _lambdaC2=new int[_nC];
85 :
86 0 : EvtComplexPtr* _HBC=new EvtComplexPtr[_nB];
87 : int ib,ic;
88 0 : for(ib=0;ib<_nB;ib++){
89 0 : _HBC[ib]=new EvtComplex[_nC];
90 : }
91 :
92 :
93 : int i;
94 : //find the allowed helicities (actually 2*times the helicity!)
95 :
96 0 : fillHelicity(_lambdaA2,_nA,_JA2);
97 0 : fillHelicity(_lambdaB2,_nB,_JB2);
98 0 : fillHelicity(_lambdaC2,_nC,_JC2);
99 :
100 0 : if (verbose()){
101 0 : report(Severity::Info,"EvtGen")<<"Helicity states of particle A:"<<endl;
102 0 : for(i=0;i<_nA;i++){
103 0 : report(Severity::Info,"EvtGen")<<_lambdaA2[i]<<endl;
104 : }
105 :
106 0 : report(Severity::Info,"EvtGen")<<"Helicity states of particle B:"<<endl;
107 0 : for(i=0;i<_nB;i++){
108 0 : report(Severity::Info,"EvtGen")<<_lambdaB2[i]<<endl;
109 : }
110 :
111 0 : report(Severity::Info,"EvtGen")<<"Helicity states of particle C:"<<endl;
112 0 : for(i=0;i<_nC;i++){
113 0 : report(Severity::Info,"EvtGen")<<_lambdaC2[i]<<endl;
114 : }
115 :
116 0 : report(Severity::Info,"EvtGen")<<"Will now figure out the valid (M_LS) states:"<<endl;
117 :
118 0 : }
119 :
120 0 : int Lmin=std::max(_JA2-_JB2-_JC2,std::max(_JB2-_JA2-_JC2,_JC2-_JA2-_JB2));
121 0 : if (Lmin<0) Lmin=0;
122 : //int Lmin=_JA2-_JB2-_JC2;
123 0 : int Lmax=_JA2+_JB2+_JC2;
124 :
125 : int L;
126 :
127 : int _nPartialWaveAmp=0;
128 :
129 0 : int _nL[50];
130 0 : int _nS[50];
131 :
132 0 : for (L=Lmin;L<=Lmax;L+=2){
133 0 : int Smin=abs(L-_JA2);
134 0 : if (Smin<abs(_JB2-_JC2)) Smin=abs(_JB2-_JC2);
135 0 : int Smax=L+_JA2;
136 0 : if (Smax>abs(_JB2+_JC2)) Smax=abs(_JB2+_JC2);
137 : int S;
138 0 : for (S=Smin;S<=Smax;S+=2){
139 0 : _nL[_nPartialWaveAmp]=L;
140 0 : _nS[_nPartialWaveAmp]=S;
141 :
142 0 : _nPartialWaveAmp++;
143 0 : if (verbose()){
144 0 : report(Severity::Info,"EvtGen")<<"M["<<L<<"]["<<S<<"]"<<endl;
145 0 : }
146 : }
147 : }
148 :
149 0 : checkNArg(_nPartialWaveAmp*2);
150 :
151 : int argcounter=0;
152 :
153 0 : EvtComplex _M[50];
154 :
155 : double partampsqtot=0.0;
156 :
157 0 : for(i=0;i<_nPartialWaveAmp;i++){
158 0 : _M[i]=getArg(argcounter)*exp(EvtComplex(0.0,getArg(argcounter+1)));;
159 0 : argcounter+=2;
160 0 : partampsqtot+=abs2(_M[i]);
161 0 : if (verbose()){
162 0 : report(Severity::Info,"EvtGen")<<"M["<<_nL[i]<<"]["<<_nS[i]<<"]="<<_M[i]<<endl;
163 0 : }
164 : }
165 :
166 : //Now calculate the helicity amplitudes
167 :
168 : double helampsqtot=0.0;
169 :
170 0 : for(ib=0;ib<_nB;ib++){
171 0 : for(ic=0;ic<_nC;ic++){
172 0 : _HBC[ib][ic]=0.0;
173 0 : if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2){
174 0 : for(i=0;i<_nPartialWaveAmp;i++){
175 0 : int L=_nL[i];
176 0 : int S=_nS[i];
177 0 : int lambda2=_lambdaB2[ib];
178 0 : int lambda3=_lambdaC2[ic];
179 : int s1=_JA2;
180 : int s2=_JB2;
181 : int s3=_JC2;
182 0 : int m1=lambda2-lambda3;
183 0 : EvtCGCoefSingle c1(s2,s3);
184 0 : EvtCGCoefSingle c2(L,S);
185 :
186 0 : if (verbose()){
187 0 : report(Severity::Info,"EvtGen") << "s2,lambda2:"<<s2<<" "<<lambda2<<endl;
188 : }
189 : //fkw changes to satisfy KCC
190 0 : double fkwTmp = (L+1.0)/(s1+1.0);
191 :
192 0 : if (S>=abs(m1)){
193 :
194 0 : EvtComplex tmp=sqrt(fkwTmp)
195 0 : *c1.coef(S,m1,s2,s3,lambda2,-lambda3)
196 0 : *c2.coef(s1,m1,L,S,0,m1)*_M[i];
197 0 : _HBC[ib][ic]+=tmp;
198 0 : }
199 0 : }
200 0 : if (verbose()){
201 0 : report(Severity::Info,"EvtGen")<<"_HBC["<<ib<<"]["<<ic<<"]="<<_HBC[ib][ic]<<endl;
202 0 : }
203 : }
204 0 : helampsqtot+=abs2(_HBC[ib][ic]);
205 : }
206 : }
207 :
208 0 : if (fabs(helampsqtot-partampsqtot)/(helampsqtot+partampsqtot)>1e-6){
209 0 : report(Severity::Error,"EvtGen")<<"In EvtPartWave for decay "
210 0 : << EvtPDL::name(getParentId()) << " -> "
211 0 : << EvtPDL::name(getDaug(0)) << " "
212 0 : << EvtPDL::name(getDaug(1)) << std::endl;
213 0 : report(Severity::Error,"EvtGen")<<"With arguments: "<<std::endl;
214 0 : for(i=0;i*2<getNArg();i++){
215 0 : report(Severity::Error,"EvtGen") <<"M("<<_nL[i]<<","<<_nS[i]<<")="
216 : // <<getArg(2*i)<<" "<<getArg(2*i+1)<<std::endl;
217 0 : <<_M[i]<<std::endl;
218 : }
219 0 : report(Severity::Error,"EvtGen")<< "The total probability in the partwave basis is: "
220 0 : << partampsqtot << std::endl;
221 0 : report(Severity::Error,"EvtGen")<< "The total probability in the helamp basis is: "
222 0 : << helampsqtot << std::endl;
223 0 : report(Severity::Error,"EvtGen")<< "Most likely this is because the specified partwave amplitudes "
224 0 : << std::endl;
225 0 : report(Severity::Error,"EvtGen")<< "project onto unphysical helicities of photons or neutrinos. "
226 0 : << std::endl;
227 0 : report(Severity::Error,"EvtGen")<< "Seriously consider if your specified amplitudes are correct. "
228 0 : << std::endl;
229 :
230 :
231 0 : }
232 :
233 0 : _evalHelAmp=new EvtEvalHelAmp(getParentId(),
234 0 : getDaug(0),
235 0 : getDaug(1),
236 : _HBC);
237 :
238 0 : }
239 :
240 :
241 : void EvtPartWave::initProbMax(){
242 :
243 0 : double maxprob=_evalHelAmp->probMax();
244 :
245 0 : if (verbose()){
246 0 : report(Severity::Info,"EvtGen")<<"Calculated probmax"<<maxprob<<endl;
247 0 : }
248 :
249 0 : setProbMax(maxprob);
250 :
251 0 : }
252 :
253 :
254 : void EvtPartWave::decay( EvtParticle *p){
255 :
256 : //first generate simple phase space
257 0 : p->initializePhaseSpace(getNDaug(),getDaugs());
258 :
259 0 : _evalHelAmp->evalAmp(p,_amp2);
260 :
261 0 : return;
262 :
263 : }
264 :
265 :
266 :
267 : void EvtPartWave::fillHelicity(int* lambda2,int n,int J2){
268 :
269 : int i;
270 :
271 : //photon is special case!
272 0 : if (n==2&&J2==2) {
273 0 : lambda2[0]=2;
274 0 : lambda2[1]=-2;
275 0 : return;
276 : }
277 :
278 0 : assert(n==J2+1);
279 :
280 0 : for(i=0;i<n;i++){
281 0 : lambda2[i]=n-i*2-1;
282 : }
283 :
284 0 : return;
285 :
286 0 : }
287 :
288 :
289 :
290 :
291 :
292 :
293 :
294 :
295 :
296 :
297 :
|
