EvtSemiLeptonicTensorAmp Class Reference

#include <EvtSemiLeptonicTensorAmp.hh>

Inheritance diagram for EvtSemiLeptonicTensorAmp:

List of all members.

Public Member Functions
void	CalcAmp (EvtParticle *parent, EvtAmp &amp, EvtSemiLeptonicFF *FormFactors)
double	CalcMaxProb (EvtId parent, EvtId meson, EvtId lepton, EvtId nudaug, EvtSemiLeptonicFF *FormFactors)

---

Detailed Description

Definition at line 30 of file EvtSemiLeptonicTensorAmp.hh.

---

Member Function Documentation

void EvtSemiLeptonicTensorAmp::CalcAmp	(	EvtParticle *	parent,
		EvtAmp &	amp,
		EvtSemiLeptonicFF *	FormFactors
	)			[virtual]

Reimplemented from EvtSemiLeptonicAmp.

Definition at line 36 of file EvtSemiLeptonicTensorAmp.cc.

References conj(), EvtTensor4C::cont2(), directProd(), dual(), EvtParticle::epsTensorParent(), calibUtil::ERROR, EvtLeptonVACurrent(), EvtParticle::getDaug(), EvtParticle::getId(), EvtPDL::getId(), EvtParticle::getP4(), EvtSemiLeptonicFF::gettensorff(), EvtVector4R::mass(), EvtParticle::mass(), EvtVector4R::mass2(), q, report(), EvtVector4R::set(), EvtParticle::spParent(), EvtParticle::spParentNeutrino(), and EvtAmp::vertex().

                                                                         {
  static EvtId EM=EvtPDL::getId("e-");
  static EvtId MUM=EvtPDL::getId("mu-");
  static EvtId TAUM=EvtPDL::getId("tau-");
  static EvtId EP=EvtPDL::getId("e+");
  static EvtId MUP=EvtPDL::getId("mu+");
  static EvtId TAUP=EvtPDL::getId("tau+");

  static EvtId D0=EvtPDL::getId("D0");
  static EvtId D0B=EvtPDL::getId("anti-D0");
  static EvtId DP=EvtPDL::getId("D+");
  static EvtId DM=EvtPDL::getId("D-");
  static EvtId DSM=EvtPDL::getId("D_s-");
  static EvtId DSP=EvtPDL::getId("D_s+");

  //Add the lepton and neutrino 4 momenta to find q2

  EvtVector4R q = parent->getDaug(1)->getP4() 
                    + parent->getDaug(2)->getP4(); 
  double q2 = (q.mass2());

  double hf,kf,bpf,bmf;

  FormFactors->gettensorff(parent->getId(),
                           parent->getDaug(0)->getId(),
                           q2,
                           parent->getDaug(0)->mass(),
                           &hf, 
                           &kf, 
                           &bpf, 
                           &bmf);


  double costhl_flag = 1.0;

  if(parent->getId()==D0||parent->getId()==D0B||
     parent->getId()==DP||parent->getId()==DM) {
    costhl_flag = -1.0;
  }
  if(parent->getId()==DSP||parent->getId()==DSM) {
    costhl_flag = -1.0;
  }
  hf = hf * costhl_flag;

  EvtVector4R p4b;
  p4b.set(parent->mass(),0.0,0.0,0.0);
 
  EvtVector4R p4meson = parent->getDaug(0)->getP4();
 
  EvtVector4C l1,l2;

  EvtId l_num = parent->getDaug(1)->getId();

  EvtVector4C ep_meson_b[5];

  ep_meson_b[0] = ((parent->getDaug(0)->epsTensorParent(0)).cont2(p4b)).conj();
  ep_meson_b[1] = ((parent->getDaug(0)->epsTensorParent(1)).cont2(p4b)).conj();
  ep_meson_b[2] = ((parent->getDaug(0)->epsTensorParent(2)).cont2(p4b)).conj();
  ep_meson_b[3] = ((parent->getDaug(0)->epsTensorParent(3)).cont2(p4b)).conj();
  ep_meson_b[4] = ((parent->getDaug(0)->epsTensorParent(4)).cont2(p4b)).conj();

  EvtVector4R pp,pm;

  pp=p4b+p4meson;
  pm=p4b-p4meson;

  //lange - October 31,2002 - try to lessen the mass dependence of probmax
  double q2max = p4b.mass2() + p4meson.mass2() - 2.0*p4b.mass()*p4meson.mass();
  double q2maxin=1.0/q2max;

  EvtComplex ep_meson_bb[5];

  ep_meson_bb[0]=ep_meson_b[0]*(p4b);
  ep_meson_bb[1]=ep_meson_b[1]*(p4b);
  ep_meson_bb[2]=ep_meson_b[2]*(p4b);
  ep_meson_bb[3]=ep_meson_b[3]*(p4b);
  ep_meson_bb[4]=ep_meson_b[4]*(p4b);


  EvtVector4C tds0,tds1,tds2,tds3,tds4;

  EvtTensor4C tds;
  if (l_num==EM||l_num==MUM||l_num==TAUM){
    EvtTensor4C tdual=EvtComplex(0.0,hf)*dual(directProd(pp,pm));
    tds0=tdual.cont2(ep_meson_b[0])
      -kf*ep_meson_b[0]
      -bpf*ep_meson_bb[0]*pp-bmf*ep_meson_bb[0]*pm;
    tds0*=q2maxin;

    tds1=tdual.cont2(ep_meson_b[1])
      -kf*ep_meson_b[1]
      -bpf*ep_meson_bb[1]*pp-bmf*ep_meson_bb[1]*pm;
    tds1*=q2maxin;

    tds2=tdual.cont2(ep_meson_b[2])
      -kf*ep_meson_b[2]
      -bpf*ep_meson_bb[2]*pp-bmf*ep_meson_bb[2]*pm;
    tds2*=q2maxin;

    tds3=tdual.cont2(ep_meson_b[3])
      -kf*ep_meson_b[3]
      -bpf*ep_meson_bb[3]*pp-bmf*ep_meson_bb[3]*pm;
    tds3*=q2maxin;

    tds4=tdual.cont2(ep_meson_b[4])
      -kf*ep_meson_b[4]
      -bpf*ep_meson_bb[4]*pp-bmf*ep_meson_bb[4]*pm;
    tds4*=q2maxin;


    l1=EvtLeptonVACurrent(parent->getDaug(1)->spParent(0),
                          parent->getDaug(2)->spParentNeutrino());
    l2=EvtLeptonVACurrent(parent->getDaug(1)->spParent(1),
                          parent->getDaug(2)->spParentNeutrino());
  }
  else{
    if (l_num==EP||l_num==MUP||l_num==TAUP){
     EvtTensor4C tdual=EvtComplex(0.0,-hf)*dual(directProd(pp,pm));
      tds0=tdual.cont2(ep_meson_b[0])
        -kf*ep_meson_b[0]
        -bpf*ep_meson_bb[0]*pp-bmf*ep_meson_bb[0]*pm;
      tds0*=q2maxin;

      tds1=tdual.cont2(ep_meson_b[1])
        -kf*ep_meson_b[1]
        -bpf*ep_meson_bb[1]*pp-bmf*ep_meson_bb[1]*pm;
      tds1*=q2maxin;

      tds2=tdual.cont2(ep_meson_b[2])
        -kf*ep_meson_b[2]
        -bpf*ep_meson_bb[2]*pp-bmf*ep_meson_bb[2]*pm;
      tds2*=q2maxin;

      tds3=tdual.cont2(ep_meson_b[3])
        -kf*ep_meson_b[3]
        -bpf*ep_meson_bb[3]*pp-bmf*ep_meson_bb[3]*pm;
      tds3*=q2maxin;

      tds4=tdual.cont2(ep_meson_b[4])
        -kf*ep_meson_b[4]
        -bpf*ep_meson_bb[4]*pp-bmf*ep_meson_bb[4]*pm;
      tds4*=q2maxin;

      l1=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
                            parent->getDaug(1)->spParent(0));
      l2=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
                            parent->getDaug(1)->spParent(1));
    }
    else{
      report(ERROR,"EvtGen") << "dfnb89agngri wrong lepton number\n";
    }
  }
 
  amp.vertex(0,0,l1*tds0);
  amp.vertex(0,1,l2*tds0);

  amp.vertex(1,0,l1*tds1);
  amp.vertex(1,1,l2*tds1);

  amp.vertex(2,0,l1*tds2);
  amp.vertex(2,1,l2*tds2);

  amp.vertex(3,0,l1*tds3);
  amp.vertex(3,1,l2*tds3);

  amp.vertex(4,0,l1*tds4);
  amp.vertex(4,1,l2*tds4);

  return;
 
}

double EvtSemiLeptonicAmp::CalcMaxProb	(	EvtId	parent,
		EvtId	meson,
		EvtId	lepton,
		EvtId	nudaug,
		EvtSemiLeptonicFF *	FormFactors
	)			[inherited]

Definition at line 38 of file EvtSemiLeptonicAmp.cc.

References boostTo(), EvtSemiLeptonicAmp::CalcAmp(), EvtParticle::deleteTree(), EvtParticle::getDaug(), EvtPDL::getMass(), EvtPDL::getMaxMass(), EvtPDL::getMeanMass(), EvtPDL::getMinMass(), EvtAmp::getSpinDensity(), EvtParticle::getSpinStates(), EvtPDL::getWidth(), genRecEmupikp::i, EvtParticle::init(), EvtAmp::init(), EvtScalarParticle::init(), ganga-rec::j, EvtParticle::makeDaughters(), EvtParticle::mass(), mass, EvtParticle::noLifeTime(), EvtSpinDensity::NormalizedProb(), EvtVector4R::set(), EvtSpinDensity::SetDiag(), and EvtParticle::setDiagonalSpinDensity().

Referenced by EvtSLPole::initProbMax(), EvtSLBKPole::initProbMax(), EvtHQET2::initProbMax(), and EvtHQET::initProbMax().

                                                      {

  //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);

  //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(root_part,amp,FormFactors);

        //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;
        } 

      }

      //report(DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" "
      //                            << probctl[0]<<" "
      //                            << probctl[1]<<" "
      //                            << probctl[2]<<endl;

      if ( prob > 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;
  
}

---