remove dependency to aliroot libraries, access of ESDEvent object through abstract...

[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtPropSLPole.cxx

//--------------------------------------------------------------------------
//
// 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 <stdlib.h>
#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 <string>

EvtPropSLPole::~EvtPropSLPole() {}

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

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

  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 = new EvtSLPoleFF(getNArg(),getArgs());
  
  if ( mesontype==EvtSpinType::SCALAR ) { 
    calcamp = new EvtSemiLeptonicScalarAmp; 
  }
  if ( mesontype==EvtSpinType::VECTOR ) { 
    calcamp = new EvtSemiLeptonicVectorAmp; 
  }
  if ( mesontype==EvtSpinType::TENSOR ) { 
    calcamp = new EvtSemiLeptonicTensorAmp; 
  }

}


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);
  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;i<par->getNDaug();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;j<nDaug;j++) {
       dauId[j]=pmeson->getDaug(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<EvtPoint1D> 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::getDecayFunc(daughter);
  if ( decayer ) {
    daughter->makeDaughters(decayer->nRealDaughters(),decayer->getDaugs());
    for(int ii=0; ii<decayer->nRealDaughters(); 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;
        }

      }

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

}