TEvtGen/EvtGenModels/EvtPartWave.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) 2000      Caltech, UCSB
  10 //
  11 // Module: EvtHelAmp.cc
  12 //
  13 // Description: Decay model for implementation of generic 2 body
  14 //              decay specified by the partial wave amplitudes
  15 //
  16 //
  17 // Modification history:
  18 //
  19 //    fkw        February 2, 2001     changes to satisfy KCC
  20 //    RYD       September 7, 2000       Module created
  21 //
  22 //------------------------------------------------------------------------
  23 //
  24 #include "EvtGenBase/EvtPatches.hh"
  25 #include <stdlib.h>
  26 #include "EvtGenBase/EvtParticle.hh"
  27 #include "EvtGenBase/EvtGenKine.hh"
  28 #include "EvtGenBase/EvtPDL.hh"
  29 #include "EvtGenBase/EvtVector4C.hh"
  30 #include "EvtGenBase/EvtTensor4C.hh"
  31 #include "EvtGenBase/EvtReport.hh"
  32 #include "EvtGenModels/EvtPartWave.hh"
  33 #include "EvtGenBase/EvtEvalHelAmp.hh"
  34 #include "EvtGenBase/EvtId.hh"
  35 #include <string>
  36 #include "EvtGenBase/EvtConst.hh"
  37 #include "EvtGenBase/EvtKine.hh"
  38 #include "EvtGenBase/EvtCGCoefSingle.hh"
  39 #include <algorithm>
  40 using std::endl;
  41 EvtPartWave::~EvtPartWave() {}
  42
  43 std::string EvtPartWave::getName(){
  44
  45   return "PARTWAVE";
  46
  47 }
  48
  49
  50 EvtDecayBase* EvtPartWave::clone(){
  51
  52   return new EvtPartWave;
  53
  54 }
  55
  56 void EvtPartWave::init(){
  57
  58   checkNDaug(2);
  59
  60   //find out how many states each particle have
  61   int _nA=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getParentId()));
  62   int _nB=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getDaug(0)));
  63   int _nC=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getDaug(1)));
  64
  65   if (verbose()){
  66     report(INFO,"EvtGen")<<"_nA,_nB,_nC:"
  67                          <<_nA<<","<<_nB<<","<<_nC<<endl;
  68   }
  69
  70   //find out what 2 times the spin is
  71   int _JA2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getParentId()));
  72   int _JB2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(0)));
  73   int _JC2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(1)));
  74
  75   if (verbose()){
  76     report(INFO,"EvtGen")<<"_JA2,_JB2,_JC2:"
  77                          <<_JA2<<","<<_JB2<<","<<_JC2<<endl;
  78   }
  79
  80
  81   //allocate memory
  82   int* _lambdaA2=new int[_nA];
  83   int* _lambdaB2=new int[_nB];
  84   int* _lambdaC2=new int[_nC];
  85
  86   EvtComplexPtr* _HBC=new EvtComplexPtr[_nB];
  87   int ib,ic;
  88   for(ib=0;ib<_nB;ib++){
  89     _HBC[ib]=new EvtComplex[_nC];
  90   }
  91
  92
  93   int i;
  94   //find the allowed helicities (actually 2*times the helicity!)
  95
  96   fillHelicity(_lambdaA2,_nA,_JA2);
  97   fillHelicity(_lambdaB2,_nB,_JB2);
  98   fillHelicity(_lambdaC2,_nC,_JC2);
  99
 100   if (verbose()){
 101     report(INFO,"EvtGen")<<"Helicity states of particle A:"<<endl;
 102     for(i=0;i<_nA;i++){
 103       report(INFO,"EvtGen")<<_lambdaA2[i]<<endl;
 104     }
 105
 106     report(INFO,"EvtGen")<<"Helicity states of particle B:"<<endl;
 107     for(i=0;i<_nB;i++){
 108       report(INFO,"EvtGen")<<_lambdaB2[i]<<endl;
 109     }
 110
 111     report(INFO,"EvtGen")<<"Helicity states of particle C:"<<endl;
 112     for(i=0;i<_nC;i++){
 113       report(INFO,"EvtGen")<<_lambdaC2[i]<<endl;
 114     }
 115
 116     report(INFO,"EvtGen")<<"Will now figure out the valid (M_LS) states:"<<endl;
 117
 118   }
 119
 120   int Lmin=std::max(_JA2-_JB2-_JC2,std::max(_JB2-_JA2-_JC2,_JC2-_JA2-_JB2));
 121   if (Lmin<0) Lmin=0;
 122   //int Lmin=_JA2-_JB2-_JC2;
 123   int Lmax=_JA2+_JB2+_JC2;
 124
 125   int L;
 126
 127   int _nPartialWaveAmp=0;
 128
 129   int _nL[50];
 130   int _nS[50];
 131
 132   for (L=Lmin;L<=Lmax;L+=2){
 133     int Smin=abs(L-_JA2);
 134     if (Smin<abs(_JB2-_JC2)) Smin=abs(_JB2-_JC2);
 135     int Smax=L+_JA2;
 136     if (Smax>abs(_JB2+_JC2)) Smax=abs(_JB2+_JC2);
 137     int S;
 138     for (S=Smin;S<=Smax;S+=2){
 139       _nL[_nPartialWaveAmp]=L;
 140       _nS[_nPartialWaveAmp]=S;
 141
 142       _nPartialWaveAmp++;
 143       if (verbose()){
 144         report(INFO,"EvtGen")<<"M["<<L<<"]["<<S<<"]"<<endl;
 145       }
 146     }
 147   }
 148
 149   checkNArg(_nPartialWaveAmp*2);
 150
 151   int argcounter=0;
 152
 153   EvtComplex _M[50];
 154
 155   double partampsqtot=0.0;
 156
 157   for(i=0;i<_nPartialWaveAmp;i++){
 158     _M[i]=getArg(argcounter)*exp(EvtComplex(0.0,getArg(argcounter+1)));;
 159     argcounter+=2;
 160     partampsqtot+=abs2(_M[i]);
 161     if (verbose()){
 162       report(INFO,"EvtGen")<<"M["<<_nL[i]<<"]["<<_nS[i]<<"]="<<_M[i]<<endl;
 163     }
 164   }
 165
 166   //Now calculate the helicity amplitudes
 167
 168   double helampsqtot=0.0;
 169
 170   for(ib=0;ib<_nB;ib++){
 171     for(ic=0;ic<_nC;ic++){
 172       _HBC[ib][ic]=0.0;
 173       if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2){
 174         for(i=0;i<_nPartialWaveAmp;i++){
 175           int L=_nL[i];
 176           int S=_nS[i];
 177           int lambda2=_lambdaB2[ib];
 178           int lambda3=_lambdaC2[ic];
 179           int s1=_JA2;
 180           int s2=_JB2;
 181           int s3=_JC2;
 182           int m1=lambda2-lambda3;
 183           EvtCGCoefSingle c1(s2,s3);
 184           EvtCGCoefSingle c2(L,S);
 185
 186           if (verbose()){
 187             report(INFO,"EvtGen") << "s2,lambda2:"<<s2<<" "<<lambda2<<endl;
 188           }
 189           //fkw changes to satisfy KCC
 190           double fkwTmp = (L+1.0)/(s1+1.0);
 191
 192           if (S>=abs(m1)){
 193
 194             EvtComplex tmp=sqrt(fkwTmp)
 195               *c1.coef(S,m1,s2,s3,lambda2,-lambda3)
 196               *c2.coef(s1,m1,L,S,0,m1)*_M[i];
 197             _HBC[ib][ic]+=tmp;
 198           }
 199         }
 200         if (verbose()){
 201           report(INFO,"EvtGen")<<"_HBC["<<ib<<"]["<<ic<<"]="<<_HBC[ib][ic]<<endl;
 202         }
 203       }
 204       helampsqtot+=abs2(_HBC[ib][ic]);
 205     }
 206   }
 207
 208   if (fabs(helampsqtot-partampsqtot)/(helampsqtot+partampsqtot)>1e-6){
 209       report(ERROR,"EvtGen")<<"In EvtPartWave for decay "
 210                             << EvtPDL::name(getParentId()) << " -> "
 211                             << EvtPDL::name(getDaug(0)) << " "
 212                             << EvtPDL::name(getDaug(1)) << std::endl;
 213       report(ERROR,"EvtGen")<<"With arguments: "<<std::endl;
 214       for(i=0;i*2<getNArg();i++){
 215           report(ERROR,"EvtGen") <<"M("<<_nL[i]<<","<<_nS[i]<<")="
 216 //                               <<getArg(2*i)<<" "<<getArg(2*i+1)<<std::endl;
 217                                  <<_M[i]<<std::endl;
 218       }
 219       report(ERROR,"EvtGen")<< "The total probability in the partwave basis is: "
 220                             << partampsqtot << std::endl;
 221       report(ERROR,"EvtGen")<< "The total probability in the helamp basis is: "
 222                             << helampsqtot << std::endl;
 223       report(ERROR,"EvtGen")<< "Most likely this is because the specified partwave amplitudes "
 224                             << std::endl;
 225       report(ERROR,"EvtGen")<< "project onto unphysical helicities of photons or neutrinos. "
 226                             << std::endl;
 227       report(ERROR,"EvtGen")<< "Seriously consider if your specified amplitudes are correct. "
 228                             << std::endl;
 229
 230
 231   }
 232
 233   _evalHelAmp=new EvtEvalHelAmp(getParentId(),
 234                                 getDaug(0),
 235                                 getDaug(1),
 236                                 _HBC);
 237
 238 }
 239
 240
 241 void EvtPartWave::initProbMax(){
 242
 243   double maxprob=_evalHelAmp->probMax();
 244
 245   if (verbose()){
 246     report(INFO,"EvtGen")<<"Calculated probmax"<<maxprob<<endl;
 247   }
 248
 249   setProbMax(maxprob);
 250
 251 }
 252
 253
 254 void EvtPartWave::decay( EvtParticle *p){
 255
 256   //first generate simple phase space
 257   p->initializePhaseSpace(getNDaug(),getDaugs());
 258
 259   _evalHelAmp->evalAmp(p,_amp2);
 260
 261   return;
 262
 263 }
 264
 265
 266
 267 void EvtPartWave::fillHelicity(int* lambda2,int n,int J2){
 268
 269   int i;
 270
 271   //photon is special case!
 272   if (n==2&&J2==2) {
 273     lambda2[0]=2;
 274     lambda2[1]=-2;
 275     return;
 276   }
 277
 278   assert(n==J2+1);
 279
 280   for(i=0;i<n;i++){
 281     lambda2[i]=n-i*2-1;
 282   }
 283
 284   return;
 285
 286 }
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297