TEvtGen/EvtGenBase/EvtDecayAmp.cxx

   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(ERROR,"EvtGen") << p->getSpinDensityForward();
  71
  72       report(DEBUG,"EvtGen") << "Decay density matrix:"<<endl;
  73       report(ERROR,"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           report(ERROR,"EvtGen")<<"Parent:"<<EvtPDL::name(p->getParent()->getId()).c_str()<<endl;
 183           report(ERROR,"EvtGen")<<"GrandParent:"<<EvtPDL::name(p->getParent()->getParent()->getId()).c_str()<<endl;
 184           report(ERROR,"EvtGen")<<"GrandGrandParent:"<<EvtPDL::name(p->getParent()->getParent()->getParent()->getId()).c_str()<<endl;
 185
 186           report(ERROR,"EvtGen") << rho;
 187
 188           _amp2.dump();
 189           for(size_t ii=0;ii<i+1;ii++){
 190             report(ERROR,"EvtGen") << rho_list[ii];
 191           }
 192           report(ERROR,"EvtGen") << "-------Done with error-------"<<endl;
 193         }
 194
 195         p->getDaug(i)->setSpinDensityForward(rho);
 196         p->getDaug(i)->decay();
 197
 198         rho_list[i+1]=p->getDaug(i)->getSpinDensityBackward();
 199
 200         if (_amp2.dstates[i]!=1){
 201         ampcont=ampcont.contract(_amp2._dnontrivial[i],rho_list[i+1]);
 202         }
 203
 204
 205       }
 206
 207       p->setSpinDensityBackward(_amp2.getBackwardSpinDensity(rho_list));
 208
 209
 210       if (!p->getSpinDensityBackward().check()) {
 211
 212         report(ERROR,"EvtGen")<<"rho_backward failed Check"<<
 213           p->getId().getId()<<" "<<p->getChannel()<<endl;
 214
 215         report(ERROR,"EvtGen") << p->getSpinDensityBackward();
 216
 217       }
 218     }
 219   }
 220
 221
 222   if (getPHOTOS() || EvtRadCorr::alwaysRadCorr()) {
 223     int n_daug_orig=p->getNDaug();
 224     EvtRadCorr::doRadCorr(p);
 225     int n_daug_new=p->getNDaug();
 226     for (int i=n_daug_orig;i<n_daug_new;i++){
 227       p->getDaug(i)->decay();
 228     }
 229   }
 230
 231 }
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241