TEvtGen/EvtGenModels/EvtPycont.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: EvtPycont.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 //    RS      October 28, 2002      copied from EvtJscont.cc
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include <stdlib.h>
  24 #include "EvtGenBase/EvtParticle.hh"
  25 #include "EvtGenBase/EvtDecayTable.hh"
  26 #include "EvtGenBase/EvtPDL.hh"
  27 #include "EvtGenModels/EvtPycont.hh"
  28 #include "EvtGenModels/EvtPythia.hh"
  29 #include "EvtGenBase/EvtId.hh"
  30 #include "EvtGenBase/EvtReport.hh"
  31 #include <string>
  32 #include <iostream>
  33 using std::endl;
  34
  35 extern "C" {
  36   extern void pystat_(int &);
  37 }
  38
  39 // extern struct
  40 // {
  41 //   int dc[18];
  42 // } decaych_;
  43
  44
  45 EvtPycont::~EvtPycont() {
  46   //  int i=1;
  47   //  pystat_(i);
  48 }
  49
  50 std::string EvtPycont::getName()
  51 {
  52   return "PYCONT";
  53 }
  54
  55 EvtDecayBase* EvtPycont::clone()
  56 {
  57   return new EvtPycont;
  58 }
  59
  60 void EvtPycont::init()
  61 {
  62   // check that there are 1 argument
  63   if ( getNArg() != 12 && getNArg() != 0 ) {
  64     report(ERROR,"EvtGen") << "EvtPYCONT expects "
  65                            << " 12 arguments (d u s c b t e nu_e mu nu_mu tau nu_tau) but found: "
  66                            << getNArg() <<endl;
  67
  68   }
  69   checkNArg(0,12);
  70
  71    // for( int i=0; i<18; i++)
  72    //   decaych_.dc[i]=0;
  73    // if ( getNArg() == 12 ) {
  74    //   decaych_.dc[0]=(int)getArg(0);
  75    //   decaych_.dc[1]=(int)getArg(1);
  76    //   decaych_.dc[2]=(int)getArg(2);
  77    //   decaych_.dc[3]=(int)getArg(3);
  78    //   decaych_.dc[4]=(int)getArg(4);
  79    //   decaych_.dc[5]=(int)getArg(5);
  80    //   decaych_.dc[10]=(int)getArg(6);
  81    //   decaych_.dc[11]=(int)getArg(7);
  82    //   decaych_.dc[12]=(int)getArg(8);
  83    //   decaych_.dc[13]=(int)getArg(9);
  84    //   decaych_.dc[14]=(int)getArg(10);
  85    //   decaych_.dc[15]=(int)getArg(11);
  86    // }
  87    // else{
  88    //   decaych_.dc[0]=1;
  89    //   decaych_.dc[1]=1;
  90    //   decaych_.dc[2]=1;
  91    //   decaych_.dc[3]=1;
  92    // }
  93
  94 }
  95
  96 void EvtPycont::initProbMax()
  97 {
  98   noProbMax();
  99 }
 100
 101 void EvtPycont::decay( EvtParticle *p)
 102 {
 103   EvtPythia::pythiaInit(0);
 104   EvtVector4R p4[100];
 105
 106   double energy=p->mass();
 107
 108   int i,more;
 109   int ndaugjs;
 110   int kf[100];
 111   EvtId id[100];
 112   int type[MAX_DAUG];
 113
 114   double px[100],py[100],pz[100],e[100];
 115
 116   if ( p->getNDaug() != 0 ) { return;}
 117   do{
 118     EvtPythia::pythiacont(&energy,&ndaugjs,kf,px,py,pz,e);
 119
 120     for(i=0;i<ndaugjs;i++)
 121       {
 122
 123         id[i]=EvtPDL::evtIdFromStdHep(kf[i]);
 124
 125         type[i]=EvtPDL::getSpinType(id[i]);
 126
 127         // have to protect against negative mass^2 for massless particles
 128         // i.e. neutrinos and photons.
 129         // this is uggly but I need to fix it right now....
 130
 131         if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i])
 132           e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001;
 133
 134         p4[i].set(e[i],px[i],py[i],pz[i]);
 135
 136       }
 137
 138     int channel=EvtDecayTable::inChannelList(p->getId(),ndaugjs,id);
 139
 140     more=((channel!=-1)&&(channel!=p->getChannel()));
 141
 142   }while(more);
 143
 144   p->makeDaughters(ndaugjs,id);
 145
 146   for(i=0;i<ndaugjs;i++)
 147     p->getDaug(i)->init( id[i], p4[i] );
 148
 149   return ;
 150 }
 151