diff --git a/GeneratorParam/TMCParticle.cxx b/GeneratorParam/TMCParticle.cxx new file mode 100644 index 0000000..425bffa --- /dev/null +++ b/GeneratorParam/TMCParticle.cxx @@ -0,0 +1,52 @@ +// @(#)root/pythia6:$Id$ +// Author: Piotr Golonka 17/09/97 + +/** \class TMCParticle + \ingroup pythia6 + +This class serves as a data storage for description of one particle. + +It is especially convenient to store information taken from LUJETS common, +which is done by interface class TPythia. + +Author: Piotr Golonka 17/09/97 +*/ + +#include "TMCParticle.h" +#include "TPrimary.h" + +#ifndef WIN32 +# define pyname pyname_ +extern "C" void pyname(const Int_t &kf, const char *name, const Int_t len); +#else +# define pyname PYNAME +extern "C" void pyname(const Int_t &kf, const char *name, const Int_t len); +#endif + +ClassImp(TMCParticle); + +//////////////////////////////////////////////////////////////////////////////// + +void TMCParticle::ls(Option_t *) const +{ + printf("(%2i,%4i) <-%3i, =>[%3i,%3i]",fKS,fKF,fParent, + fFirstChild,fLastChild); + printf(": p=(%7.3f,%7.3f,%9.3f) ;",fPx,fPy,fPz); + + printf(" E=%8.3f ; m=%7.3f ; V=(%g,%g,%g); t=%g, tau=%g\n", + fEnergy,fMass,fVx,fVy,fVz,fTime,fLifetime); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Return name of this particle via Pythia + +const char *TMCParticle::GetName() const +{ + static char name[20]; + pyname(fKF,name,16); name[15] = 0; + for (Int_t i=14;i>0;i--) { + if (name[i] != ' ') break; + name[i] = 0; + } + return name; +} diff --git a/GeneratorParam/TMCParticle.h b/GeneratorParam/TMCParticle.h new file mode 100644 index 0000000..34f2fd3 --- /dev/null +++ b/GeneratorParam/TMCParticle.h @@ -0,0 +1,117 @@ +// @(#)root/pythia6:$Id$ +// Author: Piotr Golonka 17/09/97 + +/************************************************************************* + * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +#ifndef ROOT_TMCParticle +#define ROOT_TMCParticle + +#include "TObject.h" +#include "TAttLine.h" +#include "TPrimary.h" + + +class TMCParticle : public TObject, public TAttLine { + +private: + + Int_t fKS; // status of particle ( LUJETS K[1] ) + Int_t fKF; // KF flavour code ( LUJETS K[2] ) + Int_t fParent; // parrent's id ( LUJETS K[3] ) + Int_t fFirstChild; // id of first child ( LUJETS K[4] ) + Int_t fLastChild; // id of last child ( LUJETS K[5] ) + + Float_t fPx; // X momenta [GeV/c] ( LUJETS P[1] ) + Float_t fPy; // Y momenta [GeV/c] ( LUJETS P[2] ) + Float_t fPz; // Z momenta [GeV/c] ( LUJETS P[3] ) + Float_t fEnergy; // Energy [GeV] ( LUJETS P[4] ) + Float_t fMass; // Mass [Gev/c^2] ( LUJETS P[5] ) + + Float_t fVx; // X vertex [mm] ( LUJETS V[1] ) + Float_t fVy; // Y vertex [mm] ( LUJETS V[2] ) + Float_t fVz; // Z vertex [mm] ( LUJETS V[3] ) + Float_t fTime; // time of procuction [mm/c]( LUJETS V[4] ) + Float_t fLifetime; // proper lifetime [mm/c] ( LUJETS V[5] ) + + +public: + TMCParticle() : fKS(0), fKF(0), fParent(0), fFirstChild(0), + fLastChild(0), fPx(0), fPy(0), fPz(0), fEnergy(0), fMass(0), + fVx(0), fVy(0), fVz(0), fTime(0), fLifetime(0) {} + + TMCParticle(Int_t kS, Int_t kF, Int_t parent, + Int_t firstchild, Int_t lastchild, + Float_t px, Float_t py, Float_t pz, + Float_t energy, Float_t mass, + Float_t vx, Float_t vy, Float_t vz, + Float_t time, Float_t lifetime) : + + fKS(kS), + fKF(kF), + fParent(parent), + fFirstChild(firstchild), + fLastChild(lastchild), + fPx(px), + fPy(py), + fPz(pz), + fEnergy(energy), + fMass(mass), + fVx(vx), + fVy(vy), + fVz(vz), + fTime(time), + fLifetime(lifetime) { } + + + ~TMCParticle() override { } + + Int_t GetKS() const {return fKS;} + Int_t GetKF() const {return fKF;} + Int_t GetParent() const {return fParent;} + Int_t GetFirstChild() const {return fFirstChild;} + Int_t GetLastChild() const {return fLastChild;} + + Float_t GetPx() const {return fPx;} + Float_t GetPy() const {return fPy;} + Float_t GetPz() const {return fPz;} + Float_t GetEnergy() const {return fEnergy;} + Float_t GetMass() const {return fMass;} + + Float_t GetVx() const {return fVx;} + Float_t GetVy() const {return fVy;} + Float_t GetVz() const {return fVz;} + Float_t GetTime() const {return fTime;} + Float_t GetLifetime() const {return fLifetime;} + const char *GetName() const override; + + virtual void SetKS(Int_t kS) {fKS=kS;} + virtual void SetKF(Int_t kF) {fKF=kF;} + virtual void SetParent(Int_t parent) {fParent=parent;} + virtual void SetFirstChild(Int_t first) {fFirstChild=first;} + virtual void SetLastChild(Int_t last) {fLastChild=last;} + + virtual void SetPx(Float_t px) {fPx=px;} + virtual void SetPy(Float_t py) {fPy=py;} + virtual void SetPz(Float_t pz) {fPz=pz;} + virtual void SetEnergy(Float_t energy) {fEnergy=energy;} + virtual void SetMass(Float_t mass) {fMass=mass;} + + virtual void SetVx(Float_t vx) {fVx=vx;} + virtual void SetVy(Float_t vy) {fVy=vy;} + virtual void SetVz(Float_t vz) {fVz=vz;} + virtual void SetTime(Float_t time) {fTime=time;} + virtual void SetLifetime(Float_t lifetime) {fLifetime=lifetime;} + + + void ls(Option_t* option) const override; + + ClassDefOverride(TMCParticle,1) // LUJETS particles data record. +}; + +#endif diff --git a/GeneratorParam/TPythia6.cxx b/GeneratorParam/TPythia6.cxx new file mode 100644 index 0000000..b81a8ef --- /dev/null +++ b/GeneratorParam/TPythia6.cxx @@ -0,0 +1,695 @@ +// @(#)root/pythia6:$Id$ +// Author: Rene Brun 19/10/99 + +/** \class TPythia6 + \ingroup pythia6 + +TPythia is an interface class to F77 version of Pythia 6.2 + +To use this class you must install a version of pythia6. +See the installation instructions at + http://root.cern/root/Install.html + +CERNLIB event generators, written by T.Sjostrand. +For the details about these generators look at Pythia/Jetset manual: + +\verbatim +****************************************************************************** +** ** +** ** +** PPP Y Y TTTTT H H III A JJJJ EEEE TTTTT SSS EEEE TTTTT ** +** P P Y Y T H H I A A J E T S E T ** +** PPP Y T HHHHH I AAAAA J EEE T SSS EEE T ** +** P Y T H H I A A J J E T S E T ** +** P Y T H H III A A JJ EEEE T SSS EEEE T ** +** ** +** ** +** *......* Welcome to the Lund Monte Carlo! ** +** *:::!!:::::::::::* ** +** *::::::!!::::::::::::::* This is PYTHIA version 5.720 ** +** *::::::::!!::::::::::::::::* Last date of change: 29 Nov 1995 ** +** *:::::::::!!:::::::::::::::::* ** +** *:::::::::!!:::::::::::::::::* This is JETSET version 7.408 ** +** *::::::::!!::::::::::::::::*! Last date of change: 23 Aug 1995 ** +** *::::::!!::::::::::::::* !! ** +** !! *:::!!:::::::::::* !! Main author: ** +** !! !* -><- * !! Torbjorn Sjostrand ** +** !! !! !! Dept. of theoretical physics 2 ** +** !! !! !! University of Lund ** +** !! !! Solvegatan 14A ** +** !! ep !! S-223 62 Lund, Sweden ** +** !! !! phone: +46 - 46 - 222 48 16 ** +** !! pp !! E-mail: torbjorn@thep.lu.se ** +** !! e+e- !! ** +** !! !! Copyright Torbjorn Sjostrand ** +** !! and CERN, Geneva 1993 ** +** ** +** ** +** The latest program versions and documentation is found on WWW address ** +** http://thep.lu.se/tf2/staff/torbjorn/Welcome.html ** +** ** +** When you cite these programs, priority should always be given to the ** +** latest published description. Currently this is ** +** T. Sjostrand, Computer Physics Commun. 82 (1994) 74. ** +** The most recent long description (unpublished) is ** +** T. Sjostrand, LU TP 95-20 and CERN-TH.7112/93 (revised August 1995). ** +** Also remember that the programs, to a large extent, represent original ** +** physics research. Other publications of special relevance to your ** +** studies may therefore deserve separate mention. ** +** ** +** ** +****************************************************************************** +\endverbatim +*/ + +#include "TPythia6.h" + +#include "TClonesArray.h" +#include "TMCParticle.h" +#include "TParticle.h" +#include "snprintf.h" +#include "strlcpy.h" + +TPythia6* TPythia6::fgInstance = nullptr; + + +#ifndef WIN32 +# define pydiff pydiff_ +# define pyevnt pyevnt_ +# define pyinit pyinit_ +# define pychge pychge_ +# define pycomp pycomp_ +# define pyedit pyedit_ +# define pyexec pyexec_ +# define pyhepc pyhepc_ +# define pygive pygive_ +# define pylist pylist_ +# define pymass pymass_ +# define pyname pyname_ +# define pyr pyr_ +# define pyrget pyrget_ +# define pyrset pyrset_ +# define pystat pystat_ +# define pytest pytest_ +# define pytune pytune_ +# define pyupda pyupda_ +# define py1ent py1ent_ +# ifdef PYTHIA6_DOUBLE_UNDERSCORE +# define tpythia6_open_fortran_file tpythia6_open_fortran_file__ +# define tpythia6_close_fortran_file tpythia6_close_fortran_file__ +# define pythia6_common_address pythia6_common_address__ +# elif PYTHIA6_SINGLE_UNDERSCORE +# define tpythia6_open_fortran_file tpythia6_open_fortran_file_ +# define tpythia6_close_fortran_file tpythia6_close_fortran_file_ +# define pythia6_common_address pythia6_common_address +# else +# define pythia6_common_address pythia6_common_address +# define tpythia6_open_fortran_file tpythia6_open_fortran_file_ +# define tpythia6_close_fortran_file tpythia6_close_fortran_file_ +# endif +# define type_of_call +#else +# define pydiff PYDIFF +# define pyevnt PYEVNT +# define pyinit PYINIT +# define pychge PYCHGE +# define pycomp PYCOMP +# define pyedit PYEDIT +# define pyexec PYEXEC +# define pygive PYGIVE +# define pyhepc PYHEPC +# define pylist PYLIST +# define pymass PYMASS +# define pyname PYNAME +# define pyr PYR +# define pyrget PYRGET +# define pyrset PYRSET +# define pystat PYSTAT +# define pytest PYTEST +# define pytune PYTUNE +# define pyupda PYUPDA +# define py1ent PY1ENT +# define tpythia6_open_fortran_file TPYTHIA6_OPEN_FORTRAN_FILE +# define tpythia6_close_fortran_file TPYTHIA6_CLOSE_FORTRAN_FILE +# define type_of_call _stdcall +#endif + + +extern "C" void type_of_call pyevnt(); +extern "C" void type_of_call pystat(int *key); +extern "C" void type_of_call pylist(int *key); +extern "C" void type_of_call pyedit(int *medit); +extern "C" void type_of_call pydiff(); +extern "C" void type_of_call pyexec(); +extern "C" void type_of_call pygive(const char *param, Long_t lparam); +extern "C" void type_of_call pyhepc(int *mconv); +extern "C" void type_of_call pylist(int *flag); +extern "C" int type_of_call pychge(int *kf); +extern "C" int type_of_call pycomp(int *kf); +extern "C" double type_of_call pymass(int *kf); +extern "C" void type_of_call pyname(int *kf, char *name, Long_t l_name); +extern "C" int type_of_call pyr(int *dummy); +extern "C" int type_of_call pyrget(int *lun, int *move); +extern "C" int type_of_call pyrset(int *lun, int *move); +extern "C" int type_of_call pytest(int *flag); +extern "C" int type_of_call pytune(int *itune); +extern "C" int type_of_call pyupda(int *mupda, int *lun); +extern "C" void type_of_call py1ent(Int_t&, Int_t&, Double_t&, Double_t&, Double_t&); + +#ifndef WIN32 +extern "C" void type_of_call pyinit(char *frame, char *beam, char *target, + double *win, Long_t l_frame, Long_t l_beam, + Long_t l_target); +#else +extern "C" void type_of_call pyinit(char *frame, Long_t l_frame, + char *beam, Long_t l_beam, + char *target, Long_t l_target, + double *win + ); +#endif + +extern "C" { + void* pythia6_common_address(const char*); + void type_of_call tpythia6_open_fortran_file(int* lun, char* name, int); + void type_of_call tpythia6_close_fortran_file(int* lun); +} + +ClassImp(TPythia6); + +/** \class TPythia6::TPythia6Cleaner + \ingroup pythia6 + +Utility class to manage the TPythia6 instance +*/ + +TPythia6::TPythia6Cleaner::TPythia6Cleaner() { +} + +//////////////////////////////////////////////////////////////////////////////// +///delete the TPythia6 insntance + +TPythia6::TPythia6Cleaner::~TPythia6Cleaner() { + if (TPythia6::fgInstance) { + delete TPythia6::fgInstance; + TPythia6::fgInstance = 0; + } +} + +//------------------------------------------------------------------------------ +// constructor is not supposed to be called from the outside - only +// Initialize() method +//////////////////////////////////////////////////////////////////////////////// +/// TPythia6 constructor: creates a TClonesArray in which it will store all +/// particles. Note that there may be only one functional TPythia6 object +/// at a time, so it's not use to create more than one instance of it. + +TPythia6::TPythia6() : TGenerator("TPythia6","TPythia6") { + // Protect against multiple objects. All access should be via the + // Instance member function. + if (fgInstance) + Fatal("TPythia6", "There's already an instance of TPythia6"); + + delete fParticles; // was allocated as TObjArray in TGenerator + + fParticles = new TClonesArray("TMCParticle",50); + + // initialize common-blocks + // the functions/subroutines referenced by TPythia6 can be found + // at ftp://root.cern/root/pythia6.tar.gz + + fPyjets = (Pyjets_t*) pythia6_common_address("PYJETS"); + fPydat1 = (Pydat1_t*) pythia6_common_address("PYDAT1"); + fPydat2 = (Pydat2_t*) pythia6_common_address("PYDAT2"); + fPydat3 = (Pydat3_t*) pythia6_common_address("PYDAT3"); + fPydat4 = (Pydat4_t*) pythia6_common_address("PYDAT4"); + fPydatr = (Pydatr_t*) pythia6_common_address("PYDATR"); + fPysubs = (Pysubs_t*) pythia6_common_address("PYSUBS"); + fPypars = (Pypars_t*) pythia6_common_address("PYPARS"); + fPyint1 = (Pyint1_t*) pythia6_common_address("PYINT1"); + fPyint2 = (Pyint2_t*) pythia6_common_address("PYINT2"); + fPyint3 = (Pyint3_t*) pythia6_common_address("PYINT3"); + fPyint4 = (Pyint4_t*) pythia6_common_address("PYINT4"); + fPyint5 = (Pyint5_t*) pythia6_common_address("PYINT5"); + fPyint6 = (Pyint6_t*) pythia6_common_address("PYINT6"); + fPyint7 = (Pyint7_t*) pythia6_common_address("PYINT7"); + fPyint8 = (Pyint8_t*) pythia6_common_address("PYINT8"); + fPyint9 = (Pyint9_t*) pythia6_common_address("PYINT9"); + fPymssm = (Pymssm_t*) pythia6_common_address("PYMSSM"); + fPyssmt = (Pyssmt_t*) pythia6_common_address("PYSSMT"); + fPyints = (Pyints_t*) pythia6_common_address("PYINTS"); + fPybins = (Pybins_t*) pythia6_common_address("PYBINS"); +} + +//////////////////////////////////////////////////////////////////////////////// + +TPythia6::TPythia6(const TPythia6& p6) : + TGenerator(p6), + fPyjets(p6.fPyjets), + fPydat1(p6.fPydat1), + fPydat2(p6.fPydat2), + fPydat3(p6.fPydat3), + fPydat4(p6.fPydat4), + fPydatr(p6.fPydatr), + fPysubs(p6.fPysubs), + fPypars(p6.fPypars), + fPyint1(p6.fPyint1), + fPyint2(p6.fPyint2), + fPyint3(p6.fPyint3), + fPyint4(p6.fPyint4), + fPyint5(p6.fPyint5), + fPyint6(p6.fPyint6), + fPyint7(p6.fPyint7), + fPyint8(p6.fPyint8), + fPyint9(p6.fPyint9), + fPymssm(p6.fPymssm), + fPyssmt(p6.fPyssmt), + fPyints(p6.fPyints), + fPybins(p6.fPybins) +{ } + +//////////////////////////////////////////////////////////////////////////////// +/// Destroys the object, deletes and disposes all TMCParticles currently on list. + +TPythia6::~TPythia6() +{ + if (fParticles) { + fParticles->Delete(); + delete fParticles; + fParticles = 0; + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// model of automatic memory cleanup suggested by Jim Kowalkovski: +/// destructor for local static variable `cleaner' is always called in the end +/// of the job thus deleting the only TPythia6 instance + +TPythia6* TPythia6::Instance() { + static TPythia6::TPythia6Cleaner cleaner; + return fgInstance ? fgInstance : (fgInstance=new TPythia6()) ; +} + + + + + +//////////////////////////////////////////////////////////////////////////////// +/// generate event and copy the information from /HEPEVT/ to fPrimaries + +void TPythia6::GenerateEvent() { + pyevnt(); + ImportParticles(); +} + +//////////////////////////////////////////////////////////////////////////////// +///interface with fortran i/o + +void TPythia6::OpenFortranFile(int lun, char* name) { + tpythia6_open_fortran_file(&lun, name, strlen(name)); +} + +//////////////////////////////////////////////////////////////////////////////// +///interface with fortran i/o + +void TPythia6::CloseFortranFile(int lun) { + tpythia6_close_fortran_file(&lun); +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Fills TObjArray fParticles list with particles from common LUJETS +/// Old contents of a list are cleared. This function should be called after +/// any change in common LUJETS, however GetParticles() method calls it +/// automatically - user don't need to care about it. In case you make a call +/// to LuExec() you must call this method yourself to transfer new data from +/// common LUJETS to the fParticles list. + +TObjArray *TPythia6::ImportParticles(Option_t *) +{ + fParticles->Clear(); + Int_t numpart = fPyjets->N; + TClonesArray &a = *((TClonesArray*)fParticles); + for (Int_t i = 0; iK[0][i] , + fPyjets->K[1][i] , + fPyjets->K[2][i] , + fPyjets->K[3][i] , + fPyjets->K[4][i] , + fPyjets->P[0][i] , + fPyjets->P[1][i] , + fPyjets->P[2][i] , + fPyjets->P[3][i] , + fPyjets->P[4][i] , + fPyjets->V[0][i] , + fPyjets->V[1][i] , + fPyjets->V[2][i] , + fPyjets->V[3][i] , + fPyjets->V[4][i]); + } + return fParticles; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Default primary creation method. It reads the /HEPEVT/ common block which +/// has been filled by the GenerateEvent method. If the event generator does +/// not use the HEPEVT common block, This routine has to be overloaded by +/// the subclasses. +/// The function loops on the generated particles and store them in +/// the TClonesArray pointed by the argument particles. +/// The default action is to store only the stable particles (ISTHEP = 1) +/// This can be demanded explicitly by setting the option = "Final" +/// If the option = "All", all the particles are stored. +/// + +Int_t TPythia6::ImportParticles(TClonesArray *particles, Option_t *option) +{ + if (particles == 0) return 0; + TClonesArray &clonesParticles = *particles; + clonesParticles.Clear(); + Int_t numpart = fPyjets->N; + Int_t nparts=0; + if (!strcmp(option,"") || !strcmp(option,"Final")) { + for (Int_t i = 0; iK[0][i] == 1) { + // + // Use the common block values for the TParticle constructor + // + new(clonesParticles[nparts]) TParticle( + fPyjets->K[1][i] , + fPyjets->K[0][i] , + fPyjets->K[2][i] , + -1, + fPyjets->K[3][i] , + fPyjets->K[4][i] , + fPyjets->P[0][i] , + fPyjets->P[1][i] , + fPyjets->P[2][i] , + fPyjets->P[3][i] , + fPyjets->V[0][i] , + fPyjets->V[1][i] , + fPyjets->V[2][i] , + fPyjets->V[3][i]); + + // if(gDebug) printf("%d %d %d! ",i,fPyjets->K[1][i],numpart); + nparts++; + } + } + } else if (!strcmp(option,"All")) { + for (Int_t i = 0; iK[1][i] , + fPyjets->K[0][i] , + fPyjets->K[2][i] , + -1, + fPyjets->K[3][i] , + fPyjets->K[4][i] , + fPyjets->P[0][i] , + fPyjets->P[1][i] , + fPyjets->P[2][i] , + fPyjets->P[3][i] , + fPyjets->V[0][i] , + fPyjets->V[1][i] , + fPyjets->V[2][i] , + fPyjets->V[3][i]); + } + nparts=numpart; + } + + return nparts; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Calls PyInit with the same parameters after performing some checking, +/// sets correct title. This method should preferably be called instead of PyInit. +/// PURPOSE: to initialize the generation procedure. +/// ARGUMENTS: See documentation for details. +/// - frame: - specifies the frame of the experiment: +/// "CMS","FIXT","USER","FOUR","FIVE","NONE" +/// - beam, +/// - target: - beam and target particles (with additionaly charges, tildes or "bar": +/// e,nu_e,mu,nu_mu,tau,nu_tau,gamma,pi,n,p,Lambda,Sigma,Xi,Omega, +/// pomeron,reggeon +/// - win: - related to energy system: +/// - for frame=="CMS" - total energy of system +/// - for frame=="FIXT" - momentum of beam particle +/// - for frame=="USER" - dummy - see documentation. + +void TPythia6::Initialize(const char *frame, const char *beam, const char *target, float win) +{ + char cframe[4]; + strlcpy(cframe,frame,4); + char cbeam[10]; + strlcpy(cbeam,beam,10); + char ctarget[10]; + strlcpy(ctarget,target,10); + + // For frames "3MOM", "4MOM" and "5MOM" see p. 181-182 of the version 6 manual, + // http://home.thep.lu.se/~torbjorn/pythia/lutp0613man2.pdf + // their usage may depend on the version of Pythia6 used + if ( (!strncmp(frame, "CMS" ,3)) && + (!strncmp(frame, "FIXT" ,4)) && + (!strncmp(frame, "USER" ,4)) && + (!strncmp(frame, "FOUR" ,4)) && + (!strncmp(frame, "FIVE" ,4)) && + (!strncmp(frame, "3MOM" ,4)) && + (!strncmp(frame, "4MOM" ,4)) && + (!strncmp(frame, "5MOM" ,4)) && + (!strncmp(frame, "NONE" ,4)) ) { + printf("WARNING! In TPythia6:Initialize():\n"); + printf(" specified frame=%s is neither of CMS,FIXT,USER,FOUR,FIVE,NONE,3MOM,4MOM,5MOM\n",frame); + printf(" resetting to \"CMS\" ."); + snprintf(cframe,4,"CMS"); + } + + if ( (!strncmp(beam, "e" ,1)) && + (!strncmp(beam, "nu_e" ,4)) && + (!strncmp(beam, "mu" ,2)) && + (!strncmp(beam, "nu_mu" ,5)) && + (!strncmp(beam, "tau" ,3)) && + (!strncmp(beam, "nu_tau" ,6)) && + (!strncmp(beam, "gamma" ,5)) && + (!strncmp(beam, "pi" ,2)) && + (!strncmp(beam, "n" ,1)) && + (!strncmp(beam, "p" ,1)) && + (!strncmp(beam, "Lambda" ,6)) && + (!strncmp(beam, "Sigma" ,5)) && + (!strncmp(beam, "Xi" ,2)) && + (!strncmp(beam, "Omega" ,5)) && + (!strncmp(beam, "pomeron" ,7)) && + (!strncmp(beam, "reggeon" ,7)) ) { + printf("WARNING! In TPythia6:Initialize():\n"); + printf(" specified beam=%s is unrecognized .\n",beam); + printf(" resetting to \"p+\" ."); + snprintf(cbeam,8,"p+"); + } + + if ( (!strncmp(target, "e" ,1)) && + (!strncmp(target, "nu_e" ,4)) && + (!strncmp(target, "mu" ,2)) && + (!strncmp(target, "nu_mu" ,5)) && + (!strncmp(target, "tau" ,3)) && + (!strncmp(target, "nu_tau" ,6)) && + (!strncmp(target, "gamma" ,5)) && + (!strncmp(target, "pi" ,2)) && + (!strncmp(target, "n" ,1)) && + (!strncmp(target, "p" ,1)) && + (!strncmp(target, "Lambda" ,6)) && + (!strncmp(target, "Sigma" ,5)) && + (!strncmp(target, "Xi" ,2)) && + (!strncmp(target, "Omega" ,5)) && + (!strncmp(target, "pomeron" ,7)) && + (!strncmp(target, "reggeon" ,7)) ){ + printf("WARNING! In TPythia6:Initialize():\n"); + printf(" specified target=%s is unrecognized.\n",target); + printf(" resetting to \"p+\" ."); + snprintf(ctarget,8,"p+"); + } + + Pyinit(cframe, cbeam ,ctarget, win); + + char atitle[64]; + snprintf(atitle, sizeof(atitle)," %s-%s at %g GeV", cbeam, ctarget, win); + SetTitle(atitle); +} + + +void TPythia6::Pyinit(char* frame, char* beam, char* target, double win) +{ + // Calls Pyinit with the same parameters after performing some checking, + // sets correct title. This method should preferably be called instead of PyInit. + // PURPOSE: to initialize the generation procedure. + // ARGUMENTS: See documentation for details. + // frame: - specifies the frame of the experiment: + // "CMS","FIXT","USER","FOUR","FIVE","NONE" + // beam, + // target: - beam and target particles (with additionaly charges, + // tildes or "bar": + // e,nu_e,mu,nu_mu,tau,nu_tau,gamma,pi,n,p,Lambda,Sigma,Xi,Omega, + // pomeron,reggeon + // win: - related to energy system: + // for frame=="CMS" - total energy of system + // for frame=="FIXT" - momentum of beam particle + // for frame=="USER" - dummy - see documentation. + + Double_t lwin = win; + Long_t s1 = strlen(frame); + Long_t s2 = strlen(beam); + Long_t s3 = strlen(target); +#ifndef WIN32 + pyinit(frame,beam,target,&lwin,s1,s2,s3); +#else + pyinit(frame, s1, beam , s2, target, s3, &lwin); +#endif +} + + +int TPythia6::Pycomp(int kf) { + //interface with fortran routine pycomp + return pycomp(&kf); +} + +void TPythia6::Pyedit(int medit) { + //interface with fortran routine pyedit + pyedit(&medit); + ImportParticles(); +} + +void TPythia6::Pydiff() { + //interface with fortran routine pydiff + pydiff(); +} + +void TPythia6::Pyevnt() { + //interface with fortran routine pyevnt + pyevnt(); +} + +void TPythia6::Pyexec() { + //interface with fortran routine pyexec + pyexec(); +} + +void TPythia6::Pygive(const char *param) { + //interface with fortran routine pygive + Long_t lparam = strlen(param); + pygive(param,lparam); +} + +void TPythia6::Pyhepc(int mconv) { + //interface with fortran routine pyhepc + pyhepc(&mconv); +} + +void TPythia6::Pylist(int flag) { + //interface with fortran routine pylist + pylist(&flag); +} + +void TPythia6::Pyname(int kf, char* name) { + //Note that the array name must be dimensioned in the calling program + //to at least name[16] + + pyname(&kf,name,15); + // cut trailing blanks to get C string + name[15] = 0; + //for (int i=15; (i>=0) && (name[i] == ' '); i--) { + // name[i] = 0; + // } +} + +double TPythia6::Pyr(int idummy) { + //interface with fortran routine pyr + return pyr(&idummy); +} + +void TPythia6::Pyrget(int lun, int move) { + //interface with fortran routine pyrget + pyrget(&lun,&move); +} + +void TPythia6::Pyrset(int lun, int move) { + //interface with fortran routine pyrset + pyrset(&lun,&move); +} + +void TPythia6::Pystat(int flag) { + //interface with fortran routine pystat + pystat(&flag); +} + +void TPythia6::Pytest(int flag) { + //interface with fortran routine pytest + pytest(&flag); +} + +void TPythia6::Pytune(int itune) { + //interface with fortran routine pytune + pytune(&itune); +} + +void TPythia6::Pyupda(int mupda, int lun) { + //interface with fortran routine pyupda + pyupda(&mupda,&lun); +} + +double TPythia6::Pymass(int kf) { + //interface with fortran routine pymass + return pymass(&kf); +} + +int TPythia6::Pychge(int kf) { + //interface with fortran routine pychge + return pychge(&kf); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Add one entry to the event record, i.e. either a parton or a +/// particle. +/// +/// - IP: normally line number for the parton/particle. There are two +/// exceptions: +/// - If IP = 0: line number 1 is used and PYEXEC is called. +/// - If IP < 0: line -IP is used, with status code K(-IP,2)=2 +/// rather than 1; thus a parton system may be built +/// up by filling all but the last parton of the +/// system with IP < 0. +/// - KF: parton/particle flavour code (PDG code) +/// - PE: parton/particle energy. If PE is smaller than the mass, +/// the parton/particle is taken to be at rest. +/// - THETA: +/// - PHI: polar and azimuthal angle for the momentum vector of the +/// parton/particle. + +void TPythia6::Py1ent(Int_t ip, Int_t kf, Double_t pe, Double_t theta, Double_t phi) +{ + py1ent(ip, kf, pe, theta, phi); +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Exemplary setup of Pythia parameters: +/// Switches on processes 102,123,124 (Higgs generation) and switches off +/// interactions, fragmentation, ISR, FSR... + +void TPythia6::SetupTest() +{ + SetMSEL(0); // full user controll; + + SetMSUB(102,1); // g + g -> H0 + SetMSUB(123,1); // f + f' -> f + f' + H0 + SetMSUB(124,1); // f + f' -> f" + f"' + H0 + + + SetPMAS(6,1,175.0); // mass of TOP + SetPMAS(25,1,300); // mass of Higgs + + + SetCKIN(1,290.0); // range of allowed mass + SetCKIN(2,310.0); + + SetMSTP(61, 0); // switch off ISR + SetMSTP(71, 0); // switch off FSR + SetMSTP(81, 0); // switch off multiple interactions + SetMSTP(111, 0); // switch off fragmentation/decay +} diff --git a/GeneratorParam/TPythia6Decayer.cxx b/GeneratorParam/TPythia6Decayer.cxx new file mode 100644 index 0000000..db533e8 --- /dev/null +++ b/GeneratorParam/TPythia6Decayer.cxx @@ -0,0 +1,630 @@ +// @(#)root/pythia6:$Id$ +// Author: Christian Holm Christensen 22/04/06 +// Much of this code has been lifted from AliROOT. + +/************************************************************************* + * Copyright (C) 1995-2006, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +/** \class TPythia6Decayer + \ingroup pythia6 + +This class implements the TVirtualMCDecayer interface. + +The TPythia6 singleton object is used to decay particles. Note, that since +this class modifies common blocks (global variables) it is defined as a +singleton. +*/ + +#include "TPythia6.h" +#include "TPythia6Decayer.h" +#include "TPDGCode.h" +#include "TLorentzVector.h" +#include "TClonesArray.h" + +ClassImp(TPythia6Decayer); + +TPythia6Decayer* TPythia6Decayer::fgInstance = 0; + +//////////////////////////////////////////////////////////////////////////////// +/// Get the singleton object. + +TPythia6Decayer* TPythia6Decayer::Instance() +{ + if (!fgInstance) fgInstance = new TPythia6Decayer; + return fgInstance; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Constructor + +TPythia6Decayer::TPythia6Decayer() + : fDecay(kMaxDecay), + fBraPart(501) +{ + fBraPart.Reset(1); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Initialize the decayer + +void TPythia6Decayer::Init() +{ + static Bool_t init = kFALSE; + if (init) return; + init = kTRUE; + ForceDecay(); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Decay a particle of type IDPART (PDG code) and momentum P. + +void TPythia6Decayer::Decay(Int_t idpart, TLorentzVector* p) +{ + if (!p) return; + TPythia6::Instance()->Py1ent(0, idpart, p->Energy(), p->Theta(), p->Phi()); + TPythia6::Instance()->GetPrimaries(); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Get the decay products into the passed PARTICLES TClonesArray of +/// TParticles + +Int_t TPythia6Decayer::ImportParticles(TClonesArray *particles) +{ + return TPythia6::Instance()->ImportParticles(particles,"All"); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Force a particular decay type + +void TPythia6Decayer::SetForceDecay(Int_t type) +{ + if (type > kMaxDecay) { + Warning("SetForceDecay", "Invalid decay mode: %d", type); + return; + } + fDecay = EDecayType(type); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Force a particle decay mode + +void TPythia6Decayer::ForceDecay() +{ + EDecayType decay=fDecay; + TPythia6::Instance()->SetMSTJ(21,2); + if (decay == kNoDecayHeavy) return; + + // + // select mode + Int_t products[3]; + Int_t mult[3]; + + switch (decay) { + case kHardMuons: + products[0] = 13; + products[1] = 443; + products[2] = 100443; + mult[0] = 1; + mult[1] = 1; + mult[2] = 1; + ForceParticleDecay( 511, products, mult, 3); + ForceParticleDecay( 521, products, mult, 3); + ForceParticleDecay( 531, products, mult, 3); + ForceParticleDecay( 5122, products, mult, 3); + ForceParticleDecay( 5132, products, mult, 3); + ForceParticleDecay( 5232, products, mult, 3); + ForceParticleDecay( 5332, products, mult, 3); + ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X + ForceParticleDecay( 443, 13, 2); // J/Psi -> mu+ mu- + + ForceParticleDecay( 411,13,1); // D+/- + ForceParticleDecay( 421,13,1); // D0 + ForceParticleDecay( 431,13,1); // D_s + ForceParticleDecay( 4122,13,1); // Lambda_c + ForceParticleDecay( 4132,13,1); // Xsi_c + ForceParticleDecay( 4232,13,1); // Sigma_c + ForceParticleDecay( 4332,13,1); // Omega_c + break; + case kSemiMuonic: + ForceParticleDecay( 411,13,1); // D+/- + ForceParticleDecay( 421,13,1); // D0 + ForceParticleDecay( 431,13,1); // D_s + ForceParticleDecay( 4122,13,1); // Lambda_c + ForceParticleDecay( 4132,13,1); // Xsi_c + ForceParticleDecay( 4232,13,1); // Sigma_c + ForceParticleDecay( 4332,13,1); // Omega_c + ForceParticleDecay( 511,13,1); // B0 + ForceParticleDecay( 521,13,1); // B+/- + ForceParticleDecay( 531,13,1); // B_s + ForceParticleDecay( 5122,13,1); // Lambda_b + ForceParticleDecay( 5132,13,1); // Xsi_b + ForceParticleDecay( 5232,13,1); // Sigma_b + ForceParticleDecay( 5332,13,1); // Omega_b + break; + case kDiMuon: + ForceParticleDecay( 113,13,2); // rho + ForceParticleDecay( 221,13,2); // eta + ForceParticleDecay( 223,13,2); // omega + ForceParticleDecay( 333,13,2); // phi + ForceParticleDecay( 443,13,2); // J/Psi + ForceParticleDecay(100443,13,2);// Psi' + ForceParticleDecay( 553,13,2); // Upsilon + ForceParticleDecay(100553,13,2);// Upsilon' + ForceParticleDecay(200553,13,2);// Upsilon'' + break; + case kSemiElectronic: + ForceParticleDecay( 411,11,1); // D+/- + ForceParticleDecay( 421,11,1); // D0 + ForceParticleDecay( 431,11,1); // D_s + ForceParticleDecay( 4122,11,1); // Lambda_c + ForceParticleDecay( 4132,11,1); // Xsi_c + ForceParticleDecay( 4232,11,1); // Sigma_c + ForceParticleDecay( 4332,11,1); // Omega_c + ForceParticleDecay( 511,11,1); // B0 + ForceParticleDecay( 521,11,1); // B+/- + ForceParticleDecay( 531,11,1); // B_s + ForceParticleDecay( 5122,11,1); // Lambda_b + ForceParticleDecay( 5132,11,1); // Xsi_b + ForceParticleDecay( 5232,11,1); // Sigma_b + ForceParticleDecay( 5332,11,1); // Omega_b + break; + case kDiElectron: + ForceParticleDecay( 113,11,2); // rho + ForceParticleDecay( 333,11,2); // phi + ForceParticleDecay( 221,11,2); // eta + ForceParticleDecay( 223,11,2); // omega + ForceParticleDecay( 443,11,2); // J/Psi + ForceParticleDecay(100443,11,2);// Psi' + ForceParticleDecay( 553,11,2); // Upsilon + ForceParticleDecay(100553,11,2);// Upsilon' + ForceParticleDecay(200553,11,2);// Upsilon'' + break; + case kBJpsiDiMuon: + + products[0] = 443; + products[1] = 100443; + mult[0] = 1; + mult[1] = 1; + + ForceParticleDecay( 511, products, mult, 2); // B0 -> J/Psi (Psi') X + ForceParticleDecay( 521, products, mult, 2); // B+/- -> J/Psi (Psi') X + ForceParticleDecay( 531, products, mult, 2); // B_s -> J/Psi (Psi') X + ForceParticleDecay( 5122, products, mult, 2); // Lambda_b -> J/Psi (Psi') X + ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X + ForceParticleDecay( 443,13,2); // J/Psi -> mu+ mu- + break; + case kBPsiPrimeDiMuon: + ForceParticleDecay( 511,100443,1); // B0 + ForceParticleDecay( 521,100443,1); // B+/- + ForceParticleDecay( 531,100443,1); // B_s + ForceParticleDecay( 5122,100443,1); // Lambda_b + ForceParticleDecay(100443,13,2); // Psi' + break; + case kBJpsiDiElectron: + ForceParticleDecay( 511,443,1); // B0 + ForceParticleDecay( 521,443,1); // B+/- + ForceParticleDecay( 531,443,1); // B_s + ForceParticleDecay( 5122,443,1); // Lambda_b + ForceParticleDecay( 443,11,2); // J/Psi + break; + case kBJpsi: + ForceParticleDecay( 511,443,1); // B0 + ForceParticleDecay( 521,443,1); // B+/- + ForceParticleDecay( 531,443,1); // B_s + ForceParticleDecay( 5122,443,1); // Lambda_b + break; + case kBPsiPrimeDiElectron: + ForceParticleDecay( 511,100443,1); // B0 + ForceParticleDecay( 521,100443,1); // B+/- + ForceParticleDecay( 531,100443,1); // B_s + ForceParticleDecay( 5122,100443,1); // Lambda_b + ForceParticleDecay(100443,11,2); // Psi' + break; + case kPiToMu: + ForceParticleDecay(211,13,1); // pi->mu + break; + case kKaToMu: + ForceParticleDecay(321,13,1); // K->mu + break; + case kWToMuon: + ForceParticleDecay( 24, 13,1); // W -> mu + break; + case kWToCharm: + ForceParticleDecay( 24, 4,1); // W -> c + break; + case kWToCharmToMuon: + ForceParticleDecay( 24, 4,1); // W -> c + ForceParticleDecay( 411,13,1); // D+/- -> mu + ForceParticleDecay( 421,13,1); // D0 -> mu + ForceParticleDecay( 431,13,1); // D_s -> mu + ForceParticleDecay( 4122,13,1); // Lambda_c + ForceParticleDecay( 4132,13,1); // Xsi_c + ForceParticleDecay( 4232,13,1); // Sigma_c + ForceParticleDecay( 4332,13,1); // Omega_c + break; + case kZDiMuon: + ForceParticleDecay( 23, 13,2); // Z -> mu+ mu- + break; + case kHadronicD: + ForceHadronicD(); + break; + case kPhiKK: + ForceParticleDecay(333,321,2); // Phi->K+K- + break; + case kOmega: + ForceOmega(); + case kAll: + break; + case kNoDecay: + TPythia6::Instance()->SetMSTJ(21,0); + break; + case kNoDecayHeavy: break; // cannot get here: early return above + case kMaxDecay: break; + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Get the partial branching ratio for a particle of type IPART (a +/// PDG code). + +Float_t TPythia6Decayer::GetPartialBranchingRatio(Int_t ipart) +{ + Int_t kc = TPythia6::Instance()->Pycomp(TMath::Abs(ipart)); + // return TPythia6::Instance()->GetBRAT(kc); + return fBraPart[kc]; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Get the life-time of a particle of type KF (a PDG code). + +Float_t TPythia6Decayer::GetLifetime(Int_t kf) +{ + Int_t kc=TPythia6::Instance()->Pycomp(TMath::Abs(kf)); + return TPythia6::Instance()->GetPMAS(kc,4) * 3.3333e-12; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Read in particle data from an ASCII file. The file name must +/// previously have been set using the member function +/// SetDecayTableFile. + +void TPythia6Decayer::ReadDecayTable() +{ + if (fDecayTableFile.IsNull()) { + Warning("ReadDecayTable", "No file set"); + return; + } + Int_t lun = 15; + TPythia6::Instance()->OpenFortranFile(lun, + const_cast(fDecayTableFile.Data())); + TPythia6::Instance()->Pyupda(3,lun); + TPythia6::Instance()->CloseFortranFile(lun); +} + +// =================================================================== +// BEGIN COMMENT +// +// It would be better if the particle and decay information could be +// read from the current TDatabasePDG instance. +// +// However, it seems to me that some information is missing. In +// particular +// +// - The broadning cut-off, +// - Resonance width +// - Color charge +// - MWID (?) +// +// Further more, it's not clear to me at least, what all the +// parameters Pythia needs are. +// +// Code like the below could be used to make a temporary file that +// Pythia could then read in. Ofcourse, one could also manipulate +// the data structures directly, but that's propably more dangerous. +// +#if 0 +void PrintPDG(TParticlePDG* pdg) +{ + TParticlePDG* anti = pdg->AntiParticle(); + const char* antiName = (anti ? anti->GetName() : ""); + Int_t color = 0; + switch (TMath::Abs(pdg->PdgCode())) { + case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8: // Quarks + color = 1; break; + case 21: // Gluon + color = 2; break; + case 1103: + case 2101: case 2103: case 2203: + case 3101: case 3103: case 3201: case 3203: case 3303: + case 4101: case 4103: case 4201: case 4203: case 4301: case 4303: case 4403: + case 5101: case 5103: case 5201: case 5203: case 5301: case 5303: case 5401: + case 5403: case 5503: + // Quark combinations + color = -1; break; + case 1000001: case 1000002: case 1000003: case 1000004: case 1000005: + case 1000006: // super symmetric partners to quars + color = 1; break; + case 1000021: // ~g + color = 2; break; + case 2000001: case 2000002: case 2000003: case 2000004: case 2000005: + case 2000006: // R hadrons + color = 1; break; + case 3000331: case 3100021: case 3200111: case 3100113: case 3200113: + case 3300113: case 3400113: + // Technicolor + color = 2; break; + case 4000001: case 4000002: + color = 1; break; + case 9900443: case 9900441: case 9910441: case 9900553: case 9900551: + case 9910551: + color = 2; break; + } + std::cout << std::right + << " " << std::setw(9) << pdg->PdgCode() + << " " << std::left << std::setw(16) << pdg->GetName() + << " " << std::setw(16) << antiName + << std::right + << std::setw(3) << Int_t(pdg->Charge()) + << std::setw(3) << color + << std::setw(3) << (anti ? 1 : 0) + << std::fixed << std::setprecision(5) + << std::setw(12) << pdg->Mass() + << std::setw(12) << pdg->Width() + << std::setw(12) << 0 // Broad + << std::scientific + << " " << std::setw(13) << pdg->Lifetime() + << std::setw(3) << 0 // MWID + << std::setw(3) << pdg->Stable() + << std::endl; +} + +void MakeDecayList() +{ + TDatabasePDG* pdgDB = TDatabasePDG::Instance(); + pdgDB->ReadPDGTable(); + const THashList* pdgs = pdgDB->ParticleList(); + TParticlePDG* pdg = 0; + TIter nextPDG(pdgs); + while ((pdg = static_cast(nextPDG()))) { + // std::cout << "Processing " << pdg->GetName() << std::endl; + PrintPDG(pdg); + + TObjArray* decays = pdg->DecayList(); + TDecayChannel* decay = 0; + TIter nextDecay(decays); + while ((decay = static_cast(nextDecay()))) { + // std::cout << "Processing decay number " << decay->Number() << std::endl; + } + } +} +#endif +// END COMMENT +// =================================================================== + +//////////////////////////////////////////////////////////////////////////////// +/// write particle data to an ASCII file. The file name must +/// previously have been set using the member function +/// SetDecayTableFile. +/// +/// Users can use this function to make an initial decay list file, +/// which then can be edited by hand, and re-loaded into the decayer +/// using ReadDecayTable. +/// +/// The file syntax is +/// +/// particle_list : partcle_data +/// | particle_list particle_data +/// ; +/// particle_data : particle_info +/// | particle_info '\n' decay_list +/// ; +/// particle_info : See below +/// ; +/// decay_list : decay_entry +/// | decay_list decay_entry +/// ; +/// decay_entry : See below +/// +/// The particle_info consists of 13 fields: +/// +/// PDG code int +/// Name string +/// Anti-particle name string if there's no anti-particle, +/// then this field must be the +/// empty string +/// Electic charge int in units of |e|/3 +/// Color charge int in units of quark color charges +/// Have anti-particle int 1 of there's an anti-particle +/// to this particle, or 0 +/// otherwise +/// Mass float in units of GeV +/// Resonance width float +/// Max broadning float +/// Lifetime float +/// MWID int ??? (some sort of flag) +/// Decay int 1 if it decays. 0 otherwise +/// +/// The format to write these entries in are +/// +/// " %9 %-16s %-16s%3d%3d%3d%12.5f%12.5f%12.5f%13.gf%3d%d\n" +/// +/// The decay_entry consists of 8 fields: +/// +/// On/Off int 1 for on, -1 for off +/// Matrix element type int +/// Branching ratio float +/// Product 1 int PDG code of decay product 1 +/// Product 2 int PDG code of decay product 2 +/// Product 3 int PDG code of decay product 3 +/// Product 4 int PDG code of decay product 4 +/// Product 5 int PDG code of decay product 5 +/// +/// The format for these lines are +/// +/// " %5d%5d%12.5f%10d%10d%10d%10d%10d\n" +/// + +void TPythia6Decayer::WriteDecayTable() +{ + if (fDecayTableFile.IsNull()) { + Warning("ReadDecayTable", "No file set"); + return; + } + Int_t lun = 15; + TPythia6::Instance()->OpenFortranFile(lun, + const_cast(fDecayTableFile.Data())); + TPythia6::Instance()->Pyupda(1,lun); + TPythia6::Instance()->CloseFortranFile(lun); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Count number of decay products + +Int_t TPythia6Decayer::CountProducts(Int_t channel, Int_t particle) +{ + Int_t np = 0; + for (Int_t i = 1; i <= 5; i++) + if (TMath::Abs(TPythia6::Instance()->GetKFDP(channel,i)) == particle) np++; + return np; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Force golden D decay modes + +void TPythia6Decayer::ForceHadronicD() +{ + const Int_t kNHadrons = 4; + Int_t channel; + Int_t hadron[kNHadrons] = {411, 421, 431, 4112}; + + // for D+ -> K0* (-> K- pi+) pi+ + Int_t iKstar0 = 313; + Int_t iKstarbar0 = -313; + Int_t products[2] = {kKPlus, kPiMinus}, mult[2] = {1, 1}; + ForceParticleDecay(iKstar0, products, mult, 2); + + // for Ds -> Phi pi+ + Int_t iPhi = 333; + ForceParticleDecay(iPhi,kKPlus,2); // Phi->K+K- + + Int_t decayP1[kNHadrons][3] = { + {kKMinus, kPiPlus, kPiPlus}, + {kKMinus, kPiPlus, 0 }, + {kKPlus , iKstarbar0, 0 }, + {-1 , -1 , -1 } + }; + Int_t decayP2[kNHadrons][3] = { + {iKstarbar0, kPiPlus, 0 }, + {-1 , -1 , -1 }, + {iPhi , kPiPlus, 0 }, + {-1 , -1 , -1 } + }; + + TPythia6* pyth = TPythia6::Instance(); + for (Int_t ihadron = 0; ihadron < kNHadrons; ihadron++) { + Int_t kc = pyth->Pycomp(hadron[ihadron]); + pyth->SetMDCY(kc,1,1); + Int_t ifirst = pyth->GetMDCY(kc,2); + Int_t ilast = ifirst + pyth->GetMDCY(kc,3)-1; + + for (channel = ifirst; channel <= ilast; channel++) { + if ((pyth->GetKFDP(channel,1) == decayP1[ihadron][0] && + pyth->GetKFDP(channel,2) == decayP1[ihadron][1] && + pyth->GetKFDP(channel,3) == decayP1[ihadron][2] && + pyth->GetKFDP(channel,4) == 0) || + (pyth->GetKFDP(channel,1) == decayP2[ihadron][0] && + pyth->GetKFDP(channel,2) == decayP2[ihadron][1] && + pyth->GetKFDP(channel,3) == decayP2[ihadron][2] && + pyth->GetKFDP(channel,4) == 0)) { + pyth->SetMDME(channel,1,1); + } else { + pyth->SetMDME(channel,1,0); + fBraPart[kc] -= pyth->GetBRAT(channel); + } // selected channel ? + } // decay channels + } // hadrons +} + +//////////////////////////////////////////////////////////////////////////////// +/// +/// Force decay of particle into products with multiplicity mult + +void TPythia6Decayer::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult) +{ + TPythia6* pyth = TPythia6::Instance(); + + Int_t kc = pyth->Pycomp(particle); + pyth->SetMDCY(kc,1,1); + + Int_t ifirst = pyth->GetMDCY(kc,2); + Int_t ilast = ifirst + pyth->GetMDCY(kc,3)-1; + fBraPart[kc] = 1; + + // + // Loop over decay channels + for (Int_t channel= ifirst; channel <= ilast; channel++) { + if (CountProducts(channel,product) >= mult) { + pyth->SetMDME(channel,1,1); + } else { + pyth->SetMDME(channel,1,0); + fBraPart[kc]-=pyth->GetBRAT(channel); + } + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// +/// Force decay of particle into products with multiplicity mult + +void TPythia6Decayer::ForceParticleDecay(Int_t particle, Int_t* products, + Int_t* mult, Int_t npart) +{ + TPythia6* pyth = TPythia6::Instance(); + + Int_t kc = pyth->Pycomp(particle); + pyth->SetMDCY(kc,1,1); + Int_t ifirst = pyth->GetMDCY(kc,2); + Int_t ilast = ifirst+pyth->GetMDCY(kc,3)-1; + fBraPart[kc] = 1; + // + // Loop over decay channels + for (Int_t channel = ifirst; channel <= ilast; channel++) { + Int_t nprod = 0; + for (Int_t i = 0; i < npart; i++) + nprod += (CountProducts(channel, products[i]) >= mult[i]); + if (nprod) + pyth->SetMDME(channel,1,1); + else { + pyth->SetMDME(channel,1,0); + fBraPart[kc] -= pyth->GetBRAT(channel); + } + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Force Omega -> Lambda K- Decay + +void TPythia6Decayer::ForceOmega() +{ + TPythia6* pyth = TPythia6::Instance(); + + Int_t kc = pyth->Pycomp(3334); + pyth->SetMDCY(kc,1,1); + Int_t ifirst = pyth->GetMDCY(kc,2); + Int_t ilast = ifirst + pyth->GetMDCY(kc,3)-1; + for (Int_t channel = ifirst; channel <= ilast; channel++) { + if (pyth->GetKFDP(channel,1) == kLambda0 && + pyth->GetKFDP(channel,2) == kKMinus && + pyth->GetKFDP(channel,3) == 0) + pyth->SetMDME(channel,1,1); + else + pyth->SetMDME(channel,1,0); + // selected channel ? + } // decay channels +}