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) 1998 Caltech, UCSB
10 : //
11 : // Module: EvtHelAmp.cc
12 : //
13 : // Description: Decay model for implementation of generic 2 body
14 : // decay specified by the helicity amplitudes
15 : //
16 : //
17 : // Modification history:
18 : //
19 : // RYD March 14, 1999 Module created
20 : //
21 : //------------------------------------------------------------------------
22 : //
23 : #include "EvtGenBase/EvtPatches.hh"
24 : #include <stdlib.h>
25 : #include "EvtGenBase/EvtParticle.hh"
26 : #include "EvtGenBase/EvtGenKine.hh"
27 : #include "EvtGenBase/EvtPDL.hh"
28 : #include "EvtGenBase/EvtReport.hh"
29 : #include "EvtGenModels/EvtHelAmp.hh"
30 : #include "EvtGenBase/EvtId.hh"
31 : #include <string>
32 : #include "EvtGenBase/EvtConst.hh"
33 : #include "EvtGenBase/EvtEvalHelAmp.hh"
34 : using std::endl;
35 :
36 :
37 0 : EvtHelAmp::~EvtHelAmp() {
38 :
39 0 : delete _evalHelAmp;
40 :
41 0 : }
42 :
43 : std::string EvtHelAmp::getName(){
44 :
45 0 : return "HELAMP";
46 :
47 : }
48 :
49 :
50 : EvtDecayBase* EvtHelAmp::clone(){
51 :
52 0 : return new EvtHelAmp;
53 :
54 0 : }
55 :
56 : void EvtHelAmp::init(){
57 :
58 0 : checkNDaug(2);
59 :
60 :
61 : //find out how many states each particle have
62 0 : int _nA=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getParentId()));
63 0 : int _nB=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getDaug(0)));
64 0 : int _nC=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getDaug(1)));
65 :
66 0 : if (verbose()){
67 0 : report(Severity::Info,"EvtGen")<<"_nA,_nB,_nC:"
68 0 : <<_nA<<","<<_nB<<","<<_nC<<endl;
69 0 : }
70 :
71 : //find out what 2 times the spin is
72 0 : int _JA2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getParentId()));
73 0 : int _JB2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(0)));
74 0 : int _JC2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(1)));
75 :
76 0 : if (verbose()){
77 0 : report(Severity::Info,"EvtGen")<<"_JA2,_JB2,_JC2:"
78 0 : <<_JA2<<","<<_JB2<<","<<_JC2<<endl;
79 0 : }
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 0 : for(int ib=0;ib<_nB;ib++){
88 0 : _HBC[ib]=new EvtComplex[_nC];
89 : }
90 :
91 : int i;
92 : //find the allowed helicities (actually 2*times the helicity!)
93 :
94 0 : fillHelicity(_lambdaA2,_nA,_JA2,getParentId());
95 0 : fillHelicity(_lambdaB2,_nB,_JB2,getDaug(0));
96 0 : fillHelicity(_lambdaC2,_nC,_JC2,getDaug(1));
97 :
98 0 : if (verbose()){
99 0 : report(Severity::Info,"EvtGen")<<"Helicity states of particle A:"<<endl;
100 0 : for(i=0;i<_nA;i++){
101 0 : report(Severity::Info,"EvtGen")<<_lambdaA2[i]<<endl;
102 : }
103 :
104 0 : report(Severity::Info,"EvtGen")<<"Helicity states of particle B:"<<endl;
105 0 : for(i=0;i<_nB;i++){
106 0 : report(Severity::Info,"EvtGen")<<_lambdaB2[i]<<endl;
107 : }
108 :
109 0 : report(Severity::Info,"EvtGen")<<"Helicity states of particle C:"<<endl;
110 0 : for(i=0;i<_nC;i++){
111 0 : report(Severity::Info,"EvtGen")<<_lambdaC2[i]<<endl;
112 : }
113 : }
114 :
115 : //now read in the helicity amplitudes
116 :
117 : int argcounter=0;
118 :
119 0 : for(int ib=0;ib<_nB;ib++){
120 0 : for(int ic=0;ic<_nC;ic++){
121 0 : _HBC[ib][ic]=0.0;
122 0 : if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) argcounter+=2;
123 : }
124 : }
125 :
126 0 : checkNArg(argcounter);
127 :
128 : argcounter=0;
129 :
130 0 : for(int ib=0;ib<_nB;ib++){
131 0 : for(int ic=0;ic<_nC;ic++){
132 0 : if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
133 0 : _HBC[ib][ic]=getArg(argcounter)*exp(EvtComplex(0.0,getArg(argcounter+1)));;
134 0 : argcounter+=2;
135 0 : if (verbose()){
136 0 : report(Severity::Info,"EvtGen")<<"_HBC["<<ib<<"]["<<ic<<"]="
137 0 : <<_HBC[ib][ic]<<endl;
138 0 : }
139 : }
140 : }
141 : }
142 :
143 0 : _evalHelAmp=new EvtEvalHelAmp(getParentId(),
144 0 : getDaug(0),
145 0 : getDaug(1),
146 : _HBC);
147 :
148 : // Note: these are not class data members but local variables.
149 0 : delete [] _lambdaA2;
150 0 : delete [] _lambdaB2;
151 0 : delete [] _lambdaC2;
152 0 : for(int ib=0;ib<_nB;ib++){
153 0 : delete [] _HBC[ib];
154 : }
155 0 : delete [] _HBC; // _HBC is copied in ctor of EvtEvalHelAmp above.
156 :
157 0 : }
158 :
159 :
160 : void EvtHelAmp::initProbMax(){
161 :
162 0 : double maxprob=_evalHelAmp->probMax();
163 :
164 0 : if (verbose()){
165 0 : report(Severity::Info,"EvtGen")<<"Calculated probmax"<<maxprob<<endl;
166 0 : }
167 :
168 0 : setProbMax(maxprob);
169 :
170 0 : }
171 :
172 :
173 : void EvtHelAmp::decay( EvtParticle *p){
174 :
175 : //first generate simple phase space
176 0 : p->initializePhaseSpace(getNDaug(),getDaugs());
177 :
178 0 : _evalHelAmp->evalAmp(p,_amp2);
179 :
180 0 : return ;
181 :
182 : }
183 :
184 :
185 : void EvtHelAmp::fillHelicity(int* lambda2,int n,int J2, EvtId id){
186 :
187 : int i;
188 :
189 : //photon is special case!
190 0 : if (n==2&&J2==2) {
191 0 : lambda2[0]=2;
192 0 : lambda2[1]=-2;
193 0 : return;
194 : }
195 :
196 : //and so is the neutrino!
197 0 : if (n==1&&J2==1) {
198 0 : if (EvtPDL::getStdHep(id)>0){
199 : //particle i.e. lefthanded
200 0 : lambda2[0]=-1;
201 0 : }else{
202 : //anti particle i.e. righthanded
203 0 : lambda2[0]=1;
204 : }
205 0 : return;
206 : }
207 :
208 0 : assert(n==J2+1);
209 :
210 0 : for(i=0;i<n;i++){
211 0 : lambda2[i]=n-i*2-1;
212 : }
213 :
214 0 : return;
215 :
216 0 : }
217 :
218 :
219 :
220 :
221 :
222 :
223 :
|
