Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

EvtBtoXsll Class Reference

#include <EvtBtoXsll.hh>

Inheritance diagram for EvtBtoXsll:

EvtDecayIncoherent EvtDecayBase List of all members.

Public Member Functions

void checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1)
void checkNDaug (int d1, int d2=-1)
void checkQ ()
void checkSpinDaughter (int d1, EvtSpinType::spintype sp)
void checkSpinParent (EvtSpinType::spintype sp)
EvtDecayBaseclone ()
virtual void command (std::string cmd)
virtual std::string commandName ()
void decay (EvtParticle *p)
void disableCheckQ ()
 EvtBtoXsll ()
double getArg (int j)
double * getArgs ()
std::string * getArgsStr ()
std::string getArgStr (int j)
double getBranchingFraction ()
EvtId getDaug (int i)
EvtIdgetDaugs ()
int getDSum ()
std::string getModelName ()
void getName (std::string &name)
int getNArg ()
int getNDaug ()
EvtId getParentId ()
int getPHOTOS ()
double getProbMax (double prob)
void init ()
void initProbMax ()
int isDaughterSpinDensitySet (int daughter)
void makeDecay (EvtParticle *p)
virtual bool matchingDecay (const EvtDecayBase &other) const
void noProbMax ()
virtual int nRealDaughters ()
void printSummary ()
double resetProbMax (double prob)
void saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr)
void setDaughterSpinDensity (int daughter)
void setPHOTOS ()
void setProbMax (double prbmx)
void setSummary ()
void setVerbose ()
int summary ()
int verbose ()
virtual ~EvtBtoXsll ()

Static Public Member Functions

void findMass (EvtParticle *p)
void findMasses (EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10])
double findMaxMass (EvtParticle *p)

Protected Member Functions

bool daugsDecayedByParentModel ()

Protected Attributes

bool _daugsDecayedByParentModel

Private Attributes

EvtBtoXsllUtil_calcprob
double _dGdsdupProbMax
double _dGdsProbMax
double _mb
double _mq
double _ms
double _mxmin
double _pf

Constructor & Destructor Documentation

EvtBtoXsll::EvtBtoXsll  )  [inline]
 

00032 {}

EvtBtoXsll::~EvtBtoXsll  )  [virtual]
 

00037 {}


Member Function Documentation

void EvtDecayBase::checkNArg int  a1,
int  a2 = -1,
int  a3 = -1,
int  a4 = -1
[inherited]
 

00483                                                            {
00484 
00485   if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) {
00486     report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "<<endl;
00487     report(ERROR,"EvtGen") << a1<<endl;; 
00488     if ( a2>-1) {
00489       report(ERROR,"EvtGen") << " or " << a2<<endl; 
00490     }
00491     if ( a3>-1) {
00492       report(ERROR,"EvtGen") << " or " << a3<<endl; 
00493     }
00494     if ( a4>-1) {
00495       report(ERROR,"EvtGen") << " or " << a4<<endl; 
00496     }
00497     report(ERROR,"EvtGen") << " arguments but found:"<< _narg << endl;
00498     printSummary();
00499     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00500     ::abort();
00501 
00502   } 
00503 
00504 }

void EvtDecayBase::checkNDaug int  d1,
int  d2 = -1
[inherited]
 

00505                                            {
00506 
00507   if ( _ndaug != d1 && _ndaug != d2 ) {
00508     report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected ";
00509     report(ERROR,"EvtGen") << d1; 
00510     if ( d2>-1) {
00511       report(ERROR,"EvtGen") << " or " << d2; 
00512     }
00513     report(ERROR,"EvtGen") << " daughters but found:"<< _ndaug << endl;
00514     printSummary();
00515     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00516     ::abort();
00517   } 
00518 
00519 }

void EvtDecayBase::checkQ  )  [inherited]
 

00035                           {
00036   int i;
00037   int q=0;
00038   int qpar;
00039 
00040   //If there are no daughters (jetset etc) then we do not
00041   //want to do this test.  Why?  Because sometimes the parent
00042   //will have a nonzero charge.
00043 
00044   if ( _ndaug != 0) {
00045     for(i=0; i<_ndaug; i++ ) {
00046       q += EvtPDL::chg3(_daug[i]);
00047     }
00048     qpar = EvtPDL::chg3(_parent);
00049 
00050     if ( q != qpar ) {
00051       report(ERROR,"EvtGen") <<_modelname.c_str()<< " generator expected "
00052                              << " charge to be conserved, found:"<<endl;
00053       report(ERROR,"EvtGen") << "Parent charge of "<<(qpar/3)<<endl;
00054       report(ERROR,"EvtGen") << "Sum of daughter charge of "<<(q/3)<<endl;
00055       report(ERROR,"EvtGen") << "The parent is "<< EvtPDL::name(_parent).c_str()<<endl;
00056       for(i=0; i<_ndaug; i++ ) {
00057       report(ERROR,"EvtGen") << "Daughter "<< EvtPDL::name(_daug[i]).c_str()<<endl;
00058       }
00059       report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00060       
00061       ::abort();
00062     }
00063   }
00064 }

void EvtDecayBase::checkSpinDaughter int  d1,
EvtSpinType::spintype  sp
[inherited]
 

00534                                                                    {
00535 
00536   EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getDaug(d1));
00537   if ( parenttype != sp ) {
00538     report(ERROR,"EvtGen") << _modelname.c_str() 
00539                            << " did not get the correct daughter spin d=" 
00540                            << d1 << endl;
00541     printSummary();
00542     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00543     ::abort();
00544   } 
00545 
00546 }

void EvtDecayBase::checkSpinParent EvtSpinType::spintype  sp  )  [inherited]
 

00521                                                          {
00522 
00523   EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
00524   if ( parenttype != sp ) {
00525     report(ERROR,"EvtGen") << _modelname.c_str() 
00526                            << " did not get the correct parent spin\n";
00527     printSummary();
00528     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00529     ::abort();
00530   } 
00531 
00532 }

EvtDecayBase * EvtBtoXsll::clone  )  [virtual]
 

Implements EvtDecayBase.

00045                                {
00046 
00047   return new EvtBtoXsll;
00048 
00049 }

void EvtDecayBase::command std::string  cmd  )  [virtual, inherited]
 

Reimplemented in EvtJetSet, EvtLunda, EvtLundCharm, EvtPythia, and EvtTauola.

00132                                        {
00133   report(ERROR,"EvtGen") << "Should never call EvtDecayBase::command"<<endl;
00134   ::abort();
00135 }

std::string EvtDecayBase::commandName  )  [virtual, inherited]
 

Reimplemented in EvtJetSet, EvtLunda, EvtLundCharm, EvtPythia, and EvtTauola.

00129                                    {
00130   return std::string("");
00131 }

bool EvtDecayBase::daugsDecayedByParentModel  )  [inline, protected, inherited]
 

00111 {return _daugsDecayedByParentModel;}

void EvtBtoXsll::decay EvtParticle p  )  [virtual]
 

Implements EvtDecayBase.

00196                                       {
00197 
00198   p->makeDaughters(getNDaug(),getDaugs());
00199 
00200   EvtParticle* xhadron = p->getDaug(0);
00201   EvtParticle* leptonp = p->getDaug(1);
00202   EvtParticle* leptonn = p->getDaug(2);
00203 
00204   double mass[3];
00205  
00206   findMasses( p, getNDaug(), getDaugs(), mass );
00207 
00208   double mB = p->mass();
00209   double ml = mass[1];
00210   double pb;
00211 
00212   int im = 0;
00213   static int nmsg = 0;
00214   double xhadronMass = -999.0;
00215 
00216   EvtVector4R p4xhadron;
00217   EvtVector4R p4leptonp;
00218   EvtVector4R p4leptonn;
00219 
00220   // require the hadronic system has mass greater than that of a Kaon pion pair
00221 
00222   //  while (xhadronMass < 0.6333)
00223   // the above minimum value of K+pi mass appears to be too close
00224   // to threshold as far as JETSET is concerned
00225   // (JETSET gets caught in an infinite loop)
00226   // so we choose a lightly larger value for the threshold
00227   while (xhadronMass < _mxmin)
00228   {
00229     im++;
00230 
00231     // Apply Fermi motion and determine effective b-quark mass
00232 
00233     // Old BaBar MC parameters
00234     //    double pf = 0.25;
00235     //    double ms = 0.2;
00236     //    double mq = 0.3;
00237 
00238     double mb = 0.0;
00239 
00240     double xbox, ybox;
00241 
00242     while (mb <= 0.0)
00243     {
00244       pb = _calcprob->FermiMomentum(_pf);
00245 
00246       // effective b-quark mass
00247       mb = mB*mB + _mq*_mq - 2.0*mB*sqrt(pb*pb + _mq*_mq);
00248     }
00249     mb = sqrt(mb);
00250   
00251     //    cout << "b-quark momentum = " << pb << " mass = " <<  mb << endl;
00252 
00253     // generate a dilepton invariant mass
00254 
00255     double s    = 0.0;
00256     double smin = 4.0 * ml * ml;
00257     double smax = (mb - _ms)*(mb - _ms);
00258 
00259     while (s == 0.0)
00260     {
00261       xbox = EvtRandom::Flat(smin, smax);
00262       ybox = EvtRandom::Flat(_dGdsProbMax);
00263       if (ybox < _calcprob->dGdsProb(mb, _ms, ml, xbox)) { s = xbox;}
00264     }
00265 
00266     //    cout << "dGdsProb(s) = " << _calcprob->dGdsProb(mb, _ms, ml, s)
00267     //         << " for s = " << s << endl;
00268 
00269     // two-body decay of b quark at rest into s quark and dilepton pair:
00270     // b -> s (ll)
00271 
00272     EvtVector4R p4sdilep[2];
00273 
00274     double msdilep[2];
00275     msdilep[0] = _ms;
00276     msdilep[1] = sqrt(s);
00277 
00278     EvtGenKine::PhaseSpace(2, msdilep, p4sdilep, mb);
00279 
00280     // generate dilepton decay with the expected asymmetry: (ll) -> l+ l-
00281 
00282     EvtVector4R p4ll[2];
00283 
00284     double mll[2];
00285     mll[0] = ml;
00286     mll[1] = ml;
00287 
00288     double tmp = 0.0;
00289 
00290     while (tmp == 0.0)
00291     {
00292       // (ll) -> l+ l- decay in dilepton rest frame
00293 
00294       EvtGenKine::PhaseSpace(2, mll, p4ll, msdilep[1]);
00295 
00296       // boost to b-quark rest frame
00297 
00298       p4ll[0] = boostTo(p4ll[0], p4sdilep[1]);
00299       p4ll[1] = boostTo(p4ll[1], p4sdilep[1]);
00300 
00301       // compute kinematical variable u
00302 
00303       EvtVector4R p4slp = p4sdilep[0] + p4ll[0];
00304       EvtVector4R p4sln = p4sdilep[0] + p4ll[1];
00305 
00306       double u = p4slp.mass2() - p4sln.mass2();
00307 
00308       ybox = EvtRandom::Flat(_dGdsdupProbMax);
00309 
00310       double prob = _calcprob->dGdsdupProb(mb, _ms, ml, s, u);
00311       if (prob > _dGdsdupProbMax && nmsg < 20)
00312       {
00313         report(INFO,"EvtGen") << "d2gdsdup GT d2gdsdup_max:" << prob
00314              << " " << _dGdsdupProbMax
00315              << " for s = " << s << " u = " << u << " mb = " << mb << endl;
00316         nmsg++;
00317       }
00318       if (ybox < prob)
00319       {
00320         tmp = 1.0;
00321         //        cout << "dGdsdupProb(s) = " << prob
00322         //             << " for u = " << u << endl;
00323       }
00324     }
00325 
00326     // assign 4-momenta to valence quarks inside B meson in B rest frame
00327 
00328     double phi   = EvtRandom::Flat( EvtConst::twoPi );
00329     double costh = EvtRandom::Flat( -1.0, 1.0 );
00330     double sinth = sqrt(1.0 - costh*costh);
00331 
00332     // b-quark four-momentum in B meson rest frame
00333 
00334     EvtVector4R p4b(sqrt(mb*mb + pb*pb),
00335                     pb*sinth*sin(phi),
00336                     pb*sinth*cos(phi),
00337                     pb*costh);
00338 
00339     // B meson in its rest frame
00340     //
00341     //    EvtVector4R p4B(mB, 0.0, 0.0, 0.0);
00342     //
00343     // boost B meson to b-quark rest frame
00344     //
00345     //    p4B = boostTo(p4B, p4b);
00346     //
00347     //    cout << " B meson mass in b-quark rest frame = " << p4B.mass() << endl;
00348 
00349     // boost s, l+ and l- to B meson rest frame
00350 
00351     //    EvtVector4R p4s = boostTo(p4sdilep[0], p4B);
00352     //    p4leptonp       = boostTo(p4ll[0],     p4B);
00353     //    p4leptonn       = boostTo(p4ll[1],     p4B);
00354 
00355     EvtVector4R p4s = boostTo(p4sdilep[0], p4b);
00356     p4leptonp       = boostTo(p4ll[0],     p4b);
00357     p4leptonn       = boostTo(p4ll[1],     p4b);
00358 
00359     // spectator quark in B meson rest frame
00360 
00361     EvtVector4R p4q( sqrt(pb*pb + _mq*_mq), -p4b.get(1), -p4b.get(2), -p4b.get(3) );
00362 
00363     // hadron system in B meson rest frame
00364 
00365     p4xhadron = p4s + p4q;
00366     xhadronMass = p4xhadron.mass();
00367 
00368     //    cout << "Xs mass = " << xhadronMass << " trial " << im << endl;
00369   }
00370 
00371   // initialize the decay products
00372 
00373   xhadron->init(getDaug(0), p4xhadron);
00374 
00375   // For B-bar mesons (i.e. containing a b quark) we have the normal
00376   // order of leptons
00377   if ( p->getId() == EvtPDL::getId("anti-B0") ||
00378        p->getId() == EvtPDL::getId("B-") )
00379   {
00380     leptonp->init(getDaug(1), p4leptonp);
00381     leptonn->init(getDaug(2), p4leptonn);
00382   }
00383   // For B mesons (i.e. containing a b-bar quark) we need to flip the
00384   // role of the positive and negative leptons in order to produce the
00385   // correct forward-backward asymmetry between the two leptons
00386   else
00387   {
00388     leptonp->init(getDaug(1), p4leptonn);
00389     leptonn->init(getDaug(2), p4leptonp);
00390   }
00391 
00392   return ;
00393 }

void EvtDecayBase::disableCheckQ  )  [inline, inherited]
 

00062 {_chkCharge=0;};

void EvtDecayBase::findMass EvtParticle p  )  [static, inherited]
 

00347                                           {
00348 
00349   //Need to also check that this mass does not screw
00350   //up the parent
00351   //This code assumes that for the ith daughter, 0..i-1
00352   //already have a mass
00353   double maxOkMass=findMaxMass(p);
00354 
00355   int count=0;
00356   double mass;
00357   bool massOk=false;
00358   int i;
00359   while (!massOk) { 
00360     count++;
00361     if ( count > 10000 ) {
00362       report(INFO,"EvtGen") << "Can not find a valid mass for: " << EvtPDL::name(p->getId()).c_str() <<endl;
00363       report(INFO,"EvtGen") << "Now printing parent and/or grandparent tree\n";
00364       if ( p->getParent() ) {
00365         if ( p->getParent()->getParent() ) {
00366           p->getParent()->getParent()->printTree();
00367           report(INFO,"EvtGen") << p->getParent()->getParent()->mass() <<endl;
00368           report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
00369         }
00370         else{
00371           p->getParent()->printTree();
00372           report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
00373         }
00374       }
00375       else  p->printTree();
00376       report(INFO,"EvtGen") << "maxokmass=" << maxOkMass << " " << EvtPDL::getMinMass(p->getId()) << " " << EvtPDL::getMaxMass(p->getId())<<endl;
00377       if ( p->getNDaug() ) { 
00378         for (i=0; i<p->getNDaug(); i++) {
00379           report(INFO,"EvtGen") << p->getDaug(i)->mass()<<" ";
00380         }
00381         report(INFO,"EvtGen") << endl;
00382       }
00383       if ( maxOkMass >= EvtPDL::getMinMass(p->getId()) ) {
00384         report(INFO,"EvtGen") << "taking a default value\n";
00385         p->setMass(maxOkMass);
00386         return;
00387       } 
00388       assert(0);
00389     }
00390     mass = EvtPDL::getMass(p->getId());
00391     //Just need to check that this mass is > than
00392     //the mass of all daughters
00393     double massSum=0.;
00394     if ( p->getNDaug() ) { 
00395       for (i=0; i<p->getNDaug(); i++) {
00396         massSum+= p->getDaug(i)->mass();
00397       }
00398     }
00399     //some special cases are handled with 0 (stable) or 1 (k0->ks/kl) daughters
00400     if (p->getNDaug()<2)  massOk=true;
00401     if ( p->getParent() ) {
00402       if ( p->getParent()->getNDaug()==1 ) massOk=true;
00403     }
00404     if ( !massOk ) { 
00405       if (massSum < mass) massOk=true;
00406       if ( mass> maxOkMass) massOk=false;
00407     }
00408   }
00409 
00410   p->setMass(mass);
00411   
00412 }

void EvtDecayBase::findMasses EvtParticle p,
int  ndaugs,
EvtId  daugs[10],
double  masses[10]
[static, inherited]
 

00416                                                                      {
00417 
00418   int i;
00419   double mass_sum;
00420 
00421   int count=0;
00422 
00423   if (!( p->firstornot() )) {
00424     for (i = 0; i < ndaugs; i++ ) {
00425       masses[i] = p->getDaug(i)->mass();
00426     } //for
00427   } //if
00428   else {
00429     p->setFirstOrNot();
00430     // if only one daughter do it
00431 
00432     if (ndaugs==1) {
00433       masses[0]=p->mass();
00434       return;
00435     }
00436     
00437     //until we get a combo whose masses are less than _parent mass.
00438     do {
00439       mass_sum = 0.0;
00440 
00441       for (i = 0; i < ndaugs; i++ ) {
00442         masses[i] = EvtPDL::getMass(daugs[i]);
00443         mass_sum = mass_sum + masses[i];
00444       } 
00445 
00446       count++;
00447 
00448      
00449       if(count==10000) {
00450         report(ERROR,"EvtGen") <<"Decaying particle:"<<
00451           EvtPDL::name(p->getId()).c_str()<<" (m="<<p->mass()<<")"<<endl;
00452         report(ERROR,"EvtGen") <<"To the following daugthers"<<endl;
00453         for (i = 0; i < ndaugs; i++ ) {
00454           report(ERROR,"EvtGen") <<  
00455             EvtPDL::name(daugs[i]).c_str() << endl;
00456         } 
00457         report(ERROR,"EvtGen") << "Has been rejected "<<count
00458                                << " times, will now take minimal masses "
00459                                << " of daugthers"<<endl;
00460         
00461         mass_sum=0.;
00462         for (i = 0; i < ndaugs; i++ ) {
00463           masses[i] = EvtPDL::getMinMass(daugs[i]);
00464           mass_sum = mass_sum + masses[i];
00465         } 
00466         if (mass_sum > p->mass()){
00467           report(ERROR,"EvtGen") << "Parent mass="<<p->mass()
00468                                  << "to light for daugthers."<<endl
00469                                  << "Will throw the event away."<<endl;
00470           //dont terminate - start over on the event.
00471           EvtStatus::setRejectFlag();
00472           mass_sum=0.;
00473           //      ::abort();
00474         }
00475 
00476       }
00477     } while ( mass_sum > p->mass());
00478   } //else
00479   
00480   return;
00481 }       

double EvtDecayBase::findMaxMass EvtParticle p  )  [static, inherited]
 

00312                                                {
00313 
00314   
00315   double maxOkMass=EvtPDL::getMaxMass(p->getId());
00316 
00317   //protect against vphotons
00318   if ( maxOkMass < 0.0000000001 ) return 10000000.;
00319   //and against already determined masses
00320   if ( p->hasValidP4() ) maxOkMass=p->mass();
00321 
00322   EvtParticle *par=p->getParent();
00323   if ( par ) {
00324     double maxParMass=findMaxMass(par);
00325     int i;
00326     double minDaugMass=0.;
00327     for(i=0;i<par->getNDaug();i++){
00328       EvtParticle *dau=par->getDaug(i);
00329       if ( dau!=p) {
00330         // it might already have a mass
00331         if ( dau->isInitialized() || dau->hasValidP4() )
00332           minDaugMass+=dau->mass();
00333         else
00334         //give it a bit of phase space 
00335           minDaugMass+=1.000001*EvtPDL::getMinMass(dau->getId());
00336       }
00337     }
00338     if ( maxOkMass>(maxParMass-minDaugMass)) maxOkMass=maxParMass-minDaugMass;
00339   }
00340   return maxOkMass;
00341 }

double EvtDecayBase::getArg int  j  )  [inherited]
 

00565                                  {
00566 
00567   // Verify string
00568 
00569   const char* str = _args[j].c_str();
00570   int i = 0;
00571   while(str[i]!=0){
00572     if (isalpha(str[i]) && str[i]!='e') {
00573 
00574       report(INFO,"EvtGen") << "String " << str << " is not a number" << endl;
00575       assert(0);
00576     }
00577     i++;
00578   }
00579   
00580   char** tc=0; 
00581   return strtod(_args[j].c_str(),tc);
00582 }

double * EvtDecayBase::getArgs  )  [inherited]
 

00548                               {
00549 
00550   if ( _argsD ) return _argsD;
00551   //The user has asked for a list of doubles - the arguments 
00552   //better all be doubles...
00553   if ( _narg==0 ) return _argsD;
00554 
00555   _argsD = new double[_narg];
00556 
00557   int i;
00558   char * tc;
00559   for(i=0;i<_narg;i++) { 
00560     _argsD[i] =  strtod(_args[i].c_str(),&tc);
00561   }
00562   return _argsD;
00563 }

std::string* EvtDecayBase::getArgsStr  )  [inline, inherited]
 

00073 {return _args;}

std::string EvtDecayBase::getArgStr int  j  )  [inline, inherited]
 

00075 {return _args[j];}

double EvtDecayBase::getBranchingFraction  )  [inline, inherited]
 

00061 {return _brfr;}

EvtId EvtDecayBase::getDaug int  i  )  [inline, inherited]
 

00066 {return _daug[i];}

EvtId* EvtDecayBase::getDaugs  )  [inline, inherited]
 

00065 {return _daug;}

int EvtDecayBase::getDSum  )  [inline, inherited]
 

00077 {return _dsum; }

std::string EvtDecayBase::getModelName  )  [inline, inherited]
 

00076 {return _modelname; }

void EvtBtoXsll::getName std::string &  name  )  [virtual]
 

Implements EvtDecayBase.

00039                                              {
00040 
00041   model_name="BTOXSLL";     
00042 
00043 }

int EvtDecayBase::getNArg  )  [inline, inherited]
 

00067 {return _narg;}

int EvtDecayBase::getNDaug  )  [inline, inherited]
 

00064 {return _ndaug;}

EvtId EvtDecayBase::getParentId  )  [inline, inherited]
 

00060 {return _parent;}

int EvtDecayBase::getPHOTOS  )  [inline, inherited]
 

00068 {return _photos;}

double EvtDecayBase::getProbMax double  prob  )  [inherited]
 

00067                                              {
00068 
00069   int i;
00070 
00071   //diagnostics
00072   sum_prob+=prob;
00073   if (prob>max_prob) max_prob=prob;
00074 
00075 
00076   if ( defaultprobmax && ntimes_prob<=500 ) { 
00077     //We are building up probmax with this iteration
00078      ntimes_prob += 1;
00079      if ( prob > probmax ) { probmax = prob;}
00080      if (ntimes_prob==500) { 
00081        probmax*=1.2;
00082      }
00083      return 1000000.0*prob;
00084   }
00085 
00086   if ( prob> probmax*1.0001) {
00087 
00088     report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")";
00089     report(INFO,"") << "("<<_modelname.c_str()<<") ";
00090     report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
00091     for(i=0;i<_ndaug;i++){
00092        report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
00093     }
00094     report(INFO,"") << endl;
00095 
00096     if (defaultprobmax) probmax = prob;
00097 
00098   }
00099 
00100   ntimes_prob += 1;
00101 
00102 
00103   return probmax;
00104 
00105 } //getProbMax

void EvtBtoXsll::init  )  [virtual]
 

Reimplemented from EvtDecayBase.

00052                      {
00053 
00054   // check that there are no arguments
00055 
00056   checkNArg(0,4,5);
00057 
00058   checkNDaug(3);
00059 
00060   // Check that the two leptons are the same type
00061 
00062   EvtId lepton1type = getDaug(1);
00063   EvtId lepton2type = getDaug(2);
00064 
00065   int etyp   = 0;
00066   int mutyp  = 0;
00067   int tautyp = 0;
00068   if ( lepton1type == EvtPDL::getId("e+") ||
00069        lepton1type == EvtPDL::getId("e-") ) { etyp++;}
00070   if ( lepton2type == EvtPDL::getId("e+") ||
00071        lepton2type == EvtPDL::getId("e-") ) { etyp++;}
00072   if ( lepton1type == EvtPDL::getId("mu+") ||
00073        lepton1type == EvtPDL::getId("mu-") ) { mutyp++;}
00074   if ( lepton2type == EvtPDL::getId("mu+") ||
00075        lepton2type == EvtPDL::getId("mu-") ) { mutyp++;}
00076   if ( lepton1type == EvtPDL::getId("tau+") ||
00077        lepton1type == EvtPDL::getId("tau-") ) { tautyp++;}
00078   if ( lepton2type == EvtPDL::getId("tau+") ||
00079        lepton2type == EvtPDL::getId("tau-") ) { tautyp++;}
00080 
00081   if ( etyp != 2 && mutyp != 2 && tautyp != 2 ) {
00082 
00083      report(ERROR,"EvtGen") << "Expect two leptons of the same type in EvtBtoXsll.cc\n";
00084      ::abort();
00085   }
00086 
00087   // Check that the second and third entries are leptons with positive
00088   // and negative charge, respectively
00089 
00090   int lpos = 0;
00091   int lneg = 0;
00092   if ( lepton1type == EvtPDL::getId("e+")  ||
00093        lepton1type == EvtPDL::getId("mu+") ||
00094        lepton1type == EvtPDL::getId("tau+") ) { lpos++;}
00095   if ( lepton2type == EvtPDL::getId("e-")  ||
00096        lepton2type == EvtPDL::getId("mu-") ||
00097        lepton2type == EvtPDL::getId("tau-") ) { lneg++;}
00098 
00099   if ( lpos != 1 || lneg != 1 ) {
00100 
00101      report(ERROR,"EvtGen") << "Expect 2nd and 3rd particles to be positive and negative leptons in EvtBtoXsll.cc\n";
00102      ::abort();
00103   }
00104 
00105  
00106   _mb=4.8;
00107   _ms=0.2;
00108   _mq=0.;
00109   _pf=0.41;
00110   _mxmin=1.1;
00111   if ( getNArg()==4) 
00112   {
00113     // b-quark mass
00114     _mb = getArg(0);
00115     // s-quark mass
00116     _ms = getArg(1);
00117     // spectator quark mass
00118     _mq = getArg(2);
00119     // Fermi motion parameter
00120     _pf = getArg(3);
00121   }
00122   if ( getNArg()==5) 
00123   {
00124     _mxmin = getArg(4);
00125   }
00126 
00127   _calcprob = new EvtBtoXsllUtil;
00128 
00129   double ml = EvtPDL::getMeanMass(getDaug(1));
00130 
00131   // determine the maximum probability density from dGdsProb
00132 
00133   int i, j;
00134   int nsteps  = 100;
00135   double s    = 0.0;
00136   double smin = 4.0 * ml * ml;
00137   double smax = (_mb - _ms)*(_mb - _ms);
00138   double probMax  = -10000.0;
00139   double sProbMax = -10.0;
00140   double uProbMax = -10.0;
00141 
00142   for (i=0;i<nsteps;i++)
00143   {
00144     s = smin + (i+0.0005)*(smax - smin)/(double)nsteps;
00145     double prob = _calcprob->dGdsProb(_mb, _ms, ml, s);
00146     if (prob > probMax)
00147     {
00148       sProbMax = s;
00149       probMax  = prob;
00150     }
00151   }
00152 
00153   _dGdsProbMax = probMax;
00154 
00155   if ( verbose() ) {
00156     report(INFO,"EvtGen") << "dGdsProbMax = " << probMax << " for s = "  << sProbMax << endl;
00157   }
00158 
00159   // determine the maximum probability density from dGdsdupProb
00160 
00161   probMax  = -10000.0;
00162   sProbMax = -10.0;
00163 
00164   for (i=0;i<nsteps;i++)
00165   {
00166     s = smin + (i+0.0005)*(smax - smin)/(double)nsteps;
00167     double umax = sqrt((s - (_mb + _ms)*(_mb + _ms)) * (s - (_mb - _ms)*(_mb - _ms)));
00168     for (j=0;j<nsteps;j++)
00169     {
00170       double u = -umax + (j+0.0005)*(2.0*umax)/(double)nsteps;
00171       double prob = _calcprob->dGdsdupProb(_mb, _ms, ml, s, u);
00172       if (prob > probMax)
00173       {
00174         sProbMax = s;
00175         uProbMax = u;
00176         probMax  = prob;
00177       }
00178     }
00179   }
00180 
00181   _dGdsdupProbMax = probMax;
00182 
00183   if ( verbose() ) {
00184    report(INFO,"EvtGen") << "dGdsdupProbMax = " << probMax << " for s = "  << sProbMax
00185                           << " and u = " << uProbMax << endl;
00186   }
00187 
00188 }

void EvtBtoXsll::initProbMax  )  [virtual]
 

Reimplemented from EvtDecayBase.

00190                             {
00191 
00192   noProbMax();
00193 
00194 }

int EvtDecayIncoherent::isDaughterSpinDensitySet int  daughter  )  [inline, inherited]
 

00041   {return spinDensitySet[daughter];}

void EvtDecayIncoherent::makeDecay EvtParticle p  )  [virtual, inherited]
 

Implements EvtDecayBase.

00030                                                 {
00031 
00032   int i;
00033   //initialize this the hard way..
00034   //Lange June 26, 2000
00035   for (i=0; i<MAX_DAUG; i++ ) { spinDensitySet[i]=0;}
00036   _daugsDecayedByParentModel=false;
00037 
00038   decay(p);
00039   p->setDecayProb(1.);
00040 
00041   EvtSpinDensity rho;
00042 
00043   rho.SetDiag(p->getSpinStates());
00044 
00045   p->setSpinDensityBackward(rho);
00046 
00047   if (getPHOTOS() || EvtRadCorr::alwaysRadCorr()) {
00048     EvtRadCorr::doRadCorr(p);
00049   }
00050 
00051   //Now decay the daughters.
00052 
00053   if ( !daugsDecayedByParentModel()) {
00054     
00055     for(i=0;i<p->getNDaug();i++){
00056       //Need to set the spin density of the daughters to be
00057       //diagonal.
00058       rho.SetDiag(p->getDaug(i)->getSpinStates());
00059       //if (p->getDaug(i)->getNDaug()==0){
00060       //only do this if the user has not already set the 
00061       //spin density matrix herself.
00062       //Lange June 26, 2000
00063       if ( isDaughterSpinDensitySet(i)==0 ) { 
00064         p->getDaug(i)->setSpinDensityForward(rho);
00065       }
00066       else{
00067         //report(INFO,"EvtGen") << "spinDensitymatrix already set!!!\n";
00068         EvtSpinDensity temp=p->getDaug(i)->getSpinDensityForward();
00069         //      report(INFO,"EvtGen") <<temp<<endl;
00070       }
00071       //Now decay the daughter.  Really!
00072       p->getDaug(i)->decay();
00073       //}
00074     } 
00075   }
00076                             
00077 }

bool EvtDecayBase::matchingDecay const EvtDecayBase other  )  const [virtual, inherited]
 

00588                                                                 {
00589 
00590   if ( _ndaug != other._ndaug) return false;
00591   if ( _parent != other._parent) return false;
00592   
00593   std::vector<int> useDs;
00594   for ( unsigned int i=0; i<_ndaug; i++) useDs.push_back(0);
00595 
00596   for ( unsigned int i=0; i<_ndaug; i++) {
00597     bool foundIt=false;
00598     for ( unsigned int j=0; j<_ndaug; j++) {
00599       if ( useDs[j] == 1 ) continue;
00600       if ( _daug[i] == other._daug[j] && _daug[i].getAlias() == other._daug[j].getAlias()) {
00601         foundIt=true;
00602         useDs[j]=1;
00603         break;
00604       }
00605     }
00606     if ( foundIt==false) return false;
00607   }
00608   for ( unsigned int i=0; i<_ndaug; i++) if ( useDs[i]==0) return false;
00609 
00610   return true;
00611 
00612 }

void EvtDecayBase::noProbMax  )  [inherited]
 

00305                             {
00306 
00307   defaultprobmax=0;
00308 
00309 }

virtual int EvtDecayBase::nRealDaughters  )  [inline, virtual, inherited]
 

Reimplemented in EvtBtoKD3P, and EvtVSSBMixCPT.

00105 { return _ndaug;}

void EvtDecayBase::printSummary  )  [inherited]
 

00259                                 {
00260 
00261   int i;
00262 
00263   if (ntimes_prob>0) {
00264 
00265     report(INFO,"EvtGen") << "Calls="<<ntimes_prob<<" eff:"<<
00266       sum_prob/(probmax*ntimes_prob)<<" frac. max:"<<max_prob/probmax;
00267     report(INFO,"") <<" probmax:"<<probmax<<" max:"<<max_prob<<" : ";
00268   }
00269 
00270   report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
00271   for(i=0;i<_ndaug;i++){
00272     report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
00273   }
00274   report(INFO,"") << " ("<<_modelname.c_str()<<"):"<< endl;
00275   
00276   
00277   
00278 }

double EvtDecayBase::resetProbMax double  prob  )  [inherited]
 

00108                                              {
00109   
00110   report(INFO,"EvtGen") << "Reseting prob max\n"; 
00111   report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")";
00112   report(INFO,"") << "("<<_modelname.c_str()<<")";
00113   report(INFO,"") << EvtPDL::getStdHep(_parent)<<"->";
00114   
00115   for( int i=0;i<_ndaug;i++){
00116     report(INFO,"") << EvtPDL::getStdHep(_daug[i]) << " ";
00117   }
00118   report(INFO,"") << endl;
00119   
00120   probmax = 0.0;
00121   defaultprobmax = 0;
00122   ntimes_prob = 0;
00123   
00124   return prob;
00125 
00126 }

void EvtDecayBase::saveDecayInfo EvtId  ipar,
int  ndaug,
EvtId daug,
int  narg,
std::vector< std::string > &  args,
std::string  name,
double  brfr
[inherited]
 

00170                                               {
00171 
00172   int i;
00173 
00174   _brfr=brfr;
00175   _ndaug=ndaug;
00176   _narg=narg;
00177   _parent=ipar; 
00178 
00179   _dsum=0;
00180 
00181   if (_ndaug>0) {
00182     _daug=new EvtId [_ndaug];
00183     for(i=0;i<_ndaug;i++){
00184       _daug[i]=daug[i];
00185       _dsum+=daug[i].getAlias();
00186     }
00187   }
00188   else{
00189     _daug=0;
00190   }
00191 
00192   if (_narg>0) {
00193     _args=new std::string[_narg+1];
00194     for(i=0;i<_narg;i++){
00195       _args[i]=args[i];
00196     }
00197   }
00198   else{
00199      _args = 0;
00200   }
00201 
00202   _modelname=name;
00203 
00204   this->init();
00205   this->initProbMax();
00206 
00207   if (_chkCharge){
00208     this->checkQ();
00209   }
00210 
00211 
00212   if (defaultprobmax){
00213     report(INFO,"EvtGen") << "No default probmax for ";
00214     report(INFO,"") << "("<<_modelname.c_str()<<") ";
00215     report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
00216     for(i=0;i<_ndaug;i++){
00217       report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
00218     }
00219     report(INFO,"") << endl;
00220     report(INFO,"") << "This is fine for development, but must be provided for production."<<endl;
00221     report(INFO,"EvtGen") << "Never fear though - the decay will use the \n";
00222     report(INFO,"EvtGen") << "500 iterations to build up a good probmax \n";
00223     report(INFO,"EvtGen") << "before accepting a decay. "<<endl;
00224   }
00225 
00226 }

void EvtDecayIncoherent::setDaughterSpinDensity int  daughter  )  [inline, inherited]
 

00038   { spinDensitySet[daughter]=1; return;}

void EvtDecayBase::setPHOTOS  )  [inline, inherited]
 

00069 {_photos=1;}

void EvtDecayBase::setProbMax double  prbmx  )  [inherited]
 

00298                                          {
00299 
00300   defaultprobmax=0;
00301   probmax=prbmx;
00302 
00303 }

void EvtDecayBase::setSummary  )  [inline, inherited]
 

00071 {_summary=1;}

void EvtDecayBase::setVerbose  )  [inline, inherited]
 

00070 {_verbose=1;}

int EvtDecayBase::summary  )  [inline, inherited]
 

00078 {return _summary; }

int EvtDecayBase::verbose  )  [inline, inherited]
 

00079 {return _verbose; }


Member Data Documentation

EvtBtoXsllUtil* EvtBtoXsll::_calcprob [private]
 

bool EvtDecayBase::_daugsDecayedByParentModel [protected, inherited]
 

double EvtBtoXsll::_dGdsdupProbMax [private]
 

double EvtBtoXsll::_dGdsProbMax [private]
 

double EvtBtoXsll::_mb [private]
 

double EvtBtoXsll::_mq [private]
 

double EvtBtoXsll::_ms [private]
 

double EvtBtoXsll::_mxmin [private]
 

double EvtBtoXsll::_pf [private]
 


The documentation for this class was generated from the following files:
Generated on Wed Feb 2 16:04:03 2011 for BOSS6.5.5 by  doxygen 1.3.9.1