TEvtGen/EvtGenModels/EvtPyGaGa.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/EvtDecayBase.hh"
  27 #include "EvtGenBase/EvtPDL.hh"
  28 #include "EvtGenModels/EvtPyGaGa.hh"
  29 #include "EvtGenModels/EvtPythia.hh"
  30 #include "EvtGenBase/EvtId.hh"
  31 #include "EvtGenBase/EvtReport.hh"
  32 #include <string.h>
  33 #include <iostream>
  34
  35 extern "C" {
  36   extern void pystat_(int &);
  37 }
  38
  39 // extern struct
  40 // {
  41 //   int dc[18];
  42 // } decaych_;
  43
  44 EvtPyGaGa::~EvtPyGaGa()
  45 {
  46   int i=1;
  47   pystat_(i);
  48 }
  49
  50 std::string EvtPyGaGa::getName()
  51 {
  52   return "PYGAGA";
  53 }
  54
  55 EvtDecayBase* EvtPyGaGa::clone()
  56 {
  57   return new EvtPyGaGa;
  58 }
  59
  60 void EvtPyGaGa::init()
  61 {
  62   // check that there are 1 argument
  63   checkNArg(0);
  64   // for( int i=0; i<18; i++)
  65   //   decaych_.dc[i]=0;
  66 }
  67
  68 void EvtPyGaGa::initProbMax()
  69 {
  70   noProbMax();
  71 }
  72
  73 void EvtPyGaGa::decay( EvtParticle *p)
  74 {
  75   EvtPythia::pythiaInit(1);
  76   EvtVector4R p4[100];
  77
  78   double energy=p->mass();
  79
  80   int i,more;
  81   int ndaugjs;
  82   int kf[100];
  83   EvtId id[100];
  84   int type[MAX_DAUG];
  85
  86   double px[100],py[100],pz[100],e[100];
  87
  88   if ( p->getNDaug() != 0 ) { return;}
  89   do{
  90     EvtPythia::pythiacont(&energy,&ndaugjs,kf,px,py,pz,e);
  91
  92     for(i=0;i<ndaugjs;i++)
  93       {
  94
  95         id[i]=EvtPDL::evtIdFromStdHep(kf[i]);
  96
  97         type[i]=EvtPDL::getSpinType(id[i]);
  98
  99         // have to protect against negative mass^2 for massless particles
 100         // i.e. neutrinos and photons.
 101         // this is uggly but I need to fix it right now....
 102
 103         if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i])
 104           e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000000001;
 105
 106         p4[i].set(e[i],px[i],py[i],pz[i]);
 107
 108       }
 109
 110     int channel=EvtDecayTable::inChannelList(p->getId(),ndaugjs,id);
 111
 112     more=((channel!=-1)&&(channel!=p->getChannel()));
 113
 114   }while(more);
 115
 116   p->makeDaughters(ndaugjs,id);
 117
 118   for(i=0;i<ndaugjs;i++)
 119     p->getDaug(i)->init( id[i], p4[i] );
 120
 121   return ;
 122 }
 123