EvtCalHelAmp.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// 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: EvtCalHelAmp.cc
//
// Description: Routine to decay a particle according th phase space
//
// Modification history:
//
//    RYD       January 8, 1997       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/EvtCalHelAmp.hh"
#include "EvtGenModels/EvtGlobalSet.hh"
#include <string>
 
double EvtCalHelAmp::_H2err=0;
double EvtCalHelAmp::_H1err=0;
double EvtCalHelAmp::_H0err=0;
double EvtCalHelAmp::_H2=0;
double EvtCalHelAmp::_H1=0;
double EvtCalHelAmp::_H0=0;
int EvtCalHelAmp::nevt=0;
EvtCalHelAmp::~EvtCalHelAmp() {}
 
void EvtCalHelAmp::getName(std::string& model_name){
 
  model_name="CALHELAMP";     
 
}
 
EvtDecayBase* EvtCalHelAmp::clone(){
 
  return new EvtCalHelAmp;
 
}
 
 
void EvtCalHelAmp::init(){
  _H2=0;_H1=0;_H0=0;_H2err=0;_H1err=0;_H0err=0;nevt=0;
  // check that there are 1 arguments (the number of parameters)
  checkNArg(4);
  nd=getNArg();
  VC.resize(nd); 
  for(int i=0;i<nd;i++) VC[i].resize(nd);
  ifstream coverrf("myerr.dat");
  if(coverrf==0) {std::cout<<"File of covariant error (myerr.dat)" <<" is not found"<<endl;abort();}
  for(int i=0;i<nd;i++){
    for(int j=0;j<nd;j++){
      double xr;
      coverrf>>xr;
      VC[i][j]=xr;
    }
  }
 
}
 
void EvtCalHelAmp::initProbMax(){
 
  noProbMax();
 
}
 
void EvtCalHelAmp::decay( EvtParticle *p ){
  if(getNArg()!=4) {cout<<"The model have 6 parameters: |g00| phi00 |g22| phi22 |g42| phi42"<<endl; abort();}
 
  double weight = p->initializePhaseSpace(getNDaug(),getDaugs());
  //  std::cout<<"weight= "<<weight<<std::endl;
  //std::cout<<p->getP4().mass()<<" "<<p->getDaug(0)->getP4().mass()<<" "<<p->getDaug(1)->getP4().mass()<<std::endl;
  std::string p0str=EvtPDL::name(p->getId());
  std::string p1str=EvtPDL::name(p->getDaug(0)->getId());
  std::string p2str=EvtPDL::name(p->getDaug(1)->getId());
  if(p1str=="K_2*+" || p2str=="K_2*+" ||p1str=="K_2*0" || p2str=="K_2*0" ) {flag=2;} else
  if(p1str=="K_3*+" || p2str=="K_3*+" ||p1str=="K_3*0" || p2str=="K_3*0" ) {flag=3;} else
  if(p1str=="Zc+"   || p2str=="pi-"   ||p1str=="Zc-"   || p2str=="pi^+"  ) {flag=4;} else
  if(p0str=="Zc+"   && (p1str=="J/psi" || p2str=="J/psi")  ) {flag=5;} 
  //Helicity amplitude for 2+ -> VP
  EvtVector4R p1= p->getDaug(0)->getP4();
  double r= 2*p1.d3mag();
  double R=3.0;
  double pi=3.1415926;
  EvtBlattWeisskopf B0(0,R,r/2);  double b0=B0.get(r/2);
  EvtBlattWeisskopf B1(1,R,r/2);  double b1=B1.get(r/2);
  EvtBlattWeisskopf B2(2,R,r/2);  double b2=B2.get(r/2);
  EvtBlattWeisskopf B3(3,R,r/2);  double b3=B3.get(r/2);
  EvtBlattWeisskopf B4(4,R,r/2);  double b4=B4.get(r/2);
  double g02 =   getArg(0);
  double phi02 = getArg(1)*2*pi;
  double g22 =   getArg(2);
  double phi22 = getArg(3)*2*pi;
  //double g42 =   getArg(4);
  //double phi42 = getArg(5)*2*pi;
  EvtComplex G02(g02*cos(phi02),g02*sin(phi02));
  EvtComplex G22(g22*cos(phi22),g22*sin(phi22));
  //EvtComplex G42(g42*cos(phi42),g42*sin(phi42));
  std::vector<EvtComplex> VG; VG.clear();
  VG.push_back(G02); VG.push_back(G22);// VG.push_back(G42);
  std::vector<double> VH2,VH1,VH0;
  VH2.resize(nd/2);VH1.resize(nd/2);VH0.resize(nd/2);
  if(flag==2){
    VH2[0]=sqrt(2./5.)*b1*r;  VH2[1]= 1/sqrt(10.)*r*r*r*b3; 
    VH1[0]=-1/sqrt(10.)*b1*r; VH1[1]= r*r*r*b3*sqrt(2./5.);
  }else if(flag==3){
    VH2[0]=sqrt(5./14.)*b2*r*r;  VH2[1]= 1/sqrt(7.)*r*r*r*r*b4; 
    VH1[0]=sqrt(1./7.)*b2*r*r;   VH1[1]= -sqrt(5./14.)*r*r*r*r*b4;
  }else if(flag==4 || flag==5){
    VH2[0]=sqrt(1./3.)*b0;  VH2[1]= 1/sqrt(6.)*r*r*b2; 
    VH1[0]=sqrt(1./3.)*b0;  VH1[1]=-2/sqrt(6.)*r*r*b2;
  }else{std::cout<<"Not allowed mode!"<<std::endl; abort();}
  EvtComplex H2 = G02*VH2[0] + G22*VH2[1];
  EvtComplex H1 = G02*VH1[0] + G22*VH1[1];
  
  _H2 += abs2(H2); _H1 += abs2(H1); 
  std::vector<double> DH1,DH2;
  DH2=firstder(VH2);
  DH1=firstder(VH1);
 
  // std::cout<<DH2[0]<<" "<<DH2[1]<<" "<<DH2[2]<<" "<<DH2[3]<<std::endl;
  //std::cout<<DH1[0]<<" "<<DH1[1]<<" "<<DH1[2]<<" "<<DH1[3]<<std::endl;
 
  for(int i=0;i<nd;i++){
    for(int j=0;j<nd;j++){
      _H2err += DH2[i]*DH2[j]*VC[i][j]; 
      _H1err += DH1[i]*DH1[j]*VC[i][j];
    }
  }
  
  nevt++;
 
  return ;
}
 
 
std::vector<double> EvtCalHelAmp::firstder(std::vector<double> vh){
  // vh: vector of helicity part
  double pi=3.1415926;
  std::vector<double> fd;
  double g02 =   getArg(0);
  double phi02 = getArg(1)*2*pi;
  double g22 =   getArg(2);
  double phi22 = getArg(3)*2*pi;
  std::vector<double> vpar;
  vpar.push_back(g02);vpar.push_back(phi02);vpar.push_back(g22);vpar.push_back(phi22);
  EvtComplex G02(vpar[0]*cos(vpar[1]),vpar[0]*sin(vpar[1]));
  EvtComplex G22(vpar[2]*cos(vpar[3]),vpar[2]*sin(vpar[3]));
  EvtComplex H20 = G02*vh[0] + G22*vh[1];
  double hamps0 = abs2(H20);
  for(int i=0;i<vpar.size();i++){
    vpar.clear();
    vpar.push_back(g02);vpar.push_back(phi02);vpar.push_back(g22);vpar.push_back(phi22);
    double dev=0.01;
    vpar[i] += dev;
    EvtComplex G02(vpar[0]*cos(vpar[1]),vpar[0]*sin(vpar[1]));
    EvtComplex G22(vpar[2]*cos(vpar[3]),vpar[2]*sin(vpar[3]));
    EvtComplex H20 = G02*vh[0] + G22*vh[1];
    double hamps2=abs2(H20);
    double xder=(hamps2-hamps0)/dev; 
    fd.push_back(xder);
  }
  return fd;
}