TEvtGen/EvtGenBase/EvtGenKine.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: EvtGenKine.cc
  12 //
  13 // Description: Tools for generating distributions of four vectors in
  14 //              phasespace
  15 //
  16 // Modification history:
  17 //
  18 //    DJL/RYD     September 25, 1996         Module created
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include <iostream>
  24 #include "EvtGenBase/EvtGenKine.hh"
  25 #include "EvtGenBase/EvtRandom.hh"
  26 #include "EvtGenBase/EvtVector4R.hh"
  27 #include "EvtGenBase/EvtReport.hh"
  28 #include "EvtGenBase/EvtConst.hh"
  29 #include <math.h>
  30 using std::endl;
  31
  32
  33 double EvtPawt(double a,double b,double c)
  34 {
  35   double temp=(a*a-(b+c)*(b+c))*(a*a-(b-c)*(b-c));
  36
  37   if (temp<=0) {
  38      return 0.0;
  39   }
  40
  41   return sqrt(temp)/(2.0*a);
  42 }
  43
  44
  45 double EvtGenKine::PhaseSpace( int ndaug, double mass[30], EvtVector4R p4[30],
  46                                double mp )
  47
  48 //  N body phase space routine.  Send parent with
  49 //  daughters already defined ( Number and masses )
  50 //  Returns four vectors in parent frame.
  51
  52 {
  53
  54   double energy, p3, alpha, beta;
  55
  56   if ( ndaug == 1 ) {
  57      p4[0].set(mass[0],0.0,0.0,0.0);
  58      return 1.0;
  59   }
  60
  61
  62   if ( ndaug == 2 ) {
  63
  64      //Two body phase space
  65
  66      energy = ( mp*mp + mass[0]*mass[0] -
  67               mass[1]*mass[1] ) / ( 2.0 * mp );
  68
  69      p3 = sqrt( energy*energy - mass[0]*mass[0] );
  70
  71      p4[0].set( energy, 0.0, 0.0, p3 );
  72
  73      energy = mp - energy;
  74      p3 = -1.0*p3;
  75      p4[1].set( energy, 0.0, 0.0, p3 );
  76
  77      //Now rotate four vectors.
  78
  79      alpha = EvtRandom::Flat( EvtConst::twoPi );
  80      beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
  81
  82      p4[0].applyRotateEuler( alpha, beta, -alpha );
  83      p4[1].applyRotateEuler( alpha, beta, -alpha );
  84
  85      return 1.0;
  86   }
  87
  88   if ( ndaug != 2 ) {
  89
  90     double wtmax=0.0;
  91     double pm[5][30],to[4],pmin,pmax,psum,rnd[30];
  92     double ran,wt,pa,costh,sinth,phi,p[4][30],be[4],bep,temp;
  93     int i,il,ilr,i1,il1u,il1,il2r,ilu;
  94     int il2=0;
  95
  96      for(i=0;i<ndaug;i++){
  97        pm[4][i]=0.0;
  98        rnd[i]=0.0;
  99      }
 100
 101      pm[0][0]=mp;
 102      pm[1][0]=0.0;
 103      pm[2][0]=0.0;
 104      pm[3][0]=0.0;
 105      pm[4][0]=mp;
 106
 107      to[0]=mp;
 108      to[1]=0.0;
 109      to[2]=0.0;
 110      to[3]=0.0;
 111
 112      psum=0.0;
 113      for(i=1;i<ndaug+1;i++){
 114        psum=psum+mass[i-1];
 115      }
 116
 117      pm[4][ndaug-1]=mass[ndaug-1];
 118
 119      switch (ndaug) {
 120
 121      case 1:
 122        wtmax=1.0/16.0;
 123        break;
 124      case 2:
 125        wtmax=1.0/150.0;
 126        break;
 127      case 3:
 128        wtmax=1.0/2.0;
 129        break;
 130      case 4:
 131        wtmax=1.0/5.0;
 132        break;
 133      case 5:
 134        wtmax=1.0/15.0;
 135        break;
 136      case 6:
 137        wtmax=1.0/15.0;
 138        break;
 139      case 7:
 140        wtmax=1.0/15.0;
 141        break;
 142      case 8:
 143        wtmax=1.0/15.0;
 144        break;
 145      case 9:
 146        wtmax=1.0/15.0;
 147        break;
 148      case 10:
 149        wtmax=1.0/15.0;
 150        break;
 151      case 11:
 152        wtmax=1.0/15.0;
 153        break;
 154      case 12:
 155        wtmax=1.0/15.0;
 156        break;
 157      case 13:
 158        wtmax=1.0/15.0;
 159        break;
 160      case 14:
 161        wtmax=1.0/15.0;
 162        break;
 163      case 15:
 164        wtmax=1.0/15.0;
 165        break;
 166      default:
 167        report(ERROR,"EvtGen") << "too many daughters for phase space..." << ndaug << " "<< mp <<endl;;
 168        break;
 169      }
 170
 171      pmax=mp-psum+mass[ndaug-1];
 172
 173      pmin=0.0;
 174
 175      for(ilr=2;ilr<ndaug+1;ilr++){
 176        il=ndaug+1-ilr;
 177        pmax=pmax+mass[il-1];
 178        pmin=pmin+mass[il+1-1];
 179        wtmax=wtmax*EvtPawt(pmax,pmin,mass[il-1]);
 180      }
 181
 182      do{
 183
 184        rnd[0]=1.0;
 185        il1u=ndaug-1;
 186
 187        for (il1=2;il1<il1u+1;il1++){
 188          ran=EvtRandom::Flat();
 189          for (il2r=2;il2r<il1+1;il2r++){
 190            il2=il1+1-il2r;
 191            if (ran<=rnd[il2-1]) goto two39;
 192            rnd[il2+1-1]=rnd[il2-1];
 193          }
 194        two39:
 195          rnd[il2+1-1]=ran;
 196        }
 197
 198        rnd[ndaug-1]=0.0;
 199        wt=1.0;
 200        for(ilr=2;ilr<ndaug+1;ilr++){
 201          il=ndaug+1-ilr;
 202          pm[4][il-1]=pm[4][il+1-1]+mass[il-1]+
 203            (rnd[il-1]-rnd[il+1-1])*(mp-psum);
 204          wt=wt*EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
 205        }
 206        if (wt>wtmax) {
 207          report(ERROR,"EvtGen") << "wtmax to small in EvtPhaseSpace with "
 208                                 << ndaug <<" daughters"<<endl;;
 209        }
 210      } while (wt<EvtRandom::Flat(wtmax));
 211
 212      ilu=ndaug-1;
 213
 214      for (il=1;il<ilu+1;il++){
 215        pa=EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
 216        costh=EvtRandom::Flat(-1.0,1.0);
 217        sinth=sqrt(1.0-costh*costh);
 218        phi=EvtRandom::Flat(EvtConst::twoPi);
 219        p[1][il-1]=pa*sinth*cos(phi);
 220        p[2][il-1]=pa*sinth*sin(phi);
 221        p[3][il-1]=pa*costh;
 222        pm[1][il+1-1]=-p[1][il-1];
 223        pm[2][il+1-1]=-p[2][il-1];
 224        pm[3][il+1-1]=-p[3][il-1];
 225        p[0][il-1]=sqrt(pa*pa+mass[il-1]*mass[il-1]);
 226        pm[0][il+1-1]=sqrt(pa*pa+pm[4][il+1-1]*pm[4][il+1-1]);
 227      }
 228
 229       p[0][ndaug-1]=pm[0][ndaug-1];
 230       p[1][ndaug-1]=pm[1][ndaug-1];
 231       p[2][ndaug-1]=pm[2][ndaug-1];
 232       p[3][ndaug-1]=pm[3][ndaug-1];
 233
 234       for (ilr=2;ilr<ndaug+1;ilr++){
 235         il=ndaug+1-ilr;
 236         be[0]=pm[0][il-1]/pm[4][il-1];
 237         be[1]=pm[1][il-1]/pm[4][il-1];
 238         be[2]=pm[2][il-1]/pm[4][il-1];
 239         be[3]=pm[3][il-1]/pm[4][il-1];
 240
 241         for (i1=il;i1<ndaug+1;i1++){
 242           bep=be[1]*p[1][i1-1]+be[2]*p[2][i1-1]+
 243               be[3]*p[3][i1-1]+be[0]*p[0][i1-1];
 244           temp=(p[0][i1-1]+bep)/(be[0]+1.0);
 245           p[1][i1-1]=p[1][i1-1]+temp*be[1];
 246           p[2][i1-1]=p[2][i1-1]+temp*be[2];
 247           p[3][i1-1]=p[3][i1-1]+temp*be[3];
 248           p[0][i1-1]=bep;
 249
 250         }
 251       }
 252
 253      for (ilr=0;ilr<ndaug;ilr++){
 254        p4[ilr].set(p[0][ilr],p[1][ilr],p[2][ilr],p[3][ilr]);
 255      }
 256
 257      return 1.0;
 258   }
 259
 260   return 1.0;
 261
 262 }
 263
 264
 265 double EvtGenKine::PhaseSpacePole(double M, double m1, double m2, double m3,
 266                                   double a,EvtVector4R p4[10])
 267
 268 //  generate kinematics for 3 body decays, pole for the m1,m2 mass.
 269
 270 {
 271
 272   //f1   = 1  (phasespace)
 273   //f2   = a*(1/m12sq)^2
 274
 275   double m12sqmax=(M-m3)*(M-m3);
 276   double m12sqmin=(m1+m2)*(m1+m2);
 277
 278   double m13sqmax=(M-m2)*(M-m2);
 279   double m13sqmin=(m1+m3)*(m1+m3);
 280
 281   double v1=(m12sqmax-m12sqmin)*(m13sqmax-m13sqmin);
 282   double v2= a*(1.0/m12sqmin-1.0/m12sqmax)*(m13sqmax-m13sqmin);
 283
 284   double m12sq,m13sq;
 285
 286   double r=v1/(v1+v2);
 287
 288   double m13min,m13max;
 289
 290   do{
 291
 292     m13sq=EvtRandom::Flat(m13sqmin,m13sqmax);
 293
 294     if (r>EvtRandom::Flat()){
 295       m12sq=EvtRandom::Flat(m12sqmin,m12sqmax);
 296     }
 297     else{
 298       m12sq=1.0/(1.0/m12sqmin-EvtRandom::Flat()*(1.0/m12sqmin-1.0/m12sqmax));
 299     }
 300
 301     //kinematically allowed?
 302     double E3star=(M*M-m12sq-m3*m3)/sqrt(4*m12sq);
 303     double E1star=(m12sq+m1*m1-m2*m2)/sqrt(4*m12sq);
 304     double p3star=sqrt(E3star*E3star-m3*m3);
 305     double p1star=sqrt(E1star*E1star-m1*m1);
 306     m13max=(E3star+E1star)*(E3star+E1star)-
 307       (p3star-p1star)*(p3star-p1star);
 308     m13min=(E3star+E1star)*(E3star+E1star)-
 309       (p3star+p1star)*(p3star+p1star);
 310
 311   }while(m13sq<m13min||m13sq>m13max);
 312
 313
 314   double E2=(M*M+m2*m2-m13sq)/(2.0*M);
 315   double E3=(M*M+m3*m3-m12sq)/(2.0*M);
 316   double E1=M-E2-E3;
 317   double p1mom=sqrt(E1*E1-m1*m1);
 318   double p3mom=sqrt(E3*E3-m3*m3);
 319   double cost13=(2.0*E1*E3+m1*m1+m3*m3-m13sq)/(2.0*p1mom*p3mom);
 320
 321   //report(INFO,"EvtGen") << m13sq << endl;
 322   //report(INFO,"EvtGen") << m12sq << endl;
 323   //report(INFO,"EvtGen") << E1 << endl;
 324   //report(INFO,"EvtGen") << E2 << endl;
 325   //report(INFO,"EvtGen") << E3 << endl;
 326   //report(INFO,"EvtGen") << p1mom << endl;
 327   //report(INFO,"EvtGen") << p3mom << endl;
 328   //report(INFO,"EvtGen") << cost13 << endl;
 329
 330
 331   p4[2].set(E3,0.0,0.0,p3mom);
 332   p4[0].set(E1,p1mom*sqrt(1.0-cost13*cost13),0.0,p1mom*cost13);
 333   p4[1].set(E2,-p1mom*sqrt(1.0-cost13*cost13),0.0,-p1mom*cost13-p3mom);
 334
 335   //report(INFO,"EvtGen") << "p4:"<<p4[0]<<p4[1]<<p4[2]<<endl;
 336
 337   double alpha = EvtRandom::Flat( EvtConst::twoPi );
 338   double beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
 339   double gamma = EvtRandom::Flat( EvtConst::twoPi );
 340
 341   p4[0].applyRotateEuler( alpha, beta, gamma );
 342   p4[1].applyRotateEuler( alpha, beta, gamma );
 343   p4[2].applyRotateEuler( alpha, beta, gamma );
 344
 345   return 1.0+a/(m12sq*m12sq);
 346
 347 }
 348