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