diff --git a/src/EvtGenModels/EvtPropSLPole.cpp b/src/EvtGenModels/EvtPropSLPole.cpp index 7e1f02a..9204ff8 100644 --- a/src/EvtGenModels/EvtPropSLPole.cpp +++ b/src/EvtGenModels/EvtPropSLPole.cpp @@ -1,549 +1,549 @@ //-------------------------------------------------------------------------- // // Environment: // This software is part of the EvtGen package developed jointly // for the BaBar and CLEO collaborations. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: See EvtGen/COPYRIGHT // Copyright (C) 1998 Caltech, UCSB // // Module: EvtPropSLPole.cc // // Description: Routine to implement semileptonic decays according // to light cone sum rules // // Modification history: // // DJL April 23, 1998 Module created // //------------------------------------------------------------------------ // #include "EvtGenBase/EvtPatches.hh" #include #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtGenKine.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenModels/EvtPropSLPole.hh" #include "EvtGenModels/EvtSLPoleFF.hh" #include "EvtGenBase/EvtSemiLeptonicScalarAmp.hh" #include "EvtGenBase/EvtSemiLeptonicVectorAmp.hh" #include "EvtGenBase/EvtSemiLeptonicTensorAmp.hh" #include "EvtGenBase/EvtIntervalFlatPdf.hh" #include "EvtGenBase/EvtScalarParticle.hh" #include "EvtGenBase/EvtVectorParticle.hh" #include "EvtGenBase/EvtTensorParticle.hh" #include "EvtGenBase/EvtTwoBodyVertex.hh" #include "EvtGenBase/EvtPropBreitWignerRel.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtAmpPdf.hh" #include "EvtGenBase/EvtMassAmp.hh" #include "EvtGenBase/EvtSpinType.hh" #include "EvtGenBase/EvtDecayTable.hh" #include std::string EvtPropSLPole::getName(){ return "PROPSLPOLE"; } EvtDecayBase* EvtPropSLPole::clone(){ return new EvtPropSLPole; } void EvtPropSLPole::decay( EvtParticle *p ){ if(! _isProbMaxSet){ EvtId parnum,mesnum,lnum,nunum; parnum = getParentId(); mesnum = getDaug(0); lnum = getDaug(1); nunum = getDaug(2); double mymaxprob = calcMaxProb(parnum,mesnum, lnum,nunum,SLPoleffmodel.get()); setProbMax(mymaxprob); _isProbMaxSet = true; } double minKstMass = EvtPDL::getMinMass(p->getDaug(0)->getId()); double maxKstMass = EvtPDL::getMaxMass(p->getDaug(0)->getId()); EvtIntervalFlatPdf flat(minKstMass, maxKstMass); EvtPdfGen gen(flat); EvtPoint1D point = gen(); double massKst = point.value(); p->getDaug(0)->setMass(massKst); p->initializePhaseSpace(getNDaug(),getDaugs()); // EvtVector4R p4meson = p->getDaug(0)->getP4(); calcamp->CalcAmp(p,_amp2,SLPoleffmodel.get()); EvtParticle *mesonPart = p->getDaug(0); double meson_BWAmp = calBreitWigner(mesonPart, point); int list[2]; list[0]=0; list[1]=0; _amp2.vertex(0,0,_amp2.getAmp(list)*meson_BWAmp); list[0]=0; list[1]=1; _amp2.vertex(0,1,_amp2.getAmp(list)*meson_BWAmp); list[0]=1; list[1]=0; _amp2.vertex(1,0,_amp2.getAmp(list)*meson_BWAmp); list[0]=1; list[1]=1; _amp2.vertex(1,1,_amp2.getAmp(list)*meson_BWAmp); list[0]=2; list[1]=0; _amp2.vertex(2,0,_amp2.getAmp(list)*meson_BWAmp); list[0]=2; list[1]=1; _amp2.vertex(2,1,_amp2.getAmp(list)*meson_BWAmp); return; } void EvtPropSLPole::initProbMax(){ _isProbMaxSet = false; return; } void EvtPropSLPole::init(){ checkNDaug(3); //We expect the parent to be a scalar //and the daughters to be X lepton neutrino checkSpinParent(EvtSpinType::SCALAR); checkSpinDaughter(1,EvtSpinType::DIRAC); checkSpinDaughter(2,EvtSpinType::NEUTRINO); EvtSpinType::spintype mesontype=EvtPDL::getSpinType(getDaug(0)); SLPoleffmodel = std::make_unique(getNArg(),getArgs()); switch(mesontype) { case EvtSpinType::SCALAR: calcamp = std::make_unique(); break; case EvtSpinType::VECTOR: calcamp = std::make_unique(); break; case EvtSpinType::TENSOR: calcamp = std::make_unique(); break; default: ; } } double EvtPropSLPole::calBreitWignerBasic(double maxMass){ if ( _width< 0.0001) return 1.0; //its not flat - but generated according to a BW double mMin=_massMin; double mMax=_massMax; if ( maxMass>-0.5 && maxMass< mMax) mMax=maxMass; double massGood = EvtRandom::Flat(mMin, mMax); double ampVal = sqrt(1.0/(pow(massGood-_mass, 2.0) + pow(_width, 2.0)/4.0)); return ampVal; } double EvtPropSLPole::calBreitWigner(EvtParticle *pmeson, EvtPoint1D point){ EvtId mesnum = pmeson->getId(); - double _mass = EvtPDL::getMeanMass(mesnum); - double _width = EvtPDL::getWidth(mesnum); - double _maxRange = EvtPDL::getMaxRange(mesnum); + _mass = EvtPDL::getMeanMass(mesnum); + _width = EvtPDL::getWidth(mesnum); + _maxRange = EvtPDL::getMaxRange(mesnum); EvtSpinType::spintype mesontype=EvtPDL::getSpinType(mesnum); _includeDecayFact=true; _includeBirthFact=true; _spin = mesontype; _blatt = 3.0; double maxdelta = 15.0*_width; if ( _maxRange > 0.00001 ) { _massMax=_mass+maxdelta; _massMin=_mass-_maxRange; } else{ _massMax=_mass+maxdelta; _massMin=_mass-15.0*_width; } _massMax=_mass+maxdelta; if ( _massMin< 0. ) _massMin=0.; EvtParticle* par=pmeson->getParent(); double maxMass=-1.; if ( par != 0 ) { if ( par->hasValidP4() ) maxMass=par->mass(); for ( size_t i=0;igetNDaug();i++) { EvtParticle *tDaug=par->getDaug(i); if ( pmeson != tDaug ) maxMass-=EvtPDL::getMinMass(tDaug->getId()); } } EvtId *dauId=0; double *dauMasses=0; size_t nDaug = pmeson->getNDaug(); if ( nDaug > 0) { dauId=new EvtId[nDaug]; dauMasses=new double[nDaug]; for (size_t j=0;jgetDaug(j)->getId(); dauMasses[j]=pmeson->getDaug(j)->mass(); } } EvtId *parId=0; EvtId *othDaugId=0; EvtParticle *tempPar=pmeson->getParent(); if (tempPar) { parId=new EvtId(tempPar->getId()); if ( tempPar->getNDaug()==2 ) { if ( tempPar->getDaug(0) == pmeson ) othDaugId=new EvtId(tempPar->getDaug(1)->getId()); else othDaugId=new EvtId(tempPar->getDaug(0)->getId()); } } if ( nDaug!=2) return calBreitWignerBasic(maxMass); if ( _width< 0.00001) return 1.0; //first figure out L - take the lowest allowed. EvtSpinType::spintype spinD1=EvtPDL::getSpinType(dauId[0]); EvtSpinType::spintype spinD2=EvtPDL::getSpinType(dauId[1]); int t1=EvtSpinType::getSpin2(spinD1); int t2=EvtSpinType::getSpin2(spinD2); int t3=EvtSpinType::getSpin2(_spin); int Lmin=-10; // allow for special cases. if (Lmin<-1 ) { //There are some things I don't know how to deal with if ( t3>4) return calBreitWignerBasic(maxMass); if ( t1>4) return calBreitWignerBasic(maxMass); if ( t2>4) return calBreitWignerBasic(maxMass); //figure the min and max allowwed "spins" for the daughters state Lmin=std::max(t3-t2-t1,std::max(t2-t3-t1,t1-t3-t2)); if (Lmin<0) Lmin=0; assert(Lmin==0||Lmin==2||Lmin==4); } //double massD1=EvtPDL::getMeanMass(dauId[0]); //double massD2=EvtPDL::getMeanMass(dauId[1]); double massD1=dauMasses[0]; double massD2=dauMasses[1]; // I'm not sure how to define the vertex factor here - so retreat to nonRel code. if ( (massD1+massD2)> _mass ) return calBreitWignerBasic(maxMass); //parent vertex factor not yet implemented double massOthD=-10.; double massParent=-10.; int birthl=-10; if ( othDaugId) { EvtSpinType::spintype spinOth=EvtPDL::getSpinType(*othDaugId); EvtSpinType::spintype spinPar=EvtPDL::getSpinType(*parId); int tt1=EvtSpinType::getSpin2(spinOth); int tt2=EvtSpinType::getSpin2(spinPar); int tt3=EvtSpinType::getSpin2(_spin); //figure the min and max allowwed "spins" for the daughters state if ( (tt1<=4) && ( tt2<=4) ) { birthl=std::max(tt3-tt2-tt1,std::max(tt2-tt3-tt1,tt1-tt3-tt2)); if (birthl<0) birthl=0; massOthD=EvtPDL::getMeanMass(*othDaugId); massParent=EvtPDL::getMeanMass(*parId); } } double massM=_massMax; if ( (maxMass > -0.5) && (maxMass < massM) ) massM=maxMass; //special case... if the parent mass is _fixed_ we can do a little better //and only for a two body decay as that seems to be where we have problems // Define relativistic propagator amplitude EvtTwoBodyVertex vd(massD1,massD2,_mass,Lmin/2); vd.set_f(_blatt); EvtPropBreitWignerRel bw(_mass,_width); EvtMassAmp amp(bw,vd); // if ( _fixMassForMax) amp.fixUpMassForMax(); // else std::cout << "problem problem\n"; if ( _includeDecayFact) { amp.addDeathFact(); amp.addDeathFactFF(); } if ( massParent>-1.) { if ( _includeBirthFact ) { EvtTwoBodyVertex vb(_mass,massOthD,massParent,birthl/2); amp.setBirthVtx(vb); amp.addBirthFact(); amp.addBirthFactFF(); } } EvtAmpPdf pdf(amp); double ampVal = sqrt(pdf.evaluate(point)); if ( parId) delete parId; if ( othDaugId) delete othDaugId; if ( dauId) delete [] dauId; if ( dauMasses) delete [] dauMasses; return ampVal; } double EvtPropSLPole::calcMaxProb( EvtId parent, EvtId meson, EvtId lepton, EvtId nudaug, EvtSemiLeptonicFF *FormFactors ) { //This routine takes the arguements parent, meson, and lepton //number, and a form factor model, and returns a maximum //probability for this semileptonic form factor model. A //brute force method is used. The 2D cos theta lepton and //q2 phase space is probed. //Start by declaring a particle at rest. //It only makes sense to have a scalar parent. For now. //This should be generalized later. EvtScalarParticle *scalar_part; EvtParticle *root_part; scalar_part=new EvtScalarParticle; //cludge to avoid generating random numbers! scalar_part->noLifeTime(); EvtVector4R p_init; p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0); scalar_part->init(parent,p_init); root_part=(EvtParticle *)scalar_part; // root_part->set_type(EvtSpinType::SCALAR); root_part->setDiagonalSpinDensity(); EvtParticle *daughter, *lep, *trino; EvtAmp amp; EvtId listdaug[3]; listdaug[0] = meson; listdaug[1] = lepton; listdaug[2] = nudaug; amp.init(parent,3,listdaug); root_part->makeDaughters(3,listdaug); daughter=root_part->getDaug(0); lep=root_part->getDaug(1); trino=root_part->getDaug(2); EvtDecayBase *decayer; decayer = EvtDecayTable::getInstance()->getDecayFunc(daughter); if ( decayer ) { daughter->makeDaughters(decayer->nRealDaughters(),decayer->getDaugs()); for(int ii=0; iinRealDaughters(); ii++){ daughter->getDaug(ii)->setMass(EvtPDL::getMeanMass(daughter->getDaug(ii)->getId())); } } //cludge to avoid generating random numbers! daughter->noLifeTime(); lep->noLifeTime(); trino->noLifeTime(); //Initial particle is unpolarized, well it is a scalar so it is //trivial EvtSpinDensity rho; rho.setDiag(root_part->getSpinStates()); double mass[3]; double m = root_part->mass(); EvtVector4R p4meson, p4lepton, p4nu, p4w; double q2max; double q2, elepton, plepton; int i,j; double erho,prho,costl; double maxfoundprob = 0.0; double prob = -10.0; int massiter; for (massiter=0;massiter<3;massiter++){ mass[0] = EvtPDL::getMeanMass(meson); mass[1] = EvtPDL::getMeanMass(lepton); mass[2] = EvtPDL::getMeanMass(nudaug); if ( massiter==1 ) { mass[0] = EvtPDL::getMinMass(meson); } if ( massiter==2 ) { mass[0] = EvtPDL::getMaxMass(meson); if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001; } q2max = (m-mass[0])*(m-mass[0]); //loop over q2 for (i=0;i<25;i++) { q2 = ((i+0.5)*q2max)/25.0; erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m); prho = sqrt(erho*erho-mass[0]*mass[0]); p4meson.set(erho,0.0,0.0,-1.0*prho); p4w.set(m-erho,0.0,0.0,prho); //This is in the W rest frame elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2)); plepton = sqrt(elepton*elepton-mass[1]*mass[1]); double probctl[3]; for (j=0;j<3;j++) { costl = 0.99*(j - 1.0); //These are in the W rest frame. Need to boost out into //the B frame. p4lepton.set(elepton,0.0, plepton*sqrt(1.0-costl*costl),plepton*costl); p4nu.set(plepton,0.0, -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl); EvtVector4R boost((m-erho),0.0,0.0,1.0*prho); p4lepton=boostTo(p4lepton,boost); p4nu=boostTo(p4nu,boost); //Now initialize the daughters... daughter->init(meson,p4meson); lep->init(lepton,p4lepton); trino->init(nudaug,p4nu); calcamp->CalcAmp(root_part,amp,FormFactors); EvtPoint1D *point = new EvtPoint1D(mass[0]); double meson_BWAmp = calBreitWigner(daughter, *point); int list[2]; list[0]=0; list[1]=0; amp.vertex(0,0,amp.getAmp(list)*meson_BWAmp); list[0]=0; list[1]=1; amp.vertex(0,1,amp.getAmp(list)*meson_BWAmp); list[0]=1; list[1]=0; amp.vertex(1,0,amp.getAmp(list)*meson_BWAmp); list[0]=1; list[1]=1; amp.vertex(1,1,amp.getAmp(list)*meson_BWAmp); list[0]=2; list[1]=0; amp.vertex(2,0,amp.getAmp(list)*meson_BWAmp); list[0]=2; list[1]=1; amp.vertex(2,1,amp.getAmp(list)*meson_BWAmp); //Now find the probability at this q2 and cos theta lepton point //and compare to maxfoundprob. //Do a little magic to get the probability!! prob = rho.normalizedProb(amp.getSpinDensity()); probctl[j]=prob; } //probclt contains prob at ctl=-1,0,1. //prob=a+b*ctl+c*ctl^2 double a=probctl[1]; double b=0.5*(probctl[2]-probctl[0]); double c=0.5*(probctl[2]+probctl[0])-probctl[1]; prob=probctl[0]; if (probctl[1]>prob) prob=probctl[1]; if (probctl[2]>prob) prob=probctl[2]; if (fabs(c)>1e-20){ double ctlx=-0.5*b/c; if (fabs(ctlx)<1.0){ double probtmp=a+b*ctlx+c*ctlx*ctlx; if (probtmp>prob) prob=probtmp; } } //EvtGenReport(EVTGEN_DEBUG,"EvtGen") << "prob,probctl:"< maxfoundprob ) { maxfoundprob = prob; } } if ( EvtPDL::getWidth(meson) <= 0.0 ) { //if the particle is narrow dont bother with changing the mass. massiter = 4; } } root_part->deleteTree(); maxfoundprob *=1.1; return maxfoundprob; }