]>
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) 1998 Caltech, UCSB | |
10 | // | |
11 | // Module: EvtGen/EvtDecayAmp.cc | |
12 | // | |
13 | // Description: Baseclass for models that calculates amplitudes | |
14 | // | |
15 | // Modification history: | |
16 | // | |
17 | // DJL/RYD August 11, 1998 Module created | |
18 | // | |
19 | //------------------------------------------------------------------------ | |
20 | #include "EvtGenBase/EvtPatches.hh" | |
21 | ||
22 | ||
23 | ||
24 | #include "EvtGenBase/EvtDecayBase.hh" | |
25 | #include "EvtGenBase/EvtDecayAmp.hh" | |
26 | #include "EvtGenBase/EvtParticle.hh" | |
27 | #include "EvtGenBase/EvtPDL.hh" | |
28 | #include "EvtGenBase/EvtRandom.hh" | |
29 | #include "EvtGenBase/EvtRadCorr.hh" | |
30 | #include "EvtGenBase/EvtAmp.hh" | |
31 | #include "EvtGenBase/EvtReport.hh" | |
32 | using std::endl; | |
33 | ||
34 | ||
35 | void EvtDecayAmp::makeDecay(EvtParticle* p, bool recursive){ | |
36 | ||
37 | //original default value | |
38 | int ntimes=10000; | |
39 | ||
40 | int more; | |
41 | ||
42 | EvtSpinDensity rho; | |
43 | double prob,prob_max; | |
44 | ||
45 | _amp2.init(p->getId(),getNDaug(),getDaugs()); | |
46 | ||
47 | do{ | |
48 | ||
49 | _daugsDecayedByParentModel=false; | |
50 | _weight = 1.0; | |
51 | decay(p); | |
52 | ||
53 | rho=_amp2.getSpinDensity(); | |
54 | ||
55 | prob=p->getSpinDensityForward().normalizedProb(rho); | |
56 | ||
57 | if (prob<0.0) { | |
58 | report(ERROR,"EvtGen")<<"Negative prob:"<<p->getId().getId() | |
59 | <<" "<<p->getChannel()<<endl; | |
60 | ||
61 | report(ERROR,"EvtGen") << "rho_forward:"<<endl; | |
62 | report(ERROR,"EvtGen") << p->getSpinDensityForward(); | |
63 | report(ERROR,"EvtGen") << "rho decay:"<<endl; | |
64 | report(ERROR,"EvtGen") << rho <<endl; | |
65 | } | |
66 | ||
67 | if (prob!=prob) { | |
68 | ||
69 | report(DEBUG,"EvtGen") << "Forward density matrix:"<<endl; | |
0ca57c2f | 70 | report(DEBUG,"EvtGen") << p->getSpinDensityForward(); |
da0e9ce3 | 71 | |
72 | report(DEBUG,"EvtGen") << "Decay density matrix:"<<endl; | |
0ca57c2f | 73 | report(DEBUG,"EvtGen") << rho; |
da0e9ce3 | 74 | |
75 | report(DEBUG,"EvtGen") << "prob:"<<prob<<endl; | |
76 | ||
77 | report(DEBUG,"EvtGen") << "Particle:" | |
78 | <<EvtPDL::name(p->getId()).c_str()<<endl; | |
79 | report(DEBUG,"EvtGen") << "channel :"<<p->getChannel()<<endl; | |
80 | report(DEBUG,"EvtGen") << "Momentum:" << p->getP4() << " " << p->mass() << endl; | |
81 | if( p->getParent()!=0){ | |
82 | report(DEBUG,"EvtGen") << "parent:" | |
83 | <<EvtPDL::name( | |
84 | p->getParent()->getId()).c_str()<<endl; | |
85 | report(DEBUG,"EvtGen") << "parent channel :" | |
86 | <<p->getParent()->getChannel()<<endl; | |
87 | ||
88 | size_t i; | |
89 | report(DEBUG,"EvtGen") << "parent daughters :"; | |
90 | for (i=0;i<p->getParent()->getNDaug();i++){ | |
91 | report(DEBUG,"") << EvtPDL::name( | |
92 | p->getParent()->getDaug(i)->getId()).c_str() | |
93 | << " "; | |
94 | } | |
95 | report(DEBUG,"") << endl; | |
96 | ||
97 | report(DEBUG,"EvtGen") << "daughters :"; | |
98 | for (size_t i=0;i<p->getNDaug();i++){ | |
99 | report(DEBUG,"") << EvtPDL::name( | |
100 | p->getDaug(i)->getId()).c_str() | |
101 | << " "; | |
102 | } | |
103 | report(DEBUG,"") << endl; | |
104 | ||
105 | report(DEBUG,"EvtGen") << "daughter momenta :" << endl;; | |
106 | for (size_t i=0;i<p->getNDaug();i++){ | |
107 | report(DEBUG,"") << p->getDaug(i)->getP4() << " " << p->getDaug(i)->mass(); | |
108 | report(DEBUG,"") << endl; | |
109 | } | |
110 | ||
111 | } | |
112 | } | |
113 | ||
114 | ||
115 | prob/=_weight; | |
116 | ||
117 | prob_max = getProbMax(prob); | |
118 | p->setDecayProb(prob/prob_max); | |
119 | ||
120 | more=prob<EvtRandom::Flat(prob_max); | |
121 | ||
122 | ntimes--; | |
123 | ||
124 | }while(ntimes&&more); | |
125 | ||
126 | if (ntimes==0){ | |
127 | report(DEBUG,"EvtGen") << "Tried accept/reject: 10000" | |
128 | <<" times, and rejected all the times!"<<endl; | |
129 | ||
130 | report(DEBUG,"EvtGen")<<p->getSpinDensityForward()<<endl; | |
131 | report(DEBUG,"EvtGen") << "Is therefore accepting the last event!"<<endl; | |
132 | report(DEBUG,"EvtGen") << "Decay of particle:"<< | |
133 | EvtPDL::name(p->getId()).c_str()<<"(channel:"<< | |
134 | p->getChannel()<<") with mass "<<p->mass()<<endl; | |
135 | ||
136 | for(size_t ii=0;ii<p->getNDaug();ii++){ | |
137 | report(DEBUG,"EvtGen") <<"Daughter "<<ii<<":"<< | |
138 | EvtPDL::name(p->getDaug(ii)->getId()).c_str()<<" with mass "<< | |
139 | p->getDaug(ii)->mass()<<endl; | |
140 | } | |
141 | } | |
142 | ||
143 | EvtSpinDensity rho_list[10]; | |
144 | ||
145 | rho_list[0]=p->getSpinDensityForward(); | |
146 | ||
147 | EvtAmp ampcont; | |
148 | ||
149 | if (_amp2._pstates!=1){ | |
150 | ampcont=_amp2.contract(0,p->getSpinDensityForward()); | |
151 | } | |
152 | else{ | |
153 | ampcont=_amp2; | |
154 | } | |
155 | ||
156 | ||
157 | // it may be that the parent decay model has already | |
158 | // done the decay - this should be rare and the | |
159 | // model better know what it is doing.. | |
160 | ||
161 | if ( !daugsDecayedByParentModel() ){ | |
162 | ||
163 | if(recursive) { | |
164 | ||
165 | for(size_t i=0;i<p->getNDaug();i++){ | |
166 | ||
167 | rho.setDim(_amp2.dstates[i]); | |
168 | ||
169 | if (_amp2.dstates[i]==1) { | |
170 | rho.set(0,0,EvtComplex(1.0,0.0)); | |
171 | } | |
172 | else{ | |
173 | rho=ampcont.contract(_amp2._dnontrivial[i],_amp2); | |
174 | } | |
175 | ||
176 | if (!rho.check()) { | |
177 | ||
178 | report(ERROR,"EvtGen") << "-------start error-------"<<endl; | |
179 | report(ERROR,"EvtGen")<<"forward rho failed Check:"<< | |
180 | EvtPDL::name(p->getId()).c_str()<<" "<<p->getChannel()<<" "<<i<<endl; | |
181 | ||
0ca57c2f | 182 | p->printTree(); |
183 | ||
184 | for (size_t idaug = 0; idaug < p->getNDaug(); idaug++) { | |
185 | EvtParticle* daughter = p->getDaug(idaug); | |
186 | if (daughter != 0) {daughter->printTree();} | |
187 | } | |
188 | ||
189 | EvtParticle* pParent = p->getParent(); | |
190 | if (pParent != 0) { | |
191 | report(ERROR,"EvtGen")<<"Parent:"<<EvtPDL::name(pParent->getId()).c_str()<<endl; | |
192 | ||
193 | EvtParticle* grandParent = pParent->getParent(); | |
194 | ||
195 | if (grandParent != 0) { | |
196 | report(ERROR,"EvtGen")<<"GrandParent:"<<EvtPDL::name(grandParent->getId()).c_str()<<endl; | |
197 | } | |
198 | } | |
199 | ||
200 | report(ERROR,"EvtGen") << " EvtSpinDensity rho: " << rho; | |
da0e9ce3 | 201 | |
202 | _amp2.dump(); | |
0ca57c2f | 203 | |
da0e9ce3 | 204 | for(size_t ii=0;ii<i+1;ii++){ |
0ca57c2f | 205 | report(ERROR,"EvtGen") << "rho_list[" << ii << "] = " << rho_list[ii]; |
da0e9ce3 | 206 | } |
0ca57c2f | 207 | |
da0e9ce3 | 208 | report(ERROR,"EvtGen") << "-------Done with error-------"<<endl; |
0ca57c2f | 209 | |
da0e9ce3 | 210 | } |
211 | ||
212 | p->getDaug(i)->setSpinDensityForward(rho); | |
213 | p->getDaug(i)->decay(); | |
214 | ||
215 | rho_list[i+1]=p->getDaug(i)->getSpinDensityBackward(); | |
216 | ||
217 | if (_amp2.dstates[i]!=1){ | |
218 | ampcont=ampcont.contract(_amp2._dnontrivial[i],rho_list[i+1]); | |
219 | } | |
220 | ||
221 | ||
222 | } | |
223 | ||
224 | p->setSpinDensityBackward(_amp2.getBackwardSpinDensity(rho_list)); | |
225 | ||
226 | ||
227 | if (!p->getSpinDensityBackward().check()) { | |
228 | ||
229 | report(ERROR,"EvtGen")<<"rho_backward failed Check"<< | |
230 | p->getId().getId()<<" "<<p->getChannel()<<endl; | |
231 | ||
232 | report(ERROR,"EvtGen") << p->getSpinDensityBackward(); | |
233 | ||
234 | } | |
235 | } | |
236 | } | |
237 | ||
238 | ||
239 | if (getPHOTOS() || EvtRadCorr::alwaysRadCorr()) { | |
240 | int n_daug_orig=p->getNDaug(); | |
241 | EvtRadCorr::doRadCorr(p); | |
242 | int n_daug_new=p->getNDaug(); | |
243 | for (int i=n_daug_orig;i<n_daug_new;i++){ | |
244 | p->getDaug(i)->decay(); | |
245 | } | |
246 | } | |
247 | ||
248 | } | |
249 | ||
250 | ||
251 | ||
252 | ||
253 | ||
254 | ||
255 | ||
256 | ||
257 | ||
258 |