EvtOpenCharm.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of models developed at BES collaboration
//      based on the EvtGen framework.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See CVS repository
//      Copyright (A) 2011      Ping Rong-Gang
//
// Module: EvtOpenCharm.cc
//
// Description: Routine to decay charmonium-> DD, DDpi according the
// cross section measurement by CLEO  PRD 80, 072001.
//
// Modification history:
//
//    Ping R.-G.    December, 2011       Module created
//
//------------------------------------------------------------------------
//
 
#include "EvtOpenCharm.hh"
#include "../EvtGenBase/EvtDecayTable.hh"
#include "../EvtGenBase/EvtId.hh"
#include "../EvtGenBase/EvtPDL.hh"
#include "../EvtGenBase/EvtParticle.hh"
#include "../EvtGenBase/EvtPatches.hh"
#include "../EvtGenBase/EvtReport.hh"
#include "../EvtGenBase/EvtStringParticle.hh"
#include "EvtPsi3Sdecay.hh"
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <string>
#include <unistd.h>
using std::endl;
using std::fstream;
using std::ios;
using std::ofstream;
using std::resetiosflags;
using std::setiosflags;
using std::setw;
 
int              EvtOpenCharm::nOpencharmdecays = 0;
EvtDecayBasePtr* EvtOpenCharm::Opencharmdecays  = 0;
int              EvtOpenCharm::ntable           = 0;
 
int          EvtOpenCharm::ncommand = 0;
int          EvtOpenCharm::lcommand = 0;
std::string* EvtOpenCharm::commands = 0;
int          EvtOpenCharm::nevt     = 0;
 
int                 EvtOpenCharm::myiter = 1;
std::vector<EvtId>  EvtOpenCharm::mypar;
std::vector<int>    EvtOpenCharm::vmode;
std::vector<double> EvtOpenCharm::Vcms;
 
EvtOpenCharm::EvtOpenCharm() {}
 
EvtOpenCharm::~EvtOpenCharm() {
 
  int i;
 
  // the deletion of commands is really uggly!
 
  if ( nOpencharmdecays == 0 )
  {
    delete[] commands;
    commands = 0;
    return;
  }
 
  for ( i = 0; i < nOpencharmdecays; i++ )
  {
    if ( Opencharmdecays[i] == this )
    {
      Opencharmdecays[i] = Opencharmdecays[nOpencharmdecays - 1];
      nOpencharmdecays--;
      if ( nOpencharmdecays == 0 )
      {
        delete[] commands;
        commands = 0;
      }
      return;
    }
  }
 
  report( ERROR, "EvtGen" ) << "Error in destroying OpenCharm model!" << endl;
}
 
void EvtOpenCharm::getName( std::string& model_name ) { model_name = "OPENCHARM"; }
 
EvtDecayBase* EvtOpenCharm::clone() { return new EvtOpenCharm; }
 
void EvtOpenCharm::initProbMax() { noProbMax(); }
 
void EvtOpenCharm::init() {
 
  //  checkNArg(1);
 
  if ( getParentId().isAlias() )
  {
 
    report( ERROR, "EvtGen" ) << "EvtOpenCharm finds that you are decaying the" << endl
                              << " aliased particle " << EvtPDL::name( getParentId() ).c_str()
                              << " with the OpenCharm model" << endl
                              << " this does not work, please modify decay table." << endl;
    report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
    ::abort();
  }
 
  store( this );
}
 
std::string EvtOpenCharm::commandName() { return std::string( "OpenCharmPar" ); }
 
void EvtOpenCharm::command( std::string cmd ) {
 
  if ( ncommand == lcommand )
  {
 
    lcommand = 10 + 2 * lcommand;
 
    std::string* newcommands = new std::string[lcommand];
 
    int i;
 
    for ( i = 0; i < ncommand; i++ ) { newcommands[i] = commands[i]; }
 
    delete[] commands;
 
    commands = newcommands;
  }
 
  commands[ncommand] = cmd;
 
  ncommand++;
}
 
void EvtOpenCharm::decay( EvtParticle* p ) {
 
  int istdheppar = EvtPDL::getStdHep( p->getId() );
 
  if ( istdheppar != 9000443 && istdheppar != 9010443 && istdheppar != 9030443 &&
       istdheppar != 9020443 )
  {
    std::cout << "EvtGen: EvtOpenCharm cann't not decay the particle pid= " << istdheppar
              << endl;
    ::abort();
  }
 
  double mp    = p->mass();
  float  xmp   = mp;
  double totEn = 0;
 
  // debugging
  //   std::cout<<"parent "<<EvtPDL::name(p->getId())<<"float xmp="<<xmp<<"
  //   "<<p->getP4Lab()<<std::endl;
 
  EvtVector4R p4[20];
 
  int i, more;
  int ip = EvtPDL::getStdHep( p->getId() );
  int ndaugjs;
 
  static int    myflag;
  EvtPsi3Sdecay theIni;
  EvtId         pid = p->getId();
 
  static int themode;
  if ( getNArg() == 1 )
  {
    themode = getArg( 0 );
    theIni.setMode( themode );
  }
 
  int count = 0;
  do {
 
    theIni.PHSPDecay( p );
    std::vector<EvtVector4R> v_p4 = theIni.getDaugP4();
    std::vector<EvtId>       Vid  = theIni.getDaugId();
    ndaugjs                       = Vid.size();
 
    EvtId myId[3];
 
    for ( int i = 0; i < ndaugjs; i++ ) { myId[i] = Vid[i]; }
 
    if ( p->getNDaug() != 0 ) p->resetNDaug();
    p->makeDaughters( ndaugjs, myId );
 
    for ( int i = 0; i < ndaugjs; i++ )
    {
      // std::cout<<"Vid "<<EvtPDL::name(Vid[i])<<" p4: "<<v_p4[i]<<std::endl;
      p->getDaug( i )->init( Vid[i], v_p4[i] );
    }
 
    theMode = theIni.getMode();
    p->setGeneratorFlag( theMode );
 
    totEn = 0;
    for ( i = 0; i < ndaugjs; i++ )
    {
      totEn += p->getDaug( i )->getP4().get( 0 );
      if ( p->getDaug( i )->getId() == EvtId( -1, -1 ) )
      {
        report( ERROR, "EvtGen" )
            << "OpenCharm returned particle:" << EvtPDL::name( p->getDaug( i )->getId() )
            << endl;
        report( ERROR, "EvtGen" ) << "This can not be translated to evt number" << endl;
        report( ERROR, "EvtGen" ) << "and the decay will be rejected!" << endl;
        report( ERROR, "EvtGen" ) << "The decay was of particle:" << ip << endl;
      }
    }
    int channel = EvtDecayTable::inChannelList( p->getId(), ndaugjs, myId );
 
    // std::cout<<"channel= "<<channel<<std::endl;
    more = ( channel != -1 );
    // if(more){std::cout<<"count= "<<count<<" inchannel= "<<channel<<std::endl;}
 
    count++;
 
  } while ( more && ( count < 10000 ) );
 
  if ( fabs( xmp - totEn ) > 0.01 )
  {
    std::cout << "Warning:OPENCHARM generate incomplet final state, " << mp << " " << totEn
              << endl;
    ::abort();
  }
 
  if ( count > 9999 )
  {
    report( INFO, "EvtGen" ) << "Too many loops in EvtOpenCharm!!!" << endl;
    report( INFO, "EvtGen" ) << "Parent:" << EvtPDL::name( getParentId() ).c_str() << endl;
  }
 
  fixPolarizations( p );
 
  return;
}
 
void EvtOpenCharm::fixPolarizations( EvtParticle* p ) {
 
  // special case for now to handle the J/psi polarization
 
  int ndaug = p->getNDaug();
 
  int i;
 
  static EvtId Jpsi  = EvtPDL::getId( "J/psi" );
  static EvtId psip  = EvtPDL::getId( "psi(2S)" );
  static EvtId psipp = EvtPDL::getId( "psi(3770)" );
  static EvtId psi_a = EvtPDL::getId( "psi(4040)" );
  static EvtId psi_b = EvtPDL::getId( "psi(4160)" );
  static EvtId psi_c = EvtPDL::getId( "psi(4260)" );
 
  for ( i = 0; i < ndaug; i++ )
  {
    EvtId idp   = p->getDaug( i )->getId();
    bool  bflag = ( idp == Jpsi || idp == psip || idp == psipp || idp == psi_a ||
                   idp == psi_b || idp == psi_c );
    if ( bflag )
    {
 
      EvtSpinDensity rho;
 
      rho.SetDim( 3 );
      rho.Set( 0, 0, 0.5 );
      rho.Set( 0, 1, 0.0 );
      rho.Set( 0, 2, 0.0 );
 
      rho.Set( 1, 0, 0.0 );
      rho.Set( 1, 1, 1.0 );
      rho.Set( 1, 2, 0.0 );
 
      rho.Set( 2, 0, 0.0 );
      rho.Set( 2, 1, 0.0 );
      rho.Set( 2, 2, 0.5 );
 
      EvtVector4R p4Psi = p->getDaug( i )->getP4();
 
      double alpha = atan2( p4Psi.get( 2 ), p4Psi.get( 1 ) );
      double beta  = acos( p4Psi.get( 3 ) / p4Psi.d3mag() );
 
      p->getDaug( i )->setSpinDensityForwardHelicityBasis( rho, alpha, beta, 0.0 );
      setDaughterSpinDensity( i );
    }
  }
}
 
void EvtOpenCharm::store( EvtDecayBase* jsdecay ) {
 
  if ( nOpencharmdecays == ntable )
  {
 
    EvtDecayBasePtr* newOpencharmdecays = new EvtDecayBasePtr[2 * ntable + 10];
    int              i;
    for ( i = 0; i < ntable; i++ ) { newOpencharmdecays[i] = Opencharmdecays[i]; }
    ntable = 2 * ntable + 10;
    delete[] Opencharmdecays;
    Opencharmdecays = newOpencharmdecays;
  }
 
  Opencharmdecays[nOpencharmdecays++] = jsdecay;
}
 
void EvtOpenCharm::OpencrmInit( int dummy ) {}
 
bool EvtOpenCharm::isbelong( EvtId myid ) {
  for ( int i = 0; i < mypar.size(); i++ )
  {
    if ( myid == mypar[i] )
    {
      _index = i;
      return true;
    }
  }
  return false;
}
 
int EvtOpenCharm::which_mode( EvtId myid ) {
  for ( int i = 0; i < mypar.size(); i++ )
  {
    if ( myid == mypar[i] )
    {
      _index = i;
      return i;
    }
  }
  std::cout << "EvtOpenCharm::which_mode() fails to find element" << std::endl;
  abort();
}