TEvtGen/EvtGenModels/EvtPropSLPole.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: EvtPropSLPole.cc
  12 //
  13 // Description: Routine to implement semileptonic decays according
  14 //              to light cone sum rules
  15 //
  16 // Modification history:
  17 //
  18 //    DJL       April 23, 1998       Module created
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include <stdlib.h>
  24 #include "EvtGenBase/EvtParticle.hh"
  25 #include "EvtGenBase/EvtGenKine.hh"
  26 #include "EvtGenBase/EvtPDL.hh"
  27 #include "EvtGenBase/EvtReport.hh"
  28 #include "EvtGenModels/EvtPropSLPole.hh"
  29 #include "EvtGenModels/EvtSLPoleFF.hh"
  30 #include "EvtGenBase/EvtSemiLeptonicScalarAmp.hh"
  31 #include "EvtGenBase/EvtSemiLeptonicVectorAmp.hh"
  32 #include "EvtGenBase/EvtSemiLeptonicTensorAmp.hh"
  33 #include "EvtGenBase/EvtIntervalFlatPdf.hh"
  34 #include "EvtGenBase/EvtScalarParticle.hh"
  35 #include "EvtGenBase/EvtVectorParticle.hh"
  36 #include "EvtGenBase/EvtTensorParticle.hh"
  37 #include "EvtGenBase/EvtTwoBodyVertex.hh"
  38 #include "EvtGenBase/EvtPropBreitWignerRel.hh"
  39 #include "EvtGenBase/EvtPDL.hh"
  40 #include "EvtGenBase/EvtAmpPdf.hh"
  41 #include "EvtGenBase/EvtMassAmp.hh"
  42 #include "EvtGenBase/EvtSpinType.hh"
  43 #include "EvtGenBase/EvtDecayTable.hh"
  44 #include <string>
  45
  46 EvtPropSLPole::~EvtPropSLPole() {}
  47
  48 std::string EvtPropSLPole::getName(){
  49
  50   return "PROPSLPOLE";
  51
  52 }
  53
  54
  55 EvtDecayBase* EvtPropSLPole::clone(){
  56
  57   return new EvtPropSLPole;
  58
  59 }
  60
  61 void EvtPropSLPole::decay( EvtParticle *p ){
  62
  63   if(! _isProbMaxSet){
  64
  65      EvtId parnum,mesnum,lnum,nunum;
  66
  67      parnum = getParentId();
  68      mesnum = getDaug(0);
  69      lnum = getDaug(1);
  70      nunum = getDaug(2);
  71
  72      double mymaxprob = calcMaxProb(parnum,mesnum,
  73                            lnum,nunum,SLPoleffmodel);
  74
  75      setProbMax(mymaxprob);
  76
  77      _isProbMaxSet = true;
  78
  79   }
  80
  81   double minKstMass = EvtPDL::getMinMass(p->getDaug(0)->getId());
  82   double maxKstMass = EvtPDL::getMaxMass(p->getDaug(0)->getId());
  83
  84   EvtIntervalFlatPdf flat(minKstMass, maxKstMass);
  85   EvtPdfGen<EvtPoint1D> gen(flat);
  86   EvtPoint1D point = gen();
  87
  88   double massKst = point.value();
  89
  90   p->getDaug(0)->setMass(massKst);
  91   p->initializePhaseSpace(getNDaug(),getDaugs());
  92
  93 //  EvtVector4R p4meson = p->getDaug(0)->getP4();
  94
  95   calcamp->CalcAmp(p,_amp2,SLPoleffmodel);
  96
  97   EvtParticle *mesonPart = p->getDaug(0);
  98
  99   double meson_BWAmp = calBreitWigner(mesonPart, point);
 100
 101   int list[2];
 102   list[0]=0; list[1]=0;
 103   _amp2.vertex(0,0,_amp2.getAmp(list)*meson_BWAmp);
 104   list[0]=0; list[1]=1;
 105   _amp2.vertex(0,1,_amp2.getAmp(list)*meson_BWAmp);
 106
 107   list[0]=1; list[1]=0;
 108   _amp2.vertex(1,0,_amp2.getAmp(list)*meson_BWAmp);
 109   list[0]=1; list[1]=1;
 110   _amp2.vertex(1,1,_amp2.getAmp(list)*meson_BWAmp);
 111
 112   list[0]=2; list[1]=0;
 113   _amp2.vertex(2,0,_amp2.getAmp(list)*meson_BWAmp);
 114   list[0]=2; list[1]=1;
 115   _amp2.vertex(2,1,_amp2.getAmp(list)*meson_BWAmp);
 116
 117
 118   return;
 119
 120 }
 121
 122 void EvtPropSLPole::initProbMax(){
 123
 124   _isProbMaxSet = false;
 125
 126   return;
 127
 128 }
 129
 130
 131 void EvtPropSLPole::init(){
 132
 133   checkNDaug(3);
 134
 135   //We expect the parent to be a scalar
 136   //and the daughters to be X lepton neutrino
 137
 138   checkSpinParent(EvtSpinType::SCALAR);
 139   checkSpinDaughter(1,EvtSpinType::DIRAC);
 140   checkSpinDaughter(2,EvtSpinType::NEUTRINO);
 141
 142   EvtSpinType::spintype mesontype=EvtPDL::getSpinType(getDaug(0));
 143
 144   SLPoleffmodel = new EvtSLPoleFF(getNArg(),getArgs());
 145
 146   if ( mesontype==EvtSpinType::SCALAR ) {
 147     calcamp = new EvtSemiLeptonicScalarAmp;
 148   }
 149   if ( mesontype==EvtSpinType::VECTOR ) {
 150     calcamp = new EvtSemiLeptonicVectorAmp;
 151   }
 152   if ( mesontype==EvtSpinType::TENSOR ) {
 153     calcamp = new EvtSemiLeptonicTensorAmp;
 154   }
 155
 156 }
 157
 158
 159 double EvtPropSLPole::calBreitWignerBasic(double maxMass){
 160
 161   if ( _width< 0.0001) return 1.0;
 162   //its not flat - but generated according to a BW
 163
 164   double mMin=_massMin;
 165   double mMax=_massMax;
 166   if ( maxMass>-0.5 && maxMass< mMax) mMax=maxMass;
 167
 168   double massGood = EvtRandom::Flat(mMin, mMax);
 169
 170   double ampVal = sqrt(1.0/(pow(massGood-_mass, 2.0) + pow(_width, 2.0)/4.0));
 171
 172   return ampVal;
 173
 174 }
 175
 176
 177 double EvtPropSLPole::calBreitWigner(EvtParticle *pmeson, EvtPoint1D point){
 178
 179   EvtId mesnum = pmeson->getId();
 180   double _mass = EvtPDL::getMeanMass(mesnum);
 181   double _width = EvtPDL::getWidth(mesnum);
 182   double _maxRange = EvtPDL::getMaxRange(mesnum);
 183   EvtSpinType::spintype mesontype=EvtPDL::getSpinType(mesnum);
 184   _includeDecayFact=true;
 185   _includeBirthFact=true;
 186   _spin = mesontype;
 187   _blatt = 3.0;
 188
 189   double maxdelta = 15.0*_width;
 190
 191   if ( _maxRange > 0.00001 ) {
 192     _massMax=_mass+maxdelta;
 193     _massMin=_mass-_maxRange;
 194   }
 195   else{
 196     _massMax=_mass+maxdelta;
 197     _massMin=_mass-15.0*_width;
 198   }
 199
 200   _massMax=_mass+maxdelta;
 201   if ( _massMin< 0. ) _massMin=0.;
 202
 203
 204   EvtParticle* par=pmeson->getParent();
 205   double maxMass=-1.;
 206   if ( par != 0 ) {
 207     if ( par->hasValidP4() ) maxMass=par->mass();
 208     for ( size_t i=0;i<par->getNDaug();i++) {
 209       EvtParticle *tDaug=par->getDaug(i);
 210       if ( pmeson != tDaug )
 211         maxMass-=EvtPDL::getMinMass(tDaug->getId());
 212     }
 213   }
 214
 215   EvtId *dauId=0;
 216   double *dauMasses=0;
 217   size_t nDaug = pmeson->getNDaug();
 218   if ( nDaug > 0) {
 219      dauId=new EvtId[nDaug];
 220      dauMasses=new double[nDaug];
 221      for (size_t j=0;j<nDaug;j++) {
 222        dauId[j]=pmeson->getDaug(j)->getId();
 223        dauMasses[j]=pmeson->getDaug(j)->mass();
 224      }
 225    }
 226    EvtId *parId=0;
 227    EvtId *othDaugId=0;
 228    EvtParticle *tempPar=pmeson->getParent();
 229    if (tempPar) {
 230      parId=new EvtId(tempPar->getId());
 231      if ( tempPar->getNDaug()==2 ) {
 232        if ( tempPar->getDaug(0) == pmeson ) othDaugId=new EvtId(tempPar->getDaug(1)->getId());
 233        else othDaugId=new EvtId(tempPar->getDaug(0)->getId());
 234      }
 235    }
 236
 237   if ( nDaug!=2) return calBreitWignerBasic(maxMass);
 238
 239   if ( _width< 0.00001) return 1.0;
 240
 241   //first figure out L - take the lowest allowed.
 242
 243   EvtSpinType::spintype spinD1=EvtPDL::getSpinType(dauId[0]);
 244   EvtSpinType::spintype spinD2=EvtPDL::getSpinType(dauId[1]);
 245
 246   int t1=EvtSpinType::getSpin2(spinD1);
 247   int t2=EvtSpinType::getSpin2(spinD2);
 248   int t3=EvtSpinType::getSpin2(_spin);
 249
 250   int Lmin=-10;
 251
 252   // allow for special cases.
 253   if (Lmin<-1 ) {
 254
 255     //There are some things I don't know how to deal with
 256     if ( t3>4) return calBreitWignerBasic(maxMass);
 257     if ( t1>4) return calBreitWignerBasic(maxMass);
 258     if ( t2>4) return calBreitWignerBasic(maxMass);
 259
 260     //figure the min and max allowwed "spins" for the daughters state
 261     Lmin=std::max(t3-t2-t1,std::max(t2-t3-t1,t1-t3-t2));
 262     if (Lmin<0) Lmin=0;
 263     assert(Lmin==0||Lmin==2||Lmin==4);
 264   }
 265
 266   //double massD1=EvtPDL::getMeanMass(dauId[0]);
 267   //double massD2=EvtPDL::getMeanMass(dauId[1]);
 268   double massD1=dauMasses[0];
 269   double massD2=dauMasses[1];
 270
 271   // I'm not sure how to define the vertex factor here - so retreat to nonRel code.
 272   if ( (massD1+massD2)> _mass ) return  calBreitWignerBasic(maxMass);
 273
 274   //parent vertex factor not yet implemented
 275   double massOthD=-10.;
 276   double massParent=-10.;
 277   int birthl=-10;
 278   if ( othDaugId) {
 279     EvtSpinType::spintype spinOth=EvtPDL::getSpinType(*othDaugId);
 280     EvtSpinType::spintype spinPar=EvtPDL::getSpinType(*parId);
 281
 282     int tt1=EvtSpinType::getSpin2(spinOth);
 283     int tt2=EvtSpinType::getSpin2(spinPar);
 284     int tt3=EvtSpinType::getSpin2(_spin);
 285
 286     //figure the min and max allowwed "spins" for the daughters state
 287     if ( (tt1<=4) && ( tt2<=4) ) {
 288       birthl=std::max(tt3-tt2-tt1,std::max(tt2-tt3-tt1,tt1-tt3-tt2));
 289       if (birthl<0) birthl=0;
 290
 291       massOthD=EvtPDL::getMeanMass(*othDaugId);
 292       massParent=EvtPDL::getMeanMass(*parId);
 293
 294     }
 295
 296   }
 297   double massM=_massMax;
 298   if ( (maxMass > -0.5) && (maxMass < massM) ) massM=maxMass;
 299
 300   //special case... if the parent mass is _fixed_ we can do a little better
 301   //and only for a two body decay as that seems to be where we have problems
 302
 303   // Define relativistic propagator amplitude
 304
 305   EvtTwoBodyVertex vd(massD1,massD2,_mass,Lmin/2);
 306   vd.set_f(_blatt);
 307   EvtPropBreitWignerRel bw(_mass,_width);
 308   EvtMassAmp amp(bw,vd);
 309 //  if ( _fixMassForMax) amp.fixUpMassForMax();
 310 //  else std::cout << "problem problem\n";
 311   if ( _includeDecayFact) {
 312     amp.addDeathFact();
 313     amp.addDeathFactFF();
 314   }
 315   if ( massParent>-1.) {
 316     if ( _includeBirthFact ) {
 317
 318       EvtTwoBodyVertex vb(_mass,massOthD,massParent,birthl/2);
 319       amp.setBirthVtx(vb);
 320       amp.addBirthFact();
 321       amp.addBirthFactFF();
 322     }
 323   }
 324
 325   EvtAmpPdf<EvtPoint1D> pdf(amp);
 326
 327   double ampVal = sqrt(pdf.evaluate(point));
 328
 329   if ( parId) delete parId;
 330   if ( othDaugId) delete othDaugId;
 331   if ( dauId) delete [] dauId;
 332   if ( dauMasses) delete [] dauMasses;
 333
 334   return ampVal;
 335
 336 }
 337
 338
 339 double EvtPropSLPole::calcMaxProb( EvtId parent, EvtId meson,
 340                                         EvtId lepton, EvtId nudaug,
 341                      EvtSemiLeptonicFF *FormFactors ) {
 342
 343   //This routine takes the arguements parent, meson, and lepton
 344   //number, and a form factor model, and returns a maximum
 345   //probability for this semileptonic form factor model.  A
 346   //brute force method is used.  The 2D cos theta lepton and
 347   //q2 phase space is probed.
 348
 349   //Start by declaring a particle at rest.
 350
 351   //It only makes sense to have a scalar parent.  For now.
 352   //This should be generalized later.
 353
 354   EvtScalarParticle *scalar_part;
 355   EvtParticle *root_part;
 356
 357   scalar_part=new EvtScalarParticle;
 358
 359   //cludge to avoid generating random numbers!
 360   scalar_part->noLifeTime();
 361
 362   EvtVector4R p_init;
 363
 364   p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0);
 365   scalar_part->init(parent,p_init);
 366   root_part=(EvtParticle *)scalar_part;
 367 //  root_part->set_type(EvtSpinType::SCALAR);
 368   root_part->setDiagonalSpinDensity();
 369
 370   EvtParticle *daughter, *lep, *trino;
 371
 372   EvtAmp amp;
 373
 374   EvtId listdaug[3];
 375   listdaug[0] = meson;
 376   listdaug[1] = lepton;
 377   listdaug[2] = nudaug;
 378
 379   amp.init(parent,3,listdaug);
 380
 381   root_part->makeDaughters(3,listdaug);
 382   daughter=root_part->getDaug(0);
 383   lep=root_part->getDaug(1);
 384   trino=root_part->getDaug(2);
 385
 386   EvtDecayBase *decayer;
 387   decayer = EvtDecayTable::getDecayFunc(daughter);
 388   if ( decayer ) {
 389     daughter->makeDaughters(decayer->nRealDaughters(),decayer->getDaugs());
 390     for(int ii=0; ii<decayer->nRealDaughters(); ii++){
 391       daughter->getDaug(ii)->setMass(EvtPDL::getMeanMass(daughter->getDaug(ii)->getId()));
 392     }
 393   }
 394
 395   //cludge to avoid generating random numbers!
 396   daughter->noLifeTime();
 397   lep->noLifeTime();
 398   trino->noLifeTime();
 399
 400   //Initial particle is unpolarized, well it is a scalar so it is
 401   //trivial
 402   EvtSpinDensity rho;
 403   rho.setDiag(root_part->getSpinStates());
 404
 405   double mass[3];
 406
 407   double m = root_part->mass();
 408
 409   EvtVector4R p4meson, p4lepton, p4nu, p4w;
 410   double q2max;
 411
 412   double q2, elepton, plepton;
 413   int i,j;
 414   double erho,prho,costl;
 415
 416   double maxfoundprob = 0.0;
 417   double prob = -10.0;
 418   int massiter;
 419
 420   for (massiter=0;massiter<3;massiter++){
 421
 422     mass[0] = EvtPDL::getMeanMass(meson);
 423     mass[1] = EvtPDL::getMeanMass(lepton);
 424     mass[2] = EvtPDL::getMeanMass(nudaug);
 425     if ( massiter==1 ) {
 426       mass[0] = EvtPDL::getMinMass(meson);
 427     }
 428     if ( massiter==2 ) {
 429       mass[0] = EvtPDL::getMaxMass(meson);
 430       if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001;
 431     }
 432
 433     q2max = (m-mass[0])*(m-mass[0]);
 434
 435     //loop over q2
 436
 437     for (i=0;i<25;i++) {
 438       q2 = ((i+0.5)*q2max)/25.0;
 439
 440       erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m);
 441
 442       prho = sqrt(erho*erho-mass[0]*mass[0]);
 443
 444       p4meson.set(erho,0.0,0.0,-1.0*prho);
 445       p4w.set(m-erho,0.0,0.0,prho);
 446
 447       //This is in the W rest frame
 448       elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2));
 449       plepton = sqrt(elepton*elepton-mass[1]*mass[1]);
 450
 451       double probctl[3];
 452
 453       for (j=0;j<3;j++) {
 454
 455         costl = 0.99*(j - 1.0);
 456
 457         //These are in the W rest frame. Need to boost out into
 458         //the B frame.
 459         p4lepton.set(elepton,0.0,
 460                   plepton*sqrt(1.0-costl*costl),plepton*costl);
 461         p4nu.set(plepton,0.0,
 462                  -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);
 463
 464         EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
 465         p4lepton=boostTo(p4lepton,boost);
 466         p4nu=boostTo(p4nu,boost);
 467
 468         //Now initialize the daughters...
 469
 470         daughter->init(meson,p4meson);
 471         lep->init(lepton,p4lepton);
 472         trino->init(nudaug,p4nu);
 473
 474         calcamp->CalcAmp(root_part,amp,FormFactors);
 475
 476         EvtPoint1D *point = new EvtPoint1D(mass[0]);
 477
 478         double meson_BWAmp = calBreitWigner(daughter, *point);
 479
 480         int list[2];
 481         list[0]=0; list[1]=0;
 482         amp.vertex(0,0,amp.getAmp(list)*meson_BWAmp);
 483         list[0]=0; list[1]=1;
 484         amp.vertex(0,1,amp.getAmp(list)*meson_BWAmp);
 485
 486         list[0]=1; list[1]=0;
 487         amp.vertex(1,0,amp.getAmp(list)*meson_BWAmp);
 488         list[0]=1; list[1]=1;
 489         amp.vertex(1,1,amp.getAmp(list)*meson_BWAmp);
 490
 491         list[0]=2; list[1]=0;
 492         amp.vertex(2,0,amp.getAmp(list)*meson_BWAmp);
 493         list[0]=2; list[1]=1;
 494         amp.vertex(2,1,amp.getAmp(list)*meson_BWAmp);
 495
 496         //Now find the probability at this q2 and cos theta lepton point
 497         //and compare to maxfoundprob.
 498
 499         //Do a little magic to get the probability!!
 500         prob = rho.normalizedProb(amp.getSpinDensity());
 501
 502         probctl[j]=prob;
 503       }
 504
 505       //probclt contains prob at ctl=-1,0,1.
 506       //prob=a+b*ctl+c*ctl^2
 507
 508       double a=probctl[1];
 509       double b=0.5*(probctl[2]-probctl[0]);
 510       double c=0.5*(probctl[2]+probctl[0])-probctl[1];
 511
 512       prob=probctl[0];
 513       if (probctl[1]>prob) prob=probctl[1];
 514       if (probctl[2]>prob) prob=probctl[2];
 515
 516       if (fabs(c)>1e-20){
 517         double ctlx=-0.5*b/c;
 518         if (fabs(ctlx)<1.0){
 519           double probtmp=a+b*ctlx+c*ctlx*ctlx;
 520           if (probtmp>prob) prob=probtmp;
 521         }
 522
 523       }
 524
 525       //report(DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" "
 526       //                            << probctl[0]<<" "
 527       //                            << probctl[1]<<" "
 528       //                            << probctl[2]<<endl;
 529
 530       if ( prob > maxfoundprob ) {
 531         maxfoundprob = prob;
 532       }
 533
 534     }
 535     if ( EvtPDL::getWidth(meson) <= 0.0 ) {
 536       //if the particle is narrow dont bother with changing the mass.
 537       massiter = 4;
 538     }
 539
 540   }
 541   root_part->deleteTree();
 542
 543   maxfoundprob *=1.1;
 544   return maxfoundprob;
 545
 546 }
 547