[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtJetSetCDF.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: EvtJetSet.cc
//
// Description: Routine to use JetSet for decaying particles.
//
// Modification history:
//
//    RYD     July 24, 1997        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/EvtJetSetCDF.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 EvtJetSetCDF::njetsetdecays=0;
  EvtDecayBasePtr* EvtJetSetCDF::jetsetdecays=0; 
int EvtJetSetCDF::ntable=0;

int EvtJetSetCDF::ncommand=0;
int EvtJetSetCDF::lcommand=0;
std::string* EvtJetSetCDF::commands=0;

extern "C" {
  extern void evtjetsetcdfinit_(char* fname, int len);
}


extern "C" {
  extern void jetsetcdf_(int *,double *,int *,int *,int *,
		       double *,double *,double *,double *);
}

extern "C" {
  extern void lygive_(const char *cnfgstr,int length);
}

extern "C" {
  extern int lycomp_(int* kf);
}


EvtJetSetCDF::EvtJetSetCDF(){}

EvtJetSetCDF::~EvtJetSetCDF(){


  int i;


  //the deletion of commands is really uggly!

  if (njetsetdecays==0) {
    delete [] commands;
    commands=0;
    return;
  }

  for(i=0;i<njetsetdecays;i++){
    if (jetsetdecays[i]==this){
      jetsetdecays[i]=jetsetdecays[njetsetdecays-1];
      njetsetdecays--;
      if (njetsetdecays==0) {
	delete [] commands;
	commands=0;
      }
      return;
    }
  }

  report(ERROR,"EvtGen") << "Error in destroying JetSet model!"<<endl;
 
}


std::string EvtJetSetCDF::getName(){

  return "JETSETCDF";     

}

EvtDecayBase* EvtJetSetCDF::clone(){

  return new EvtJetSetCDF;

}


void EvtJetSetCDF::initProbMax(){

  noProbMax();

}


void EvtJetSetCDF::init(){

  checkNArg(1);


  if (getParentId().isAlias()){

    report(ERROR,"EvtGen") << "EvtJetSetCDF finds that you are decaying the"<<endl
                           << " aliased particle "
			   << EvtPDL::name(getParentId()).c_str()
			   << " with the JetSet model"<<endl
			   << " this does not work, please modify decay table."
			   << endl;
    report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
    ::abort();

  }

  store(this);

}


std::string EvtJetSetCDF::commandName(){

  return std::string("JetSetCDFPar");
  
}


void EvtJetSetCDF::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 EvtJetSetCDF::decay( EvtParticle *p){


  //added by Lange Jan4,2000
  static EvtId STRNG=EvtPDL::getId("string");

  int istdheppar=EvtPDL::getStdHep(p->getId());

  if (lycomp_(&istdheppar)==0){
    report(ERROR,"EvtGen") << "Jetset can not decay:"
      <<EvtPDL::name(p->getId()).c_str()<<endl;
    return;
  }

  double mp=p->mass();

  EvtVector4R p4[20];
  
  int i,more;
  int ip=EvtPDL::getStdHep(p->getId());
  int ndaugjs;
  int kf[100];
  EvtId evtnumstable[100],evtnumparton[100];
  int stableindex[100],partonindex[100];
  int numstable;
  int numparton;
  int km[100];
  EvtId type[MAX_DAUG];

  jetSetInit();

  double px[100],py[100],pz[100],e[100];

  if ( p->getNDaug() != 0 ) { p->deleteDaughters(true);}

  int count=0;

  do{
    //report(INFO,"EvtGen") << "calling jetset " << ip<< " " << mp <<endl;
    jetsetcdf_(&ip,&mp,&ndaugjs,kf,km,px,py,pz,e);
    

    numstable=0;
    numparton=0;
    //report(INFO,"EvtGen") << "found some daughters " << ndaugjs << endl;
    for(i=0;i<ndaugjs;i++){

      if (EvtPDL::evtIdFromStdHep(kf[i])==EvtId(-1,-1)) {
	report(ERROR,"EvtGen") << "JetSet returned particle:"<<kf[i]<<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;

      }

      //sort out the partons
      if (abs(kf[i])<=6||kf[i]==21){
	partonindex[numparton]=i;
	evtnumparton[numparton]=EvtPDL::evtIdFromStdHep(kf[i]);
	numparton++;
      }
      else{
	stableindex[numstable]=i;
	evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(kf[i]); 
	numstable++;
      }


      // have to protect against negative mass^2 for massless particles
      // i.e. neutrinos and photons.
      // this is uggly but I need to fix it right now....

      if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){

        e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001;

      }

      p4[i].set(e[i],px[i],py[i],pz[i]);


    }

    int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);


   more=(channel!=-1);

   
  count++;

  }while( more && (count<10000) );

  if (count>9999) {
    report(INFO,"EvtGen") << "Too many loops in EvtJetSetCDF!!!"<<endl;
    report(INFO,"EvtGen") << "Parent:"<<EvtPDL::name(getParentId()).c_str()<<endl;
    for(i=0;i<numstable;i++){
      report(INFO,"EvtGen") << "Daug("<<i<<")"<<EvtPDL::name(evtnumstable[i]).c_str()<<endl;
    }

  }


  if (numparton==0){

    p->makeDaughters(numstable,evtnumstable);
    int ndaugFound=0;
    for(i=0;i<numstable;i++){
      p->getDaug(i)->init(evtnumstable[i],p4[stableindex[i]]);
      ndaugFound++;
    }
    if ( ndaugFound == 0 ) {
      report(ERROR,"EvtGen") << "Jetset has failed to do a decay ";
      report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass()<<endl;
      assert(0);
    }

    fixPolarizations(p);

    return ;
   
  }
  else{

    //have partons in JETSET

    EvtVector4R p4string(0.0,0.0,0.0,0.0);

    for(i=0;i<numparton;i++){
      p4string+=p4[partonindex[i]];
    }

    int nprimary=1;
    type[0]=STRNG;
    for(i=0;i<numstable;i++){
      if (km[stableindex[i]]==0){
	type[nprimary++]=evtnumstable[i];
      }
    }

    p->makeDaughters(nprimary,type);

    p->getDaug(0)->init(STRNG,p4string);

    EvtVector4R p4partons[10];

    for(i=0;i<numparton;i++){
      p4partons[i]=p4[partonindex[i]];
    }

    ((EvtStringParticle*)p->getDaug(0))->initPartons(numparton,p4partons,evtnumparton);


    nprimary=1;

    for(i=0;i<numstable;i++){

      if (km[stableindex[i]]==0){
	p->getDaug(nprimary++)->init(evtnumstable[i],p4[stableindex[i]]);
      }
    }


    int nsecond=0;
    for(i=0;i<numstable;i++){
      if (km[stableindex[i]]!=0){
	type[nsecond++]=evtnumstable[i];
      }
    }


    p->getDaug(0)->makeDaughters(nsecond,type);

    EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
			      -p4string.get(2),-p4string.get(3));

    nsecond=0;
    for(i=0;i<numstable;i++){
      if (km[stableindex[i]]!=0){
	p4[stableindex[i]]=boostTo(p4[stableindex[i]],p4stringboost);
	p->getDaug(0)->getDaug(nsecond)->init(evtnumstable[i],p4[stableindex[i]]);
	p->getDaug(0)->getDaug(nsecond)->setDiagonalSpinDensity();
	p->getDaug(0)->getDaug(nsecond)->decay();
	nsecond++;
      }
    }

    if ( nsecond == 0 ) {
      report(ERROR,"EvtGen") << "Jetset has failed to do a decay ";
      report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass() <<endl;
      assert(0);
    }

    fixPolarizations(p);

    return ;

  }

}

void EvtJetSetCDF::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");

  for(i=0;i<ndaug;i++){
    if(p->getDaug(i)->getId()==Jpsi){
  
      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 EvtJetSetCDF::store(EvtDecayBase* jsdecay){

  if (njetsetdecays==ntable){

    EvtDecayBasePtr* newjetsetdecays=new EvtDecayBasePtr[2*ntable+10];
    int i;
    for(i=0;i<ntable;i++){
      newjetsetdecays[i]=jetsetdecays[i];
    }
    ntable=2*ntable+10;
    delete [] jetsetdecays;
    jetsetdecays=newjetsetdecays;
  }

  jetsetdecays[njetsetdecays++]=jsdecay;


}


void EvtJetSetCDF::WriteJetSetEntryHeader(ofstream &outdec, int lundkc,
			       EvtId evtnum,std::string name,
			       int chg, int cchg, int spin2,double mass,
			       double width, double maxwidth,double ctau,
			       int stable,double rawbrfrsum){


  char sname[100];

  int namelength=8;

  int i,j;
  int temp;
  temp = spin2;

  if (ctau>1000000.0) ctau=0.0;

  strcpy(sname,name.c_str());

  i=0;

  while (sname[i]!=0){
    i++;
  }

  // strip up to two + or -
 
  if(evtnum.getId()>=0) {
    if (sname[i-1]=='+'||sname[i-1]=='-'){ 
      sname[i-1]=0;
      i--;
    }
    if (sname[i-1]=='+'||sname[i-1]=='-'){ 
      sname[i-1]=0;
      i--;
    }
    // strip 0 except for _0 and chi...0
    if (sname[i-1]=='0' && sname[i-2]!='_' && !(sname[0]=='c' && sname[1]=='h')){
      sname[i-1]=0;
      i--;
    }
  }

  if (i>namelength) {
    for(j=1;j<namelength;j++){
      sname[j]=sname[j+i-namelength];
    }
    sname[namelength]=0;
  }

  cchg=0;

  if(evtnum.getId()>=0) {
    if (abs(EvtPDL::getStdHep(evtnum))==21) cchg=2;
    if (abs(EvtPDL::getStdHep(evtnum))==90) cchg=-1;
    if ((abs(EvtPDL::getStdHep(evtnum))<=8)&&
	(abs(EvtPDL::getStdHep(evtnum))!=0)) cchg=1;

  }

  outdec << setw(5) << lundkc << "  ";
  outdec.width(namelength);
  outdec << setiosflags(ios::left) << sname << resetiosflags(ios::left);
  outdec << setw(3) << chg;
  outdec << setw(3) << cchg;
  outdec.width(3);
  if (evtnum.getId()>=0) {
    if (EvtPDL::chargeConj(evtnum)==evtnum) {
      outdec << 0;
    }
    else{
      outdec << 1;
    }
  }
  else{
    outdec << 0;
  }
  outdec.setf(ios::fixed);
  outdec.precision(5);
  outdec << setw(12) << mass;
  outdec << setw(12) << width;
  outdec.width(12);
  if (fabs(width)<0.0000000001) {
    outdec << 0.0 ;
  }
  else{
    outdec << maxwidth;
  }
  outdec << setw(14) << ctau;
  outdec.width(3);
  if (evtnum.getId()>=0) {
    if (ctau>1.0 || rawbrfrsum<0.000001) {  
      stable=0;
    }
  }
  outdec << stable;
  outdec << endl;
  outdec.width(0);

}

void EvtJetSetCDF::WriteJetSetParticle(ofstream &outdec,EvtId ipar,
				    EvtId iparname,int &first){

  int ijetset;

  double br_sum=0.0;

  for(ijetset=0;ijetset<njetsetdecays;ijetset++){
   
    if (jetsetdecays[ijetset]->getParentId()==ipar){
      br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
    }
    if (jetsetdecays[ijetset]->getParentId()!=
	EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())&&
	EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())==ipar){
      br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
    }


  }

  double br_sum_true=br_sum;

  if (br_sum<0.000001) br_sum=1.0;

  for(ijetset=0;ijetset<njetsetdecays;ijetset++){
    if (jetsetdecays[ijetset]->getParentId()==ipar){

      double br=jetsetdecays[ijetset]->getBranchingFraction();
    
      int i,daugs[5];
      EvtId cdaugs[5];
    
      for(i=0;i<5;i++){
      
	if(i<jetsetdecays[ijetset]->getNDaug()){
	  daugs[i]=EvtPDL::getStdHep(
			 jetsetdecays[ijetset]->getDaugs()[i]);
	  cdaugs[i]=EvtPDL::chargeConj(jetsetdecays[ijetset]->getDaugs()[i]);
	}
	else{
	  daugs[i]=0;
	}
      }

      int channel;

      channel=EvtDecayTable::findChannel(EvtPDL::chargeConj(ipar),
			     jetsetdecays[ijetset]->getModelName(),
			     jetsetdecays[ijetset]->getNDaug(),
			     cdaugs,
			     jetsetdecays[ijetset]->getNArg(),
			     jetsetdecays[ijetset]->getArgsStr());     

      if (jetsetdecays[ijetset]->getModelName()=="JETSET"){
	
	if (first) {
	  first=0;      
	  WriteJetSetEntryHeader(outdec,
				 EvtPDL::getLundKC(iparname),
				 iparname,
				 EvtPDL::name(iparname), 
				 EvtPDL::chg3(iparname),
				 0,0,EvtPDL::getMeanMass(ipar),
				 EvtPDL::getWidth(ipar),
				 EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar),
				 EvtPDL::getctau(ipar),1,br_sum_true);
	}
	
	int dflag=2;

	if (EvtPDL::getStdHep(ipar)<0) {
	  dflag=3;
	  for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
	    daugs[i]=EvtPDL::getStdHep(cdaugs[i]);
	  }

	}

	//now lets check to make sure that jetset, lycomp, knows
        //about all particles!
	int unknown=0;
	for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
	  if (lycomp_(&daugs[i])==0) {
	    unknown=1;
	    report(ERROR,"EvtGen") << "JetSet (lycomp) does not "
				  << "know the particle:"<<
	      EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i]).c_str()<<endl;
	  }
	}

	int istdheppar=EvtPDL::getStdHep(ipar);

	if (lycomp_(&istdheppar)==0){
	  unknown=1;
	  report(ERROR,"EvtGen") << "JetSet (lycomp) does not "
	          << "know the particle:"<<
	      EvtPDL::name(ipar).c_str()<<endl;
	}


	if (unknown){
	  report(ERROR,"EvtGen") << "Therfore the decay:"<<endl;
	  report(ERROR,"EvtGen") << EvtPDL::name(jetsetdecays[ijetset]->getParentId()).c_str()<<" -> ";
	  for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
	    report(ERROR,"") << EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i]).c_str()<<" ";
	  }
	  report(ERROR,"")<<endl;
	  report(ERROR,"EvtGen")<<"Will not be generated."<<endl;
	  return;
	}


	if (EvtPDL::chargeConj(ipar)==ipar) {
	  dflag=1;
	  //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because C(ipar)=ipar!"<<endl;
	}


	//if (channel>=0) {
	//  dflag=1;
	  //report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because channel>=0"<<endl;
	//}
 
	//	if (!(EvtPDL::getStdHep(ipar)<0&&channel>=0)){
	if (1){

	  outdec.width(10);
	  outdec <<dflag;
	  outdec.width(5);
	  outdec <<(int)jetsetdecays[ijetset]->getArgs()[0];
	  outdec.width(12);
	  if (fabs(br)<0.000000001) {
	    outdec <<"0.00000";
	  }
	  else{
	    outdec <<br/br_sum;
	  }
	  outdec.width(8);
	  outdec <<daugs[0];
	  outdec.width(8);
	  outdec <<daugs[1];
	  outdec.width(8);
	  outdec <<daugs[2];
	  outdec.width(8);
	  outdec <<daugs[3];
	  outdec.width(8);
	  outdec <<daugs[4];
	  outdec<<endl;
	  outdec.width(0);
	}
      }
    }
  }
}


void EvtJetSetCDF::MakeJetSetFile(char* fname){
  
  EvtId ipar;
  int lundkc;
  
  //int part_list[MAX_PART];
  
  ofstream outdec;
 
  outdec.open(fname);
  
  //outdec << ";"<<endl;
  //outdec << ";This decayfile has been automatically created by"<<endl;
  //outdec << ";EvtGen from the DECAY.DEC file"<<endl;
  //outdec << ";"<<endl;

  int nokcentry;

  for(lundkc=1;lundkc<=500;lundkc++){

    nokcentry=1;

    for(size_t iipar=0;iipar<EvtPDL::entries();iipar++){

      ipar=EvtId(iipar,iipar);
      //no aliased particles!
      std::string tempStr = EvtPDL::name(ipar);
      EvtId realId = EvtPDL::getId(tempStr);
      if ( realId.isAlias() != 0 ) continue;
      if (lundkc==EvtPDL::getLundKC(ipar)){
        
        nokcentry=0;

	int first=1;
    
	WriteJetSetParticle(outdec,ipar,ipar,first);

	
	EvtId ipar2=EvtPDL::chargeConj(ipar);


        if (ipar2!=ipar){
	  WriteJetSetParticle(outdec,ipar2,ipar,first);
	}

        if (first){
	  WriteJetSetEntryHeader(outdec, 
				    EvtPDL::getLundKC(ipar),
				    ipar,
				    EvtPDL::name(ipar),
				    EvtPDL::chg3(ipar),
				    0,0,EvtPDL::getMeanMass(ipar),
				    EvtPDL::getWidth(ipar),
				    EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar),
				    EvtPDL::getctau(ipar),0,0.0);

	}
      }
    }
    if (nokcentry){

      WriteJetSetEntryHeader(outdec, 
				lundkc,EvtId(-1,-1),"  ",
				0,0,0,EvtPDL::getMeanMass(ipar),0.0,0.0,
				EvtPDL::getctau(ipar),0,0.0);

    }
  }
    outdec.close();
}

void EvtJetSetCDF::jetSetInit(){

  static int first=1;

  if (first){

    first=0;

    report(INFO,"EvtGen") << "Will initialize JetSet."<<endl;

    char fname[200];

    char hostBuffer[100];
    
    if ( gethostname( hostBuffer, 100 ) != 0 ){
      report(ERROR,"EvtGen") << " couldn't get hostname." << endl;
      strncpy( hostBuffer, "hostnameNotFound", 100 );
    }

    char pid[100];

    int thePid=getpid();

    if ( sprintf( pid, "%d", thePid ) == 0 ){
      report(ERROR,"EvtGen") << " couldn't get process ID." << endl;
      strncpy( pid, "666", 100 );
    }

    strcpy(fname,"jet.d-");
    strcat(fname,hostBuffer);
    strcat(fname,"-");
    strcat(fname,pid);
    
    MakeJetSetFile(fname);
    evtjetsetcdfinit_(fname,strlen(fname));

    if (0==getenv("EVTSAVEJETD")){
      char delcmd[300];
      strcpy(delcmd,"rm -f ");
      strcat(delcmd,fname);
      system(delcmd);
    }

    int i;

    for(i=0;i<ncommand;i++){
      lygive_(commands[i].c_str(),strlen(commands[i].c_str()));

    }

    report(INFO,"EvtGen") << "Done initializing JetSet."<<endl;

    //int itest = 5122;
    //report(ERROR,"EvtGen TEST:") << lycomp_(&itest) << endl;


  }

}
Commit	Line	Data
da0e9ce3	1	//--------------------------------------------------------------------------
	2	//
	3	// Environment:
	4	// This software is part of the EvtGen package developed jointly
	5	// for the BaBar and CLEO collaborations. If you use all or part
	6	// of it, please give an appropriate acknowledgement.
	7	//
	8	// Copyright Information: See EvtGen/COPYRIGHT
	9	// Copyright (C) 1998 Caltech, UCSB
	10	//
	11	// Module: EvtJetSet.cc
	12	//
	13	// Description: Routine to use JetSet for decaying particles.
	14	//
	15	// Modification history:
	16	//
	17	// RYD July 24, 1997 Module created
	18	//
	19	//------------------------------------------------------------------------
	20	//
	21	#include "EvtGenBase/EvtPatches.hh"
	22	#include "EvtGenBase/EvtPatches.hh"
	23	#include "EvtGenBase/EvtParticle.hh"
	24	#include "EvtGenBase/EvtStringParticle.hh"
	25	#include "EvtGenBase/EvtDecayTable.hh"
	26	#include "EvtGenBase/EvtPDL.hh"
	27	#include "EvtGenModels/EvtJetSetCDF.hh"
	28	#include "EvtGenBase/EvtReport.hh"
	29	#include <string>
	30	#include "EvtGenBase/EvtId.hh"
	31	#include <iostream>
	32	#include <iomanip>
	33	#include <fstream>
	34	#include <string.h>
	35	#include <stdlib.h>
	36	#include <unistd.h>
	37	#include <stdio.h>
	38	using std::endl;
	39	using std::fstream;
	40	using std::ios;
	41	using std::ofstream;
	42	using std::resetiosflags;
	43	using std::setiosflags;
	44	using std::setw;
	45
	46
	47	int EvtJetSetCDF::njetsetdecays=0;
	48	EvtDecayBasePtr* EvtJetSetCDF::jetsetdecays=0;
	49	int EvtJetSetCDF::ntable=0;
	50
	51	int EvtJetSetCDF::ncommand=0;
	52	int EvtJetSetCDF::lcommand=0;
	53	std::string* EvtJetSetCDF::commands=0;
	54
	55	extern "C" {
	56	extern void evtjetsetcdfinit_(char* fname, int len);
	57	}
	58
	59
	60	extern "C" {
	61	extern void jetsetcdf_(int ,double ,int ,int ,int *,
	62	double ,double ,double ,double );
	63	}
	64
65	extern "C" {
66	extern void lygive_(const char *cnfgstr,int length);
67	}
68
69	extern "C" {
70	extern int lycomp_(int* kf);
71	}
72
73
74	EvtJetSetCDF::EvtJetSetCDF(){}
75
76	EvtJetSetCDF::~EvtJetSetCDF(){
77
78
79	int i;
80
81
82	//the deletion of commands is really uggly!
83
84	if (njetsetdecays==0) {
85	delete [] commands;
86	commands=0;
87	return;
88	}
89
90	for(i=0;i<njetsetdecays;i++){
91	if (jetsetdecays[i]==this){
92	jetsetdecays[i]=jetsetdecays[njetsetdecays-1];
93	njetsetdecays--;
94	if (njetsetdecays==0) {
95	delete [] commands;
96	commands=0;
97	}
98	return;
99	}
100	}
101
102	report(ERROR,"EvtGen") << "Error in destroying JetSet model!"<<endl;
103
104	}
105
106
107	std::string EvtJetSetCDF::getName(){
108
109	return "JETSETCDF";
110
111	}
112
113	EvtDecayBase* EvtJetSetCDF::clone(){
114
115	return new EvtJetSetCDF;
116
117	}
118
119
120	void EvtJetSetCDF::initProbMax(){
121
122	noProbMax();
123
124	}
125
126
127	void EvtJetSetCDF::init(){
128
129	checkNArg(1);
130
131
132	if (getParentId().isAlias()){
133
134	report(ERROR,"EvtGen") << "EvtJetSetCDF finds that you are decaying the"<<endl
135	<< " aliased particle "
136	<< EvtPDL::name(getParentId()).c_str()
137	<< " with the JetSet model"<<endl
138	<< " this does not work, please modify decay table."
139	<< endl;
140	report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
141	::abort();
142
143	}
144
145	store(this);
146
147	}
148
149
150	std::string EvtJetSetCDF::commandName(){
151
152	return std::string("JetSetCDFPar");
153
154	}
155
156
157	void EvtJetSetCDF::command(std::string cmd){
158
159	if (ncommand==lcommand){
160
161	lcommand=10+2*lcommand;
162
163	std::string* newcommands=new std::string[lcommand];
164
165	int i;
166
167	for(i=0;i<ncommand;i++){
168	newcommands[i]=commands[i];
169	}
170
171	delete [] commands;
172
173	commands=newcommands;
174
175	}
176
177	commands[ncommand]=cmd;
178
179	ncommand++;
180
181	}
182
183
184
185	void EvtJetSetCDF::decay( EvtParticle *p){
186
187
188	//added by Lange Jan4,2000
189	static EvtId STRNG=EvtPDL::getId("string");
190
191	int istdheppar=EvtPDL::getStdHep(p->getId());
192
193	if (lycomp_(&istdheppar)==0){
194	report(ERROR,"EvtGen") << "Jetset can not decay:"
195	<<EvtPDL::name(p->getId()).c_str()<<endl;
196	return;
197	}
198
199	double mp=p->mass();
200
201	EvtVector4R p4[20];
202
203	int i,more;
204	int ip=EvtPDL::getStdHep(p->getId());
205	int ndaugjs;
206	int kf[100];
207	EvtId evtnumstable[100],evtnumparton[100];
208	int stableindex[100],partonindex[100];
209	int numstable;
210	int numparton;
211	int km[100];
212	EvtId type[MAX_DAUG];
213
214	jetSetInit();
215
216	double px[100],py[100],pz[100],e[100];
217
218	if ( p->getNDaug() != 0 ) { p->deleteDaughters(true);}
219
220	int count=0;
221
222	do{
223	//report(INFO,"EvtGen") << "calling jetset " << ip<< " " << mp <<endl;
224	jetsetcdf_(&ip,&mp,&ndaugjs,kf,km,px,py,pz,e);
225
226
227	numstable=0;
228	numparton=0;
229	//report(INFO,"EvtGen") << "found some daughters " << ndaugjs << endl;
230	for(i=0;i<ndaugjs;i++){
231
232	if (EvtPDL::evtIdFromStdHep(kf[i])==EvtId(-1,-1)) {
233	report(ERROR,"EvtGen") << "JetSet returned particle:"<<kf[i]<<endl;
234	report(ERROR,"EvtGen") << "This can not be translated to evt number"<<endl;
235	report(ERROR,"EvtGen") << "and the decay will be rejected!"<<endl;
236	report(ERROR,"EvtGen") << "The decay was of particle:"<<ip<<endl;
237
238	}
239
240	//sort out the partons
241	if (abs(kf[i])<=6\|\|kf[i]==21){
242	partonindex[numparton]=i;
243	evtnumparton[numparton]=EvtPDL::evtIdFromStdHep(kf[i]);
244	numparton++;
245	}
246	else{
247	stableindex[numstable]=i;
248	evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(kf[i]);
249	numstable++;
250	}
251
252
253	// have to protect against negative mass^2 for massless particles
254	// i.e. neutrinos and photons.
255	// this is uggly but I need to fix it right now....
256
257	if (px[i]px[i]+py[i]py[i]+pz[i]pz[i]>=e[i]e[i]){
258
259	e[i]=sqrt(px[i]px[i]+py[i]py[i]+pz[i]*pz[i])+0.0000000000001;
260
261	}
262
263	p4[i].set(e[i],px[i],py[i],pz[i]);
264
265
266	}
267
268	int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
269
270
271	more=(channel!=-1);
272
273
274
275
276	count++;
277
278	}while( more && (count<10000) );
279
280	if (count>9999) {
281	report(INFO,"EvtGen") << "Too many loops in EvtJetSetCDF!!!"<<endl;
282	report(INFO,"EvtGen") << "Parent:"<<EvtPDL::name(getParentId()).c_str()<<endl;
283	for(i=0;i<numstable;i++){
284	report(INFO,"EvtGen") << "Daug("<<i<<")"<<EvtPDL::name(evtnumstable[i]).c_str()<<endl;
285	}
286
287	}
288
289
290
291	if (numparton==0){
292
293	p->makeDaughters(numstable,evtnumstable);
294	int ndaugFound=0;
295	for(i=0;i<numstable;i++){
296	p->getDaug(i)->init(evtnumstable[i],p4[stableindex[i]]);
297	ndaugFound++;
298	}
299	if ( ndaugFound == 0 ) {
300	report(ERROR,"EvtGen") << "Jetset has failed to do a decay ";
301	report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass()<<endl;
302	assert(0);
303	}
304
305	fixPolarizations(p);
306
307	return ;
308
309	}
310	else{
311
312	//have partons in JETSET
313
314	EvtVector4R p4string(0.0,0.0,0.0,0.0);
315
316	for(i=0;i<numparton;i++){
317	p4string+=p4[partonindex[i]];
318	}
319
320	int nprimary=1;
321	type[0]=STRNG;
322	for(i=0;i<numstable;i++){
323	if (km[stableindex[i]]==0){
324	type[nprimary++]=evtnumstable[i];
325	}
326	}
327
328	p->makeDaughters(nprimary,type);
329
330	p->getDaug(0)->init(STRNG,p4string);
331
332	EvtVector4R p4partons[10];
333
334	for(i=0;i<numparton;i++){
335	p4partons[i]=p4[partonindex[i]];
336	}
337
338	((EvtStringParticle*)p->getDaug(0))->initPartons(numparton,p4partons,evtnumparton);
339
340
341
342	nprimary=1;
343
344	for(i=0;i<numstable;i++){
345
346	if (km[stableindex[i]]==0){
347	p->getDaug(nprimary++)->init(evtnumstable[i],p4[stableindex[i]]);
348	}
349	}
350
351
352	int nsecond=0;
353	for(i=0;i<numstable;i++){
354	if (km[stableindex[i]]!=0){
355	type[nsecond++]=evtnumstable[i];
356	}
357	}
358
359
360	p->getDaug(0)->makeDaughters(nsecond,type);
361
362	EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
363	-p4string.get(2),-p4string.get(3));
364
365	nsecond=0;
366	for(i=0;i<numstable;i++){
367	if (km[stableindex[i]]!=0){
368	p4[stableindex[i]]=boostTo(p4[stableindex[i]],p4stringboost);
369	p->getDaug(0)->getDaug(nsecond)->init(evtnumstable[i],p4[stableindex[i]]);
370	p->getDaug(0)->getDaug(nsecond)->setDiagonalSpinDensity();
371	p->getDaug(0)->getDaug(nsecond)->decay();
372	nsecond++;
373	}
374	}
375
376	if ( nsecond == 0 ) {
377	report(ERROR,"EvtGen") << "Jetset has failed to do a decay ";
378	report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass() <<endl;
379	assert(0);
380	}
381
382	fixPolarizations(p);
383
384	return ;
385
386	}
387
388	}
389
390	void EvtJetSetCDF::fixPolarizations(EvtParticle *p){
391
392	//special case for now to handle the J/psi polarization
393
394	int ndaug=p->getNDaug();
395
396	int i;
397
398	static EvtId Jpsi=EvtPDL::getId("J/psi");
399
400	for(i=0;i<ndaug;i++){
401	if(p->getDaug(i)->getId()==Jpsi){
402
403	EvtSpinDensity rho;
404
405	rho.setDim(3);
406	rho.set(0,0,0.5);
407	rho.set(0,1,0.0);
408	rho.set(0,2,0.0);
409
410	rho.set(1,0,0.0);
411	rho.set(1,1,1.0);
412	rho.set(1,2,0.0);
413
414	rho.set(2,0,0.0);
415	rho.set(2,1,0.0);
416	rho.set(2,2,0.5);
417
418	EvtVector4R p4Psi=p->getDaug(i)->getP4();
419
420	double alpha=atan2(p4Psi.get(2),p4Psi.get(1));
421	double beta=acos(p4Psi.get(3)/p4Psi.d3mag());
422
423
424	p->getDaug(i)->setSpinDensityForwardHelicityBasis(rho,alpha,beta,0.0);
425	setDaughterSpinDensity(i);
426
427	}
428	}
429
430	}
431
432	void EvtJetSetCDF::store(EvtDecayBase* jsdecay){
433
434	if (njetsetdecays==ntable){
435
436	EvtDecayBasePtr* newjetsetdecays=new EvtDecayBasePtr[2*ntable+10];
437	int i;
438	for(i=0;i<ntable;i++){
439	newjetsetdecays[i]=jetsetdecays[i];
440	}
441	ntable=2*ntable+10;
442	delete [] jetsetdecays;
443	jetsetdecays=newjetsetdecays;
444	}
445
446	jetsetdecays[njetsetdecays++]=jsdecay;
447
448
449
450	}
451
452
453	void EvtJetSetCDF::WriteJetSetEntryHeader(ofstream &outdec, int lundkc,
454	EvtId evtnum,std::string name,
455	int chg, int cchg, int spin2,double mass,
456	double width, double maxwidth,double ctau,
457	int stable,double rawbrfrsum){
458
459
460	char sname[100];
461
462	int namelength=8;
463
464	int i,j;
465	int temp;
466	temp = spin2;
467
468	if (ctau>1000000.0) ctau=0.0;
469
470	strcpy(sname,name.c_str());
471
472	i=0;
473
474	while (sname[i]!=0){
475	i++;
476	}
477
478	// strip up to two + or -
479
480	if(evtnum.getId()>=0) {
481	if (sname[i-1]=='+'\|\|sname[i-1]=='-'){
482	sname[i-1]=0;
483	i--;
484	}
485	if (sname[i-1]=='+'\|\|sname[i-1]=='-'){
486	sname[i-1]=0;
487	i--;
488	}
489	// strip 0 except for _0 and chi...0
490	if (sname[i-1]=='0' && sname[i-2]!='_' && !(sname[0]=='c' && sname[1]=='h')){
491	sname[i-1]=0;
492	i--;
493	}
494	}
495
496	if (i>namelength) {
497	for(j=1;j<namelength;j++){
498	sname[j]=sname[j+i-namelength];
499	}
500	sname[namelength]=0;
501	}
502
503	cchg=0;
504
505	if(evtnum.getId()>=0) {
506	if (abs(EvtPDL::getStdHep(evtnum))==21) cchg=2;
507	if (abs(EvtPDL::getStdHep(evtnum))==90) cchg=-1;
508	if ((abs(EvtPDL::getStdHep(evtnum))<=8)&&
509	(abs(EvtPDL::getStdHep(evtnum))!=0)) cchg=1;
510
511	}
512
513	outdec << setw(5) << lundkc << " ";
514	outdec.width(namelength);
515	outdec << setiosflags(ios::left) << sname << resetiosflags(ios::left);
516	outdec << setw(3) << chg;
517	outdec << setw(3) << cchg;
518	outdec.width(3);
519	if (evtnum.getId()>=0) {
520	if (EvtPDL::chargeConj(evtnum)==evtnum) {
521	outdec << 0;
522	}
523	else{
524	outdec << 1;
525	}
526	}
527	else{
528	outdec << 0;
529	}
530	outdec.setf(ios::fixed);
531	outdec.precision(5);
532	outdec << setw(12) << mass;
533	outdec << setw(12) << width;
534	outdec.width(12);
535	if (fabs(width)<0.0000000001) {
536	outdec << 0.0 ;
537	}
538	else{
539	outdec << maxwidth;
540	}
541	outdec << setw(14) << ctau;
542	outdec.width(3);
543	if (evtnum.getId()>=0) {
544	if (ctau>1.0 \|\| rawbrfrsum<0.000001) {
545	stable=0;
546	}
547	}
548	outdec << stable;
549	outdec << endl;
550	outdec.width(0);
551
552	}
553
554	void EvtJetSetCDF::WriteJetSetParticle(ofstream &outdec,EvtId ipar,
555	EvtId iparname,int &first){
556
557	int ijetset;
558
559	double br_sum=0.0;
560
561	for(ijetset=0;ijetset<njetsetdecays;ijetset++){
562
563	if (jetsetdecays[ijetset]->getParentId()==ipar){
564	br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
565	}
566	if (jetsetdecays[ijetset]->getParentId()!=
567	EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())&&
568	EvtPDL::chargeConj(jetsetdecays[ijetset]->getParentId())==ipar){
569	br_sum+=jetsetdecays[ijetset]->getBranchingFraction();
570	}
571
572
573	}
574
575	double br_sum_true=br_sum;
576
577	if (br_sum<0.000001) br_sum=1.0;
578
579	for(ijetset=0;ijetset<njetsetdecays;ijetset++){
580	if (jetsetdecays[ijetset]->getParentId()==ipar){
581
582	double br=jetsetdecays[ijetset]->getBranchingFraction();
583
584	int i,daugs[5];
585	EvtId cdaugs[5];
586
587	for(i=0;i<5;i++){
588
589	if(i<jetsetdecays[ijetset]->getNDaug()){
590	daugs[i]=EvtPDL::getStdHep(
591	jetsetdecays[ijetset]->getDaugs()[i]);
592	cdaugs[i]=EvtPDL::chargeConj(jetsetdecays[ijetset]->getDaugs()[i]);
593	}
594	else{
595	daugs[i]=0;
596	}
597	}
598
599	int channel;
600
601	channel=EvtDecayTable::findChannel(EvtPDL::chargeConj(ipar),
602	jetsetdecays[ijetset]->getModelName(),
603	jetsetdecays[ijetset]->getNDaug(),
604	cdaugs,
605	jetsetdecays[ijetset]->getNArg(),
606	jetsetdecays[ijetset]->getArgsStr());
607
608	if (jetsetdecays[ijetset]->getModelName()=="JETSET"){
609
610	if (first) {
611	first=0;
612	WriteJetSetEntryHeader(outdec,
613	EvtPDL::getLundKC(iparname),
614	iparname,
615	EvtPDL::name(iparname),
616	EvtPDL::chg3(iparname),
617	0,0,EvtPDL::getMeanMass(ipar),
618	EvtPDL::getWidth(ipar),
619	EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar),
620	EvtPDL::getctau(ipar),1,br_sum_true);
621	}
622
623	int dflag=2;
624
625	if (EvtPDL::getStdHep(ipar)<0) {
626	dflag=3;
627	for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
628	daugs[i]=EvtPDL::getStdHep(cdaugs[i]);
629	}
630
631	}
632
633	//now lets check to make sure that jetset, lycomp, knows
634	//about all particles!
635	int unknown=0;
636	for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
637	if (lycomp_(&daugs[i])==0) {
638	unknown=1;
639	report(ERROR,"EvtGen") << "JetSet (lycomp) does not "
640	<< "know the particle:"<<
641	EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i]).c_str()<<endl;
642	}
643	}
644
645	int istdheppar=EvtPDL::getStdHep(ipar);
646
647	if (lycomp_(&istdheppar)==0){
648	unknown=1;
649	report(ERROR,"EvtGen") << "JetSet (lycomp) does not "
650	<< "know the particle:"<<
651	EvtPDL::name(ipar).c_str()<<endl;
652	}
653
654
655
656	if (unknown){
657	report(ERROR,"EvtGen") << "Therfore the decay:"<<endl;
658	report(ERROR,"EvtGen") << EvtPDL::name(jetsetdecays[ijetset]->getParentId()).c_str()<<" -> ";
659	for(i=0;i<jetsetdecays[ijetset]->getNDaug();i++){
660	report(ERROR,"") << EvtPDL::name(jetsetdecays[ijetset]->getDaugs()[i]).c_str()<<" ";
661	}
662	report(ERROR,"")<<endl;
663	report(ERROR,"EvtGen")<<"Will not be generated."<<endl;
664	return;
665	}
666
667
668	if (EvtPDL::chargeConj(ipar)==ipar) {
669	dflag=1;
670	//report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because C(ipar)=ipar!"<<endl;
671	}
672
673
674	//if (channel>=0) {
675	// dflag=1;
676	//report(INFO,"EvtGen") << EvtPDL::name(iparname) << " dflag=1 because channel>=0"<<endl;
677	//}
678
679	// if (!(EvtPDL::getStdHep(ipar)<0&&channel>=0)){
680	if (1){
681
682	outdec.width(10);
683	outdec <<dflag;
684	outdec.width(5);
685	outdec <<(int)jetsetdecays[ijetset]->getArgs()[0];
686	outdec.width(12);
687	if (fabs(br)<0.000000001) {
688	outdec <<"0.00000";
689	}
690	else{
691	outdec <<br/br_sum;
692	}
693	outdec.width(8);
694	outdec <<daugs[0];
695	outdec.width(8);
696	outdec <<daugs[1];
697	outdec.width(8);
698	outdec <<daugs[2];
699	outdec.width(8);
700	outdec <<daugs[3];
701	outdec.width(8);
702	outdec <<daugs[4];
703	outdec<<endl;
704	outdec.width(0);
705	}
706	}
707	}
708	}
709	}
710
711
712	void EvtJetSetCDF::MakeJetSetFile(char* fname){
713
714	EvtId ipar;
715	int lundkc;
716
717	//int part_list[MAX_PART];
718
719	ofstream outdec;
720
721	outdec.open(fname);
722
723	//outdec << ";"<<endl;
724	//outdec << ";This decayfile has been automatically created by"<<endl;
725	//outdec << ";EvtGen from the DECAY.DEC file"<<endl;
726	//outdec << ";"<<endl;
727
728	int nokcentry;
729
730	for(lundkc=1;lundkc<=500;lundkc++){
731
732	nokcentry=1;
733
734	for(size_t iipar=0;iipar<EvtPDL::entries();iipar++){
735
736	ipar=EvtId(iipar,iipar);
737	//no aliased particles!
738	std::string tempStr = EvtPDL::name(ipar);
739	EvtId realId = EvtPDL::getId(tempStr);
740	if ( realId.isAlias() != 0 ) continue;
741	if (lundkc==EvtPDL::getLundKC(ipar)){
742
743	nokcentry=0;
744
745	int first=1;
746
747	WriteJetSetParticle(outdec,ipar,ipar,first);
748
749
750	EvtId ipar2=EvtPDL::chargeConj(ipar);
751
752
753	if (ipar2!=ipar){
754	WriteJetSetParticle(outdec,ipar2,ipar,first);
755	}
756
757	if (first){
758	WriteJetSetEntryHeader(outdec,
759	EvtPDL::getLundKC(ipar),
760	ipar,
761	EvtPDL::name(ipar),
762	EvtPDL::chg3(ipar),
763	0,0,EvtPDL::getMeanMass(ipar),
764	EvtPDL::getWidth(ipar),
765	EvtPDL::getMeanMass(ipar)-EvtPDL::getMinMass(ipar),
766	EvtPDL::getctau(ipar),0,0.0);
767
768	}
769	}
770	}
771	if (nokcentry){
772
773	WriteJetSetEntryHeader(outdec,
774	lundkc,EvtId(-1,-1)," ",
775	0,0,0,EvtPDL::getMeanMass(ipar),0.0,0.0,
776	EvtPDL::getctau(ipar),0,0.0);
777
778	}
779	}
780	outdec.close();
781	}
782
783	void EvtJetSetCDF::jetSetInit(){
784
785	static int first=1;
786
787	if (first){
788
789	first=0;
790
791	report(INFO,"EvtGen") << "Will initialize JetSet."<<endl;
792
793	char fname[200];
794
795	char hostBuffer[100];
796
797	if ( gethostname( hostBuffer, 100 ) != 0 ){
798	report(ERROR,"EvtGen") << " couldn't get hostname." << endl;
799	strncpy( hostBuffer, "hostnameNotFound", 100 );
800	}
801
802	char pid[100];
803
804	int thePid=getpid();
805
806	if ( sprintf( pid, "%d", thePid ) == 0 ){
807	report(ERROR,"EvtGen") << " couldn't get process ID." << endl;
808	strncpy( pid, "666", 100 );
809	}
810
811	strcpy(fname,"jet.d-");
812	strcat(fname,hostBuffer);
813	strcat(fname,"-");
814	strcat(fname,pid);
815
816	MakeJetSetFile(fname);
817	evtjetsetcdfinit_(fname,strlen(fname));
818
819	if (0==getenv("EVTSAVEJETD")){
820	char delcmd[300];
821	strcpy(delcmd,"rm -f ");
822	strcat(delcmd,fname);
823	system(delcmd);
824	}
825
826	int i;
827
828	for(i=0;i<ncommand;i++){
829	lygive_(commands[i].c_str(),strlen(commands[i].c_str()));
830
831	}
832
833	report(INFO,"EvtGen") << "Done initializing JetSet."<<endl;
834
835	//int itest = 5122;
836	//report(ERROR,"EvtGen TEST:") << lycomp_(&itest) << endl;
837
838
839	}
840
841	}