TEvtGen/EvtGenModels/EvtJscontCDF.cxx

   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: EvtJscontCDF.cc
  12 //
  13 // Description: Routine to generate e+e- --> q\barq  via Jetset
  14 //
  15 // Modification history:
  16 //
  17 //    PCK     August 4, 1997        Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <stdlib.h>
  23 #include <stdio.h>
  24 #include "EvtGenBase/EvtParticle.hh"
  25 #include "EvtGenBase/EvtDecayTable.hh"
  26 #include "EvtGenBase/EvtPDL.hh"
  27 #include "EvtGenModels/EvtJscontCDF.hh"
  28 #include "EvtGenModels/EvtJetSetCDF.hh"
  29 #include "EvtGenBase/EvtId.hh"
  30 #include "EvtGenBase/EvtReport.hh"
  31 #include <string>
  32 using namespace std;
  33
  34 extern "C" {
  35   extern void continuumcdf_(double *,int *,int *,int *,
  36                        double *,double *,double *,double *);
  37 }
  38
  39 extern "C" {
  40   extern void lygive_(const char *cnfgstr,int length);
  41 }
  42
  43
  44 EvtJscontCDF::~EvtJscontCDF() {}
  45
  46 std::string EvtJscontCDF::getName(){
  47
  48   return "JSCONTCDF";
  49
  50 }
  51
  52 EvtDecayBase* EvtJscontCDF::clone(){
  53
  54   return new EvtJscontCDF;
  55
  56 }
  57
  58 void EvtJscontCDF::init(){
  59
  60   // check that there are 1 argument
  61
  62   checkNArg(1,2);
  63
  64 }
  65
  66
  67 void EvtJscontCDF::initProbMax(){
  68
  69   noProbMax();
  70
  71 }
  72
  73
  74 void EvtJscontCDF::decay( EvtParticle *p){
  75
  76   EvtJetSetCDF::jetSetInit();
  77   static int first=1;
  78
  79   if (first){
  80     first=0;
  81
  82     float val=0.6;
  83     if ( getNArg()>1) {
  84       val=getArg(1);
  85     }
  86     char vak[20];
  87     sprintf(vak,"PARJ(13)=%f",val);
  88     std::string temp(vak);
  89     lygive_(temp.c_str(),strlen(temp.c_str()));
  90   }
  91   EvtVector4R p4[100];
  92
  93   double energy=p->mass();
  94
  95   int flavor;
  96
  97   int i,more;
  98   int ndaugjs;
  99   int kf[100];
 100   EvtId id[100];
 101   int type[MAX_DAUG];
 102
 103   flavor=(int)getArg(0);
 104
 105   double px[100],py[100],pz[100],e[100];
 106
 107   if ( p->getNDaug() != 0 ) { return;}
 108   do{
 109
 110     continuumcdf_(&energy,&flavor,&ndaugjs,kf,px,py,pz,e);
 111
 112     for(i=0;i<ndaugjs;i++){
 113
 114       id[i]=EvtPDL::evtIdFromStdHep(kf[i]);
 115
 116       type[i]=EvtPDL::getSpinType(id[i]);
 117
 118       // have to protect against negative mass^2 for massless particles
 119       // i.e. neutrinos and photons.
 120       // this is uggly but I need to fix it right now....
 121
 122       if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
 123
 124         e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001;
 125
 126       }
 127
 128       p4[i].set(e[i],px[i],py[i],pz[i]);
 129
 130     }
 131
 132     int channel=EvtDecayTable::inChannelList(p->getId(),ndaugjs,id);
 133
 134     more=((channel!=-1)&&(channel!=p->getChannel()));
 135
 136
 137   }while(more);
 138
 139   p->makeDaughters(ndaugjs,id);
 140
 141   for(i=0;i<ndaugjs;i++){
 142     p->getDaug(i)->init( id[i], p4[i] );
 143   }
 144   return ;
 145 }
 146
 147
 148
 149