Fix for definitions for CINT
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenBase / EvtDecayAmp.cpp
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;
70       report(DEBUG,"EvtGen") << p->getSpinDensityForward();
71
72       report(DEBUG,"EvtGen") << "Decay density matrix:"<<endl;
73       report(DEBUG,"EvtGen") << rho;
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           
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;
201           
202           _amp2.dump();
203
204           for(size_t ii=0;ii<i+1;ii++){
205             report(ERROR,"EvtGen") << "rho_list[" << ii << "] = " << rho_list[ii];
206           }
207
208           report(ERROR,"EvtGen") << "-------Done with error-------"<<endl;  
209
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