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: 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 0 : EvtSpinDensity rho;
43 : double prob,prob_max;
44 :
45 0 : _amp2.init(p->getId(),getNDaug(),getDaugs());
46 :
47 : do{
48 :
49 0 : _daugsDecayedByParentModel=false;
50 0 : _weight = 1.0;
51 0 : decay(p);
52 :
53 0 : rho=_amp2.getSpinDensity();
54 :
55 0 : prob=p->getSpinDensityForward().normalizedProb(rho);
56 :
57 0 : if (prob<0.0) {
58 0 : report(Severity::Error,"EvtGen")<<"Negative prob:"<<p->getId().getId()
59 0 : <<" "<<p->getChannel()<<endl;
60 :
61 0 : report(Severity::Error,"EvtGen") << "rho_forward:"<<endl;
62 0 : report(Severity::Error,"EvtGen") << p->getSpinDensityForward();
63 0 : report(Severity::Error,"EvtGen") << "rho decay:"<<endl;
64 0 : report(Severity::Error,"EvtGen") << rho <<endl;
65 : }
66 :
67 0 : if (prob!=prob) {
68 :
69 0 : report(Severity::Debug,"EvtGen") << "Forward density matrix:"<<endl;
70 0 : report(Severity::Debug,"EvtGen") << p->getSpinDensityForward();
71 :
72 0 : report(Severity::Debug,"EvtGen") << "Decay density matrix:"<<endl;
73 0 : report(Severity::Debug,"EvtGen") << rho;
74 :
75 0 : report(Severity::Debug,"EvtGen") << "prob:"<<prob<<endl;
76 :
77 0 : report(Severity::Debug,"EvtGen") << "Particle:"
78 0 : <<EvtPDL::name(p->getId()).c_str()<<endl;
79 0 : report(Severity::Debug,"EvtGen") << "channel :"<<p->getChannel()<<endl;
80 0 : report(Severity::Debug,"EvtGen") << "Momentum:" << p->getP4() << " " << p->mass() << endl;
81 0 : if( p->getParent()!=0){
82 0 : report(Severity::Debug,"EvtGen") << "parent:"
83 0 : <<EvtPDL::name(
84 0 : p->getParent()->getId()).c_str()<<endl;
85 0 : report(Severity::Debug,"EvtGen") << "parent channel :"
86 0 : <<p->getParent()->getChannel()<<endl;
87 :
88 : size_t i;
89 0 : report(Severity::Debug,"EvtGen") << "parent daughters :";
90 0 : for (i=0;i<p->getParent()->getNDaug();i++){
91 0 : report(Severity::Debug,"") << EvtPDL::name(
92 0 : p->getParent()->getDaug(i)->getId()).c_str()
93 0 : << " ";
94 : }
95 0 : report(Severity::Debug,"") << endl;
96 :
97 0 : report(Severity::Debug,"EvtGen") << "daughters :";
98 0 : for (size_t i=0;i<p->getNDaug();i++){
99 0 : report(Severity::Debug,"") << EvtPDL::name(
100 0 : p->getDaug(i)->getId()).c_str()
101 0 : << " ";
102 : }
103 0 : report(Severity::Debug,"") << endl;
104 :
105 0 : report(Severity::Debug,"EvtGen") << "daughter momenta :" << endl;;
106 0 : for (size_t i=0;i<p->getNDaug();i++){
107 0 : report(Severity::Debug,"") << p->getDaug(i)->getP4() << " " << p->getDaug(i)->mass();
108 0 : report(Severity::Debug,"") << endl;
109 : }
110 :
111 0 : }
112 : }
113 :
114 :
115 0 : prob/=_weight;
116 :
117 0 : prob_max = getProbMax(prob);
118 0 : p->setDecayProb(prob/prob_max);
119 :
120 0 : more=prob<EvtRandom::Flat(prob_max);
121 :
122 0 : ntimes--;
123 :
124 0 : }while(ntimes&&more);
125 :
126 0 : if (ntimes==0){
127 0 : report(Severity::Debug,"EvtGen") << "Tried accept/reject: 10000"
128 0 : <<" times, and rejected all the times!"<<endl;
129 :
130 0 : report(Severity::Debug,"EvtGen")<<p->getSpinDensityForward()<<endl;
131 0 : report(Severity::Debug,"EvtGen") << "Is therefore accepting the last event!"<<endl;
132 0 : report(Severity::Debug,"EvtGen") << "Decay of particle:"<<
133 0 : EvtPDL::name(p->getId()).c_str()<<"(channel:"<<
134 0 : p->getChannel()<<") with mass "<<p->mass()<<endl;
135 :
136 0 : for(size_t ii=0;ii<p->getNDaug();ii++){
137 0 : report(Severity::Debug,"EvtGen") <<"Daughter "<<ii<<":"<<
138 0 : EvtPDL::name(p->getDaug(ii)->getId()).c_str()<<" with mass "<<
139 0 : p->getDaug(ii)->mass()<<endl;
140 : }
141 0 : }
142 :
143 0 : EvtSpinDensity rho_list[10];
144 :
145 0 : rho_list[0]=p->getSpinDensityForward();
146 :
147 0 : EvtAmp ampcont;
148 :
149 0 : if (_amp2._pstates!=1){
150 0 : ampcont=_amp2.contract(0,p->getSpinDensityForward());
151 0 : }
152 : else{
153 0 : 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 0 : if ( !daugsDecayedByParentModel() ){
162 :
163 0 : if(recursive) {
164 :
165 0 : for(size_t i=0;i<p->getNDaug();i++){
166 :
167 0 : rho.setDim(_amp2.dstates[i]);
168 :
169 0 : if (_amp2.dstates[i]==1) {
170 0 : rho.set(0,0,EvtComplex(1.0,0.0));
171 : }
172 : else{
173 0 : rho=ampcont.contract(_amp2._dnontrivial[i],_amp2);
174 : }
175 :
176 0 : if (!rho.check()) {
177 :
178 0 : report(Severity::Error,"EvtGen") << "-------start error-------"<<endl;
179 0 : report(Severity::Error,"EvtGen")<<"forward rho failed Check:"<<
180 0 : EvtPDL::name(p->getId()).c_str()<<" "<<p->getChannel()<<" "<<i<<endl;
181 :
182 0 : p->printTree();
183 :
184 0 : for (size_t idaug = 0; idaug < p->getNDaug(); idaug++) {
185 0 : EvtParticle* daughter = p->getDaug(idaug);
186 0 : if (daughter != 0) {daughter->printTree();}
187 : }
188 :
189 0 : EvtParticle* pParent = p->getParent();
190 0 : if (pParent != 0) {
191 0 : report(Severity::Error,"EvtGen")<<"Parent:"<<EvtPDL::name(pParent->getId()).c_str()<<endl;
192 :
193 0 : EvtParticle* grandParent = pParent->getParent();
194 :
195 0 : if (grandParent != 0) {
196 0 : report(Severity::Error,"EvtGen")<<"GrandParent:"<<EvtPDL::name(grandParent->getId()).c_str()<<endl;
197 0 : }
198 0 : }
199 :
200 0 : report(Severity::Error,"EvtGen") << " EvtSpinDensity rho: " << rho;
201 :
202 0 : _amp2.dump();
203 :
204 0 : for(size_t ii=0;ii<i+1;ii++){
205 0 : report(Severity::Error,"EvtGen") << "rho_list[" << ii << "] = " << rho_list[ii];
206 : }
207 :
208 0 : report(Severity::Error,"EvtGen") << "-------Done with error-------"<<endl;
209 :
210 0 : }
211 :
212 0 : p->getDaug(i)->setSpinDensityForward(rho);
213 0 : p->getDaug(i)->decay();
214 :
215 0 : rho_list[i+1]=p->getDaug(i)->getSpinDensityBackward();
216 :
217 0 : if (_amp2.dstates[i]!=1){
218 0 : ampcont=ampcont.contract(_amp2._dnontrivial[i],rho_list[i+1]);
219 : }
220 :
221 :
222 : }
223 :
224 0 : p->setSpinDensityBackward(_amp2.getBackwardSpinDensity(rho_list));
225 :
226 :
227 0 : if (!p->getSpinDensityBackward().check()) {
228 :
229 0 : report(Severity::Error,"EvtGen")<<"rho_backward failed Check"<<
230 0 : p->getId().getId()<<" "<<p->getChannel()<<endl;
231 :
232 0 : report(Severity::Error,"EvtGen") << p->getSpinDensityBackward();
233 :
234 0 : }
235 : }
236 : }
237 :
238 :
239 0 : if (getPHOTOS() || EvtRadCorr::alwaysRadCorr()) {
240 0 : int n_daug_orig=p->getNDaug();
241 0 : EvtRadCorr::doRadCorr(p);
242 0 : int n_daug_new=p->getNDaug();
243 0 : for (int i=n_daug_orig;i<n_daug_new;i++){
244 0 : p->getDaug(i)->decay();
245 : }
246 0 : }
247 :
248 0 : }
249 :
250 :
251 :
252 :
253 :
254 :
255 :
256 :
257 :
258 :
|
