/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Generator/BesEvtGen/BesEvtGen-00-03-58/src/EvtGen/EvtGenModels/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 "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtStringParticle.hh"
#include "EvtGenBase/EvtDecayTable.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenModels/EvtOpenCharm.hh"
#include "EvtGenModels/EvtPsi3Sdecay.hh"
#include "EvtGenBase/EvtReport.hh"
#include <string>
#include "EvtGenBase/EvtId.hh"
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <stdio.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();
}
  

---