//////////////////////////////////////////////////////////////////
//
// JSFHadronizer
//
// Hadronize parton information in JSFParton class and
// create JSFGeneratorParticle class.
//
//$Id: JSFHadronizer.cxx,v 1.14 2003/02/03 00:09:25 miyamoto Exp $
//////////////////////////////////////////////////////////////////
#include "TRandom.h"
#include "JSFSteer.h"
#include "JSFLCFULL.h"
#include "JSFHadronizer.h"
#if __PYTHIA_VERSION__ >= 6
extern "C" {
extern void tauint_(int *inut, int *iout, int *jak1,
int *jak2, int *itdkrc, int *keya1, float *xk0dec);
extern void lutuhl_(int *one, int *ihlon);
extern void taudec_(int *kto, int *npnt, float *heltau, float p4tau[4]);
extern void pyjoin_(int *njoin, int ijoin[]);
extern void pyshow_(int *ip1, int *ip2, double *qmax);
extern void pyprep_(int *ip);
extern int pychge_(int *kf);
extern void pylist_(int *flag);
extern void pyhepc_(int *flag);
extern void pyexec_();
};
#else
extern "C" {
extern void luhadr_(int *idrec, int *level, int *idebug, int *iret);
extern void tauint_(int *inut, int *iout, int *jak1,
int *jak2, int *itdkrc, int *keya1, float *xk0dec);
extern void lutuhl_(int *one, int *ihlon);
extern void lugive_(char* opt);
};
#endif
Int_t JAK1, JAK2, ITDKRC, KEYA1;
Float_t XK0DEC;
ClassImp(JSFHadronizer)
#if __PYTHIA_VERSION__ >= 6
typedef struct { // Common for JETSET random variables
Int_t mrpy[6];
Float_t rrpy[100];
} COMMON_PYDATR;
extern COMMON_PYDATR pydatr_;
#else
typedef struct { // Flag for parton shower on/off
Int_t NDoPartonShower;
} COMMON_LufragFlag;
extern COMMON_LufragFlag lufragflag_;
typedef struct {
Int_t mdcy[3][500], mdme[2][2000];
Float_t brat[2000];
Int_t kfdp[5][2000];
} COMMON_LUDAT3;
extern COMMON_LUDAT3 ludat3_;
typedef struct { // Common for JETSET random variables
Int_t mrlu[6];
Float_t rrlu[100];
} COMMON_LUDATR;
extern COMMON_LUDATR ludatr_;
#endif
#ifndef __LCLIBRAN_USE_RANMAR__
typedef struct { // Commo for Tauola random variables
Float_t u[98];
Int_t ij97[2];
} COMMON_RASET1;
extern COMMON_RASET1 raset1_;
#endif
//_____________________________________________________________________________
JSFHadronizer::JSFHadronizer(const char *name, const char *title)
: JSFFULLGenerator(name,title)
{
// Constructor of JSFHadronizer. In this constructor, JSFSpring
// module is searched from the defined JSF module. last JSFModule
// which is inherited from JSFSpring class is a class whose Parton
// data is hadronized.
fSpring=0;
fCopySpringClassDataToBank=kTRUE;
TList *list=gJSF->Modules();
TIter next(list);
JSFModule *module;
while( (module=(JSFModule*)next()) ){
if( module->InheritsFrom("JSFSpring") ){
if( fSpring ) {
printf("Find module %s , inherited from JSFSpringn",fSpring->GetName());
printf("More than one JSFSpring are defined.");
printf(" Specify correspinding Hadronizer explicitylyn");
}
fSpring=(JSFSpring*)module;
}
}
if( !fSpring ){ Error("JSFHadronizer","No JSFSpring class was found"); }
JAK1 = gJSF->Env()->GetValue("JSFHadronizer.JAK1",0);
JAK2 = gJSF->Env()->GetValue("JSFHadronizer.JAK2",0);
ITDKRC = gJSF->Env()->GetValue("JSFHadronizer.ITDKRC",1);
KEYA1 = gJSF->Env()->GetValue("JSFHadronizer.KEYA1",1);
fIHLON = gJSF->Env()->GetValue("JSFHadronizer.IHLON",1);
sscanf(gJSF->Env()->GetValue("JSFHadronizer.XK0DEC","0.001"),"%g",&XK0DEC);
fDebug = gJSF->Env()->GetValue("JSFHadronizer.DebugFlag",0);
#if __PYTHIA_VERSION__ >= 6
fDoesShower = gJSF->Env()->GetValue("JSFHadronizer.PartonShower",1);
fPythia = new TPythia6();
#else
lufragflag_.NDoPartonShower = gJSF->Env()->GetValue("JSFHadronizer.PartonShower",1);
#endif
}
//_____________________________________________________________________________
JSFHadronizer::~JSFHadronizer()
{
#if __PYTHIA_VERSION__ >= 6
delete fPythia;
#endif
}
//_____________________________________________________________________________
Bool_t JSFHadronizer::Initialize()
{
// Initialize tauola
Int_t inut=5;
Int_t iout=6;
tauint_(&inut, &iout, &JAK1, &JAK2, &ITDKRC, &KEYA1, &XK0DEC);
#if __PYTHIA_VERSION__ >= 6
fPythia->SetMDCY(23,1,0);
fPythia->SetMDCY(24,1,0);
fPythia->SetMDCY(33,1,0);
#else
Int_t one=1;
lutuhl_(&one, &fIHLON);
ludat3_.mdcy[0][22]=0;
ludat3_.mdcy[0][23]=0;
ludat3_.mdcy[0][32]=0;
#endif
return kTRUE ;
}
//_____________________________________________________________________________
Bool_t JSFHadronizer::Process(Int_t ev)
{
// Create JSFGeneratorParticle from JSFSpringParton.
// For a time bing, first put JSFSpringParton in to TBS bank,
// calls luhadr fortran program, and put results into JSFGeneratorParticle
//
JSFFULLGenerator::Process(ev);
Int_t i;
#if __PYTHIA_VERSION__ >= 6
for(i=0;i<6;i++){ fMRPY[i]=pydatr_.mrpy[i];}
for(i=0;i<100;i++){ fRRPY[i]=pydatr_.rrpy[i];}
#else
for(i=0;i<6;i++){ fMRLU[i]=ludatr_.mrlu[i];}
for(i=0;i<100;i++){ fRRLU[i]=ludatr_.rrlu[i];}
#endif
#ifndef __LCLIBRAN_USE_RANMAR__
for(i=0;i<98;i++){ fRASET1U[i]=raset1_.u[i];}
for(i=0;i<2;i++){ fRASET1IJ97[i]=raset1_.ij97[i];}
#endif
Int_t iret;
#if __PYTHIA_VERSION__ >= 6
Hadronize(fSpring, iret);
#else
if( fCopySpringClassDataToBank ) TBPUT(fSpring);
Int_t idrec=1;
Int_t level=1;
Int_t idebug=0;
luhadr_(&idrec, &level, &idebug, &iret);
TBGET(); // Put hadronized results into Generator class.
// Remove Generator banks, since they are not necessary anymore.
gJSFLCFULL->TBDELB(1,"Spring:Parton_List",iret);
gJSFLCFULL->TBDELB(1,"Spring:Header",iret);
gJSFLCFULL->TBDELB(1,"Generator:Particle_List",iret);
gJSFLCFULL->TBDELB(1,"Generator:Header",iret);
#endif
if ( iret < 0 ) return kFALSE;
if( fDebug > 0 ) {
printf(" End of JSFHadronizer::Process()\n");
}
return kTRUE;
}
//_____________________________________________________________________________
Bool_t JSFHadronizer::EndRun()
{
Int_t i;
#if __PYTHIA_VERSION__ >= 6
for(i=0;i<6;i++){ fMRPY[i]=pydatr_.mrpy[i];}
for(i=0;i<100;i++){ fRRPY[i]=pydatr_.rrpy[i];}
#else
for(i=0;i<6;i++){ fMRLU[i]=ludatr_.mrlu[i];}
for(i=0;i<100;i++){ fRRLU[i]=ludatr_.rrlu[i];}
#endif
#ifndef __LCLIBRAN_USE_RANMAR__
for(i=0;i<98;i++){ fRASET1U[i]=raset1_.u[i];}
for(i=0;i<2;i++){ fRASET1IJ97[i]=raset1_.ij97[i];}
#endif
if( fFile->IsWritable() ) {
if( !JSFFULLGenerator::EndRun() ) return kFALSE;
Write();
}
return kTRUE;
}
//_____________________________________________________________________________
Bool_t JSFHadronizer::GetLastRunInfo()
{
// Read seed of previous run
Read(GetName());
printf("Random seeds for JSFHadronizer were reset by ");
printf("values from a file.n");
return kTRUE;
}
#if __ROOT_FULLVERSION__ >=30000
//___________________________________________________________________________
void JSFHadronizer::Streamer(TBuffer &R__b)
{
// Stream an object of class JSFHadronizer.
if (R__b.IsReading()) {
UInt_t R__s, R__c;
Version_t R__v=R__b.ReadVersion(&R__s, &R__c);
if( R__v > 1 ) {
JSFHadronizer::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c);
return;
}
JSFFULLGenerator::Streamer(R__b);
#if __PYTHIA_VERSION__ >= 6
R__b.ReadStaticArray(fMRPY);
R__b.ReadStaticArray(fRRPY);
#else
R__b.ReadStaticArray(fMRLU);
R__b.ReadStaticArray(fRRLU);
#endif
#ifndef __LCLIBRAN_USE_RANMAR__
R__b.ReadStaticArray(fRASET1U);
R__b.ReadStaticArray(fRASET1IJ97);
#endif
R__b.CheckByteCount(R__s, R__c, JSFHadronizer::IsA());
} else {
JSFHadronizer::Class()->WriteBuffer(R__b, this);
}
}
#else
//______________________________________________________________________________
void JSFHadronizer::Streamer(TBuffer &R__b)
{
// Stream an object of class JSFHadronizer.
if (R__b.IsReading()) {
Version_t R__v = R__b.ReadVersion(); if (R__v) { }
JSFFULLGenerator::Streamer(R__b);
#if __PYTHIA_VERSION__ >= 6
R__b.ReadStaticArray(fMRPY);
R__b.ReadStaticArray(fRRPY);
#else
R__b.ReadStaticArray(fMRLU);
R__b.ReadStaticArray(fRRLU);
#endif
#ifndef __LCLIBRAN_USE_RANMAR__
R__b.ReadStaticArray(fRASET1U);
R__b.ReadStaticArray(fRASET1IJ97);
#endif
} else {
R__b.WriteVersion(JSFHadronizer::IsA());
JSFFULLGenerator::Streamer(R__b);
#if __PYTHIA_VERSION__ >= 6
R__b.WriteArray(fMRPY, 6);
R__b.WriteArray(fRRPY, 100);
#else
R__b.WriteArray(fMRLU, 6);
R__b.WriteArray(fRRLU, 100);
#endif
#ifndef __LCLIBRAN_USE_RANMAR__
R__b.WriteArray(fRASET1U, 98);
R__b.WriteArray(fRASET1IJ97, 2);
#endif
}
}
#endif
#if __PYTHIA_VERSION__ >= 6
// ____________________________________________________________
void JSFHadronizer::Hadronize(JSFSpring *spring, Int_t &nret)
{
//(Function)
// Hadronize partons in JSFSpringParton objects
//(Input)
// spring: pointer to the input spring objects
//(Output)
// nret : =0 on sucess
// : =-1 on error
//CC**********************************************************************
//C*
//C*===========================================----===
//C* Subroutine LUHADR( IDREC, LEVEL, IDEBUG, NRET )
//C*===========================================----===
//C*
//C* (Function)
//C* LUND hadronizer.
//C* Hadronize particle in 'Spring:Particle_List' bank.
//C* (Inputs)
//C* IDREC ; record ID, use only ID=1
//C* LEVEL ; not used.
//C* IDEBUG ; not used.
//C* (Output)
//C* NRET ; Return code.
//C* (Update Record)
//C* 91/08/27 K.Fujii Derived from LUHADR by TAYM.
//C* Heavily modified.
//C* 95/01/11 K.Fujii Another big changes for implementation
//C* of tau polarization.
//C* ID=220000 is the LSP.
//C*
//CC**********************************************************************
// SUBROUTINE LUHADR( IDREC, LEVEL, IDEBUG, NRET )
// INTEGER*4 IDREC, LEVEL, IDEBUG, NRET
// COMMON /SSWORK/ IHEAD(20), NOELM(30), RBUF(20,30), RINLST(10,30),
// . ROTLST(20, 400), KSTAT(30), RTMP(20),
// . ISHUFL(30), ISHPR(2,30), ITMP(20),
// . NSGRP(10), ISGRP(2,30,10), JSTAT(30),
// . NPTLV(10), IPTLV(30,10), IDIDIT(30),
// . NDKPTN(10), IDKPTN(30,10), INOSHW(30)
// COMMON /LUSEED/ NLSEED
// DATA ILEVT/0/
//C
//C========< Entry Point >================================================
//C
//C--
//C Skip if non-event record.
//C--
// NRET = 0
// IF ( IDREC.NE.1 ) RETURN
// ILSAVE = NLSEED
//C--
//C Get 'Spring:Parton_List' into local common.
//C--
//C The parton showering is handled level by level, since showers
//C from parent parton pairs might modify the momenta of their
//C daughter partons.
//C The level is defined by
//C RBUF(19,*) = ISLEV*100 + ISPAIR
//C where partons with the same ISLEV are passed to LUFRAG at one
//C time (that means they comprise a set of color singlet states and
//C that there is no parent-daughter relation among the set) and
//C partons with the same ISPAIR are paired in LUSHOW calls in
//C LUFRAG.
//C Lepton pairs are required to have ISLEV = 0 and NDAU = 0 for
//C proper polarization handling when the pair contains a tau.
//C If you are satisfied with the default tau decay treatment, you
//C do not have to distinguish leptons from quarks in the decay
//C daughters from W or Z, etc.
//C Color singlet groups are distinguished by
//C RBUF(18,*) = ICF
//C where ICF must have different values for different color singlet
//C groups.
//C Helicities of the primary partons should be stored in
//C RBUF(17,*) = IHEL
//C when necessary as with tau leptons.
//C--
// CALL TBNOEL(1, 'Spring:Parton_List', NPART, NOELM )
// IF ( NPART.GT.30 ) THEN
// PRINT *,'%Fatal Error in LUHADR...'
// PRINT *,'# of partons in ''SPring:Parton_list'' ',
// . ' bank exceeds buffer size.'
// STOP
// ENDIF
// DO 100 IP = 1, NPART
// IE = NOELM(IP)
// CALL TBGET(1,'Spring:Parton_List',IE, NW, RBUF(1,IP), IRET )
//100 CONTINUE
JSFGeneratorBuf *gbuf=(JSFGeneratorBuf*)EventBuf();
TClonesArray &part=*(gbuf->GetParticles());
nret=0;
Int_t npgen=0;
JSFSpringBuf *spevt=(JSFSpringBuf*)spring->EventBuf();
TClonesArray *ps=spevt->GetPartons();
Int_t npart=spevt->GetNpartons();
Float_t rbuf[30][20];
for(Int_t j=0;j<npart;j++){
JSFSpringParton *p=(JSFSpringParton*)ps->At(j);
rbuf[j][0]=p->fSer;
rbuf[j][1]=p->fID;
rbuf[j][2]=p->fMass;
rbuf[j][3]=p->fCharge;
rbuf[j][4]=p->fP[1];
rbuf[j][5]=p->fP[2];
rbuf[j][6]=p->fP[3];
rbuf[j][7]=p->fP[0];
rbuf[j][11]=p->fNdaughter;
rbuf[j][12]=p->fFirstDaughter;
rbuf[j][13]=p->fMother;
rbuf[j][16]=p->fHelicity;
rbuf[j][17]=p->fColorID;
rbuf[j][18]=p->fShowerInfo;
}
//C--
//C Scan through Spring:Parton_List and classify them according to
//C required showing scheme.
//C--
//C NSLVL = # levels of showering
//C NPTLV(ISLV) = # partons in the level ISLV(=1,..,NSLVL)
//C IPTLV(IQ,ISLV) = IP pointer to RBUF of IQ-th parton in
//C the ISLV-th level (IQ=1,...,NPTLV(ISLV))
//C NSGRP(ISLV) = # showering pairs in the ISLV-th level
//C ISGRP(1,IPSHW,ISLV) = IQ pointer of the 1st parton
//C (2,IPSHW,ISLV) = IQ pointer of the 2nd parton
//C in IPSHW-th pair (IPSHW=1,..,NSGRP(ISLV))
//C NDKPTN(ISLV) = # unstable partons (NDAU>0) in the
//C ISLV-th level.
//C IDKPTN(IMTH,ISLV) = IQ pointer of IMTH-th unstable parton
//C in the ISLV-th level
//C--
//C NNOSHW = # non-showering partons
//C INOSHW(IQ) = IP pointer of IQ-th non-showering
//C partons (IQ=1,..,NNOSHW)
//C Notice that non-showering partons might be affected by parton
//C showering, if they are daughters of some parton in a showering
//C pair: their momenta have to be readjusted when the parent
//C momentum is modified by the parton showering.
//C--
Int_t nslvl = -999;
Int_t nnoshw=0 ;
Int_t nptlv[10], iptlv[10][30], nsgrp[10], idkptn[10][30];
Int_t isgrp[10][30][2], ndkptn[10];
for(Int_t i=0;i<10;i++){
nptlv[i]=0;
nsgrp[i]=0;
ndkptn[i]=0;
for(Int_t k=0;k<30;k++){
iptlv[i][k]=0;
isgrp[i][k][0]=0;
isgrp[i][k][1]=0;
}
}
Int_t inoshw[30];
//C--
for(Int_t ip=1;ip<=npart;ip++){
Int_t ndau = (Int_t)rbuf[ip-1][11];
Int_t ifshw = (Int_t)rbuf[ip-1][18];
Int_t islv = ifshw/100;
Int_t ipshw = ifshw%100;
if( fDebug > 0 ) {
printf(" Looking for shower pairn");
printf(" ip=%d ndau=%d ifshw=%d islv=%d ipshw=%dn",
ip, ndau, ifshw, islv, ipshw);
}
nslvl = TMath::Max(nslvl, islv);
if( islv <= 0 ) {
if ( ndau <= 0 ) {
nnoshw = nnoshw + 1 ;
inoshw[nnoshw-1] = ip;
}
}
else {
nptlv[islv-1]++;
iptlv[islv-1][nptlv[islv-1]-1] = ip;
if( ipshw > 0 ) {
if( nsgrp[islv-1] > 0 && isgrp[islv-1][ipshw-1][0] > 0 ) {
isgrp[islv-1][ipshw-1][1] = nptlv[islv-1];
}
else {
nsgrp[islv-1]++;
isgrp[islv-1][ipshw-1][0]=nptlv[islv-1];
}
if( ndau > 0 ) {
ndkptn[islv-1] ++;
idkptn[islv-1][ndkptn[islv-1]-1]=nptlv[islv-1];
}
}
}
}
if( fDebug > 0 ) {
printf(" nslvl=%d n",nslvl);
for(Int_t klvl=0;klvl<nslvl;klvl++){
printf(" nsgrp[%d]=%d",klvl,nsgrp[klvl]);
printf(" isgrp[%d][0][0:1]=%d,%d",klvl,isgrp[klvl][0][0],isgrp[klvl][0][1]);
}
printf("n");
}
//C--
//C Now start creating Generator:Particle_List bank.
//C--
//C Data format of Generator:Particle_List bank.
//C Generator:Particle_List
//C Elm# n RBUF( 1) = particle serial #
//C ( 2) = Particle ID (PDG).
//C ( 3) = mass (GeV)
//C ( 4) = charge
//C ( 5) = Px (GeV)
//C ( 6) = Py (GeV)
//C ( 7) = Pz (GeV)
//C ( 8) = E (GeV)
//C ( 9) = X (cm) of vertex point.
//C (10) = Y (cm)
//C (11) = Z (cm)
//C (12) = # of daughter particle (=0 in this hadronizer)
//C (13) = Pointer to the daughter (=0)
//C (14) = Pointer to the mother (-n n=element# of
//C SPRING:Parton_List )
//C (15) = trigger timing.
//C (16) = particle life time (c*nsec)
//C (18) = not used
//C (19) = not used
//C (20) = not used.
//C--
//C Loop over levels of parton showers.
//C--
Double_t rinlst[30][10];
Double_t rotlst[400][20];
Int_t inelm=0;
Int_t ididit[30];
Int_t ishpr1[30], ishpr2[30];
Int_t kstat[30], jstat[30], ishufl[30];
Int_t nout;
for(Int_t islv=1;islv<=nslvl;islv++){
if( nptlv[islv-1] < 0 ) { continue; }
Int_t nin=0;
for(Int_t k=0;k<30;k++){ ididit[k]=0; }
//C--
//C Loop over partons in this level and store their information
//C in RINLST and invoke LUFRAG to carry out parton showering
//C and fragmentation and decay.
//C--
//C NIN = # partons in RINLST
//C KSTAT(IIN) = color flow flag
//C NSG = NSGRP(ISLV) = # showering pairs in the
//C current showering level
//C ISHPR(1,IS) = IIN pointer of the 1st parton
//C (2,IS) = IIN pointer of the 2nd parton
//C in IS-th showering pair (IS=1,..,NSG)
//C RINLST(1,IIN) = ID
//C (2,IIN) = Px
//C (3,IIN) = Py
//C (4,IIN) = Pz
//C (5,IIN) = E
//C (6,IIN) = mass
//C (7,IIN) = helicity
//C--
//C Partons in RINLST must be orderd to make a planar tree in
//C terms of color flow. (e.g.) When this level consists of b,
//C W+, bbar, and W- from a t-tbar pair, they must be ordered,
//C for instance, as b, bbar, W+, and W- and
//C KSTAT(1-4) = 1 0 0 0
//C JSTAT(1-4) = 0 0 2 2
//C--
//C ISHUFL(IQ) = Pointer to RINLST of IQ-th parton in the
//C current showering level (IIN=1,..,NIN)
//C--
for(Int_t iq=1;iq<=nptlv[islv-1];iq++){
if( ididit[iq-1] > 0 ) { continue; }
Int_t ip = iptlv[islv-1][iq-1];
Int_t id = (Int_t)rbuf[ip-1][1];
Int_t ida = TMath::Abs(id);
Int_t ndau = (Int_t)rbuf[ip-1][11];
nin++;
Int_t icf = (Int_t)rbuf[ip-1][17];
for(Int_t k=0;k<10;k++){ rinlst[nin-1][k]=0.0; }
rinlst[nin-1][0] = rbuf[ip-1][1];
rinlst[nin-1][1] = rbuf[ip-1][4];
rinlst[nin-1][2] = rbuf[ip-1][5];
rinlst[nin-1][3] = rbuf[ip-1][6];
rinlst[nin-1][4] = rbuf[ip-1][7];
rinlst[nin-1][5] = rbuf[ip-1][2];
rinlst[nin-1][6] = rbuf[ip-1][16];
ishufl[iq-1] = nin;
ididit[iq-1] = 1;
kstat[nin-1] = 0;
jstat[nin-1] = ndau;
if( ida > 9 && ida != 21 ) { continue; }
kstat[nin-1] = 1;
//C--
//C Look for color singlet partner.
//C--
for(Int_t jq=iq+1;jq<=nptlv[islv-1];jq++){
if( ididit[jq-1] > 0 ) { continue; }
Int_t jp = iptlv[islv-1][jq-1];
Int_t jd = (Int_t)rbuf[jp-1][1];
Int_t jda = TMath::Abs(jd);
Int_t ndau= (Int_t)rbuf[jp-1][11];
if( jda > 9 && jda != 21 ) { continue; }
Int_t jcf = (Int_t)rbuf[jp-1][17];
if( jcf != icf ) { continue; }
nin++;
for(Int_t k=0;k<10;k++){ rinlst[nin-1][k]=0; }
rinlst[nin-1][0] = rbuf[jp-1][1];
rinlst[nin-1][1] = rbuf[jp-1][4];
rinlst[nin-1][2] = rbuf[jp-1][5];
rinlst[nin-1][3] = rbuf[jp-1][6];
rinlst[nin-1][4] = rbuf[jp-1][7];
rinlst[nin-1][5] = rbuf[jp-1][2];
rinlst[nin-1][6] = rbuf[jp-1][16];
ishufl[jq-1] = nin ;
ididit[jq-1] = 1;
kstat[nin-1] = 1;
jstat[nin-1] = ndau;
}
kstat[nin-1] = 0;
}
//C--
//C Do parton shower and fragmentation.
//C--
Int_t nsg=nsgrp[islv-1];
for(Int_t is=1;is<=nsg;is++){
ishpr1[is-1]=ishufl[isgrp[islv-1][is-1][0]-1];
ishpr2[is-1]=ishufl[isgrp[islv-1][is-1][1]-1];
}
Int_t maxout=400;
Int_t iret;
if( fDebug > 0 ) {
printf(" Before JSFHadronize::Fragmentaion call n");
printf(" nsg=%d ishpr=",nsg);
for(Int_t k=0;k<nsg;k++){
printf(" (%d,%d) ",ishpr1[k],ishpr2[k]);
}
printf("n");
}
Fragmentation(nin, rinlst, maxout, nsg, ishpr1, ishpr2,
kstat, jstat, nout, rotlst, iret);
if( iret < 0 ) {
nret = -1;
return;
}
if( fDebug > 0 ) {
printf(" After Fragmentation .. nout=%dn",nout);
}
//C--
//C Put particles in this level to bank.
//C--
for(Int_t i=1;i<=nout;i++){
rotlst[i-1][0] = i + inelm;
rotlst[i-1][12] = (Int_t)(rotlst[i-1][12]+0.1) + inelm;
if( rotlst[i-1][13] > 0 ) {
rotlst[i-1][13] = (Int_t)(rotlst[i-1][13]+0.1) + inelm;
}
else {
for(Int_t j=1;j<=nptlv[islv-1];j++){
if( jstat[j-1] == 0 ) {
rotlst[i-1][13] = -iptlv[islv-1][j-1];
break;
}
}
}
Int_t ie=i+inelm;
new(part[ie-1]) JSFGeneratorParticle(&rotlst[i-1][0]);
npgen++;
}
inelm += nout;
//C--
//C Modify daughter momenta if parent momentum got modified.
//C--
for(Int_t imth=1;imth<=ndkptn[islv-1];imth++){
Int_t iq = idkptn[islv-1][imth-1];
Int_t ip = iptlv[islv-1][iq-1];
Int_t in = ishufl[iq-1];
Int_t ndau= (Int_t)rbuf[ip-1][11];
Int_t idau= (Int_t)rbuf[ip-1][12]-1;
TLorentzVector ppar(rinlst[in-1][1], rinlst[in-1][2],
rinlst[in-1][3], rinlst[in-1][4]);
TVector3 bpar=-ppar.BoostVector();
for(Int_t jdau=1;jdau<=ndau;jdau++){
TLorentzVector pdau(rbuf[idau][4], rbuf[idau][5],
rbuf[idau][6], rbuf[idau][7]);
TLorentzVector pref(rbuf[ip-1][4], rbuf[ip-1][5],
rbuf[ip-1][6], rbuf[ip-1][7]);
TVector3 b=pref.BoostVector();
pdau.Boost(b);
pdau.Boost(bpar);
rbuf[4][idau]=pdau.Px();
rbuf[5][idau]=pdau.Py();
rbuf[6][idau]=pdau.Pz();
rbuf[7][idau]=pdau.E();
idau = idau+1;
}
}
}
//C--
//C Then output non-showering particles.
//C Notice that IDA = 220000 is the LSP(chi^0_1) which is
//C assumed to be stable.
//C--
Double_t rtmp[20];
for(Int_t jp=1;jp<=nnoshw;jp++){
Int_t ip = inoshw[jp-1];
Int_t id = (Int_t)rbuf[ip-1][1];
Int_t ida = TMath::Abs(id);
if( ida > 24 && ida != 220000 && ida != 1000022 ) {
printf("Warning in JSFHadronizer::Hadronize");
printf(" Particle ID=%d is not recognized.n",id);
continue;
}
//C--
//C Neutrinos, e, mu, gamma, and LSP.
//C--
if( ida==12 || ida==14 || ida==16 ||
ida==11 || ida==13 || ida==22 ||
ida==220000 || ida==1000022 ) {
inelm++;
for(Int_t k=0;k<20;k++){ rtmp[k]=0; }
rtmp[0]=inelm;
rtmp[1]=id;
rtmp[2]=rbuf[ip-1][2];
rtmp[3]=-id%2;
rtmp[4]=rbuf[ip-1][4];
rtmp[5]=rbuf[ip-1][5];
rtmp[6]=rbuf[ip-1][6];
rtmp[7]=rbuf[ip-1][7];
rtmp[13]=-ip;
new(part[inelm-1]) JSFGeneratorParticle(rtmp);
npgen++;
}
//C--
//C Tau, Z, and W.
//C--
else if( ida==15 || ida==23 || ida==24 ) {
for(Int_t k=0;k<10;k++){ rinlst[k][0]=0; }
rinlst[0][0] = rbuf[ip-1][1];
rinlst[0][1] = rbuf[ip-1][4];
rinlst[0][2] = rbuf[ip-1][5];
rinlst[0][3] = rbuf[ip-1][6];
rinlst[0][4] = rbuf[ip-1][7];
rinlst[0][5] = rbuf[ip-1][2];
rinlst[0][6] = rbuf[ip-1][16];
kstat[0] = 0 ;
jstat[0] = 0 ;
Int_t nin = 1;
Int_t n400=400;
Int_t zero=0;
Int_t iret;
// When 4th argument is zero, parton shower is not performed.
// 5th and 6th argumnet is dummy
Fragmentation(nin, rinlst, n400, zero, ishpr1, ishpr2,
kstat, jstat, nout, rotlst, iret);
if( iret < 0 ) {
nret = -1;
return;
}
for(Int_t i=1;i<=nout;i++){
rotlst[i-1][0] = i+inelm ;
rotlst[i-1][12] = (Int_t)(rotlst[i-1][12]+0.1) + inelm;
if( rotlst[i-1][13] > 0.0 ) {
rotlst[i-1][13] = (Int_t)(rotlst[i-1][13]+0.1) + inelm;
}
else {
rotlst[i-1][13] = -ip;
}
Int_t ie=i+inelm;
new(part[ie-1]) JSFGeneratorParticle(&rotlst[i-1][0]);
npgen++;
}
inelm = inelm + nout;
}
}
gbuf->SetNparticles(npgen);
if( fDebug > 0 ) {
printf(" JSFHadronizer::Hadronize() .. Number of generator particles");
printf(" is %dn",npgen);
}
}
//______________________________________________________________
void JSFHadronizer::Fragmentation(const Int_t nin, Double_t inlist[][10],
const Int_t mxxout, const Int_t nspar,
const Int_t ispar1[], const Int_t ispar2[],
const Int_t kstat[], const Int_t jstat[],
Int_t &nout, Double_t outlst[][20],
Int_t &nret)
{
/*
CC**********************************************************************
C*
C*==============================================================
C* Subroutine LUFRAG(NIN,INLIST,MXxOUT,NSPAR,ISPAR,KSTAT,JSTAT,
C* NOUT,OUTLST, NRET )
C*===================------------------=========================
C*
C* (Function)
C* Parton-shower, decay, and hadronize partons with
C* LUND-PS and string fragmentation.
C* (Inputs)
C* NIN :(I*4): # of input partons.
C* INLIST( 1,i) :(R*4): Particle ID (PDG format).
C* ( 2,i) :(R*4): Px (GeV).
C* ( 3,i) :(R*4): Py (GeV).
C* ( 4,i) :(R*4): Pz (GeV).
C* ( 5,i) :(R*4): E (GeV).
C* ( 6,i) :(R*4): m (GeV).
C* ( 7,i) :(R*4): helicity.
C* MXxOUT :(I*4): Size of OUTLST.
C* NSPAR :(I*4): # showering pairs.
C* ISPAR ( 1,j) :(I*4): pointer to 1-st parton in j-th pair.
C* ( 2,j) :(I*4): pointer to 2-nd parton in j-th pair.
C* KSTAT ( i) :(I*4): color singlet group flag.
C* JSTAT ( i) :(I*4): # daughters.
C* (Outputs)
C* NOUT :(I*4): # output particles.
C* OUTLST( 1,i) :(R*4): 0.
C* ( 2,i) :(R*4): Particle ID ( PDG format)
C* ( 3,i) :(R*4): mass (GeV)
C* ( 4,i) :(R*4): charge
C* ( 5,i) :(R*4): Px
C* ( 6,i) :(R*4): Py
C* ( 7,i) :(R*4): Pz
C* ( 8,i) :(R*4): E
C* ( 9,i) :(R*4): 0.
C* (10,i) :(R*4): 0.
C* (11,i) :(R*4): 0.
C* (12,i) :(R*4): # of daughter
C* (13,i) :(R*4): Pointer to the daughter (relative)
C* (14,i) :(R*4): Pointer to the mother (relative)
C* (15,i) :(R*4): 0.
C* (16,i) :(R*4): life time (c*t (cm))
C* (17,i) :(R*4): 0.
C* (18,i) :(R*4): 0.
C* (19,i) :(R*4): 0.
C* (20,i) :(R*4): 0.
C* NRET :(I*4): return flag.
C* (Update Records)
C* 95/02/10 K.Fujii New version for JETSET7.3.
C*
CC**********************************************************************
*/
Pyjets_t *pyjets=fPythia->GetPyjets();
//C--
//C Pointer from /LUJETS/ to OUTLST.
//C IDHIST ( i ) : Address in OUTLST for particles in /LUJETS/
//C--
Int_t idhist[4000];
for(Int_t k=0;k<4000;k++){ idhist[k]=0; }
Int_t njoin, ijoin[4000];
//C
//C ======< Entry Point >=============================================
//C
//C--
//C Reset #particles in /LUJETS/.
//C--
pyjets->N=0;
//C--
//C Check NSPAR to decide whether parton shower is required or not.
//C--
Int_t kf=0;
Int_t kfa=0;
Int_t two=2;
Int_t one=1;
if( nspar <= 0 ) {
kf=(Int_t)inlist[0][0];
kfa=TMath::Abs(kf);
if( nin==1 && kfa==15 ) {
//C--
//C Decay taus using TAUOLA.
//C Notice that TAUOLA initializatin is carried out in FLNPAR
//C and TAUOLA summary is printed out in FLUERU.
//C TAUOLA parameters should be read in by FLNPAR from FT05.
//C--
Int_t kto;
if( kf <= 0 ) {
kto=1;
}
else {
kto=2;
}
Float_t pol=inlist[0][6]*fIHLON;
Float_t tinlst[10];
for(Int_t i=0;i<10;i++){ tinlst[i]=(Float_t)inlist[0][i]; }
taudec_(&kto,&one, &pol, &tinlst[1]);
pyhepc_(&two);
if( fDebug > 0 ) {
printf(" Tau decay .. n");
pylist_(&one);
}
}
//C--
//C If input particles do not require parton shower generations,
//C fragment them according to the standard procedure.
//C--
else {
for(Int_t i=0;i<nin;i++){
pyjets->N++;
kf=(Int_t)inlist[i][0];
kfa=TMath::Abs(kf);
if( kfa == 220000) { kf=TMath::Sign(33,kf); }
Int_t ip=pyjets->N - 1;
pyjets->K[0][ip]=kstat[i]+1;
pyjets->K[1][ip]=kf;
pyjets->K[2][ip]=0;
pyjets->K[3][ip]=0;
pyjets->K[4][ip]=0;
pyjets->P[0][ip]=inlist[i][1];
pyjets->P[1][ip]=inlist[i][2];
pyjets->P[2][ip]=inlist[i][3];
pyjets->P[3][ip]=inlist[i][4];
pyjets->P[4][ip]=inlist[i][5];
pyjets->V[0][ip]=0;
pyjets->V[1][ip]=0;
pyjets->V[2][ip]=0;
pyjets->V[3][ip]=0;
pyjets->V[4][ip]=0;
}
}
}
else {
//C--
//C If input particles require parton shower generations,
//C Do PARTON SHOWER first.
//C--
njoin=0;
for(Int_t i=1;i<=nin;i++){
pyjets->N++;
kf=(Int_t)inlist[i-1][0];
kfa=TMath::Abs(kf);
if( kfa == 220000) { kf=TMath::Sign(33,kf); }
if( kstat[i-1] != 0 ) {
njoin++;
ijoin[njoin-1]=i;
}
else if ( njoin > 0 ) {
if( kstat[i-2] != 0 ) {
njoin++;
ijoin[njoin-1] = i;
}
}
Int_t ip=pyjets->N - 1;
pyjets->K[0][ip]=kstat[i-1]+1;
pyjets->K[1][ip]=kf;
pyjets->K[2][ip]=0;
pyjets->K[3][ip]=0;
pyjets->K[4][ip]=0;
pyjets->P[0][ip]=inlist[i-1][1];
pyjets->P[1][ip]=inlist[i-1][2];
pyjets->P[2][ip]=inlist[i-1][3];
pyjets->P[3][ip]=inlist[i-1][4];
pyjets->P[4][ip]=inlist[i-1][5];
pyjets->V[0][ip]=0;
pyjets->V[1][ip]=0;
pyjets->V[2][ip]=0;
pyjets->V[3][ip]=0;
pyjets->V[4][ip]=0;
}
//C--
//C Do LUJOIN.
//C--
pyjoin_(&njoin, ijoin);
if( fDebug > 0 ) {
printf(" After join njoin=%d ijoin=",njoin);
for(Int_t k=0;k>njoin;k++){
printf(" %d ",ijoin[k]);
}
printf("n");
}
//C--
//C Do LUSHOW for each showering pair.
//C--
if( fDoesShower) {
for(Int_t is=0;is<nspar;is++){
Int_t ip1=ispar1[is];
Int_t ip2=ispar2[is];
Double_t qmx=TMath::Power(pyjets->P[3][ip1-1]+pyjets->P[3][ip2-1],2)
- TMath::Power(pyjets->P[0][ip1-1]+pyjets->P[0][ip2-1],2)
- TMath::Power(pyjets->P[1][ip1-1]+pyjets->P[1][ip2-1],2)
- TMath::Power(pyjets->P[2][ip1-1]+pyjets->P[2][ip2-1],2);
qmx = TMath::Sqrt(qmx)-pyjets->P[4][ip1-1]-pyjets->P[4][ip2-1];
qmx = TMath::Max(qmx, 1.);
pyshow_(&ip1, &ip2, &qmx);
if( fPythia->GetMSTU(23) != 0 || fPythia->GetMSTU(24) != 0
|| fPythia->GetMSTU(28) != 0 ) {
printf(" Warning .. JSFHadronizer::Fragmentationn");
printf(" Possible error in PYSHOW detected.");
printf(" MSTU(23)=%d MSTU(24)=%d MSTU(28)=%dn",
fPythia->GetMSTU(23), fPythia->GetMSTU(24), fPythia->GetMSTU(28)),
printf(" ip1=%d,ip2=%dqmx=%g",ip1,ip2,qmx);
printf(" This event will be skipped.n");
fPythia->SetMSTU(23,0);
fPythia->SetMSTU(24,0);
fPythia->SetMSTU(28,0);
nret=-1;
return;
}
}
}
//C--
//C Store modified momenta for unstable particles.
//C When one such particle has a daughter of the same flavor,
//C take it as the one after gluon emission and test if further
//C gluon emission occured or not. If not, consider the last one
//C as the parent parton after gluon radiations.
//C--
//C MDF = the line# of the 1-st daughter
//C MDL = the line# of the last daughter
//C--
for(Int_t i=1;i<=nin;i++){
if( jstat[i-1] > 0 ) {
kf=pyjets->K[1][i-1];
Int_t mp=i;
label35:
Int_t mdf=pyjets->K[mp-1][3];
Int_t mdl=pyjets->K[mp-1][4];
if( mdf*mdl != 0 ) {
for(Int_t m=mdf;m<=mdl;m++){
if( pyjets->K[1][m-1] == kf ) {
mp=m;
if( pyjets->K[0][m-1] > 3 ) { goto label35; }
}
}
}
inlist[i-1][1]=pyjets->P[0][mp-1];
inlist[i-1][2]=pyjets->P[1][mp-1];
inlist[i-1][3]=pyjets->P[2][mp-1];
inlist[i-1][4]=pyjets->P[3][mp-1];
}
} // end of loop for(Int_t i=1;i<=nin;i++){
Int_t one=1;
pyprep_(&one);
}
//C--
//C Do LUEXEC for fragmentation.
//C-
pyexec_();
if( fDebug > 0 ) {
printf(" After pyexec...n");
pylist_(&one);
}
if( fPythia->GetMSTU(23) != 0 || fPythia->GetMSTU(24) != 0
|| fPythia->GetMSTU(28) != 0 ) {
printf(" Warning .. JSFHadronizer::Fragmentationn");
printf(" Possible error in PYEXEC detected.");
printf(" MSTU(23)=%d MSTU(24)=%d MSTU(28)=%d",
fPythia->GetMSTU(23), fPythia->GetMSTU(24), fPythia->GetMSTU(28)),
printf(" This event will be skipped.n");
fPythia->SetMSTU(23,0);
fPythia->SetMSTU(24,0);
fPythia->SetMSTU(28,0);
nret=-1;
return;
}
//C--
//C Output particles after fragmentations and decays.
//C--
nout=0;
if( fDebug > 0 ) {
printf(" After fragmentation...pyjets->N=%dn",pyjets->N);
}
for(Int_t i=1;i<=pyjets->N;i++){
kf=pyjets->K[1][i-1];
Int_t kh=pyjets->K[2][i-1];
kfa=TMath::Abs(kf);
idhist[i-1]=0;
//C--
//C Skip quarks and gluons.
//C--
if( ((kfa >= 1 && kfa <= 10)) || kfa==21 ) { continue; }
//C--
//C Skip fundamental bosons except for photons.
//C--
if( kfa >= 23 && kfa <= 100 ) { continue;}
//C--
//C Now store this particle in OUTLST.
//C--
nout++;
if( nout > mxxout ) {
printf(" Warning in JSFHadronizer::Fragment");
printf("Output buffer overflowed. Trancated");
printf(" nout=%dn",nout);
Int_t two=2;
pylist_(&two);
nout=mxxout;
goto label900;
}
for(Int_t k=0;k<20;k++){ outlst[nout-1][k]=0; }
idhist[i-1]=nout;
//C--
Double_t chrg=0.0;
Double_t xctau=0.0;
if( kfa == 33 ) {
chrg=0.0;
kf=TMath::Sign(220000,kf);
xctau=0.0;
}
else {
chrg = ((Double_t)pychge_(&kf))/3.0;
Int_t kc=fPythia->Pycomp(kf);
xctau=fPythia->GetPMAS(kc,4)*0.1;
}
//C--
outlst[nout-1][1]=kf;
outlst[nout-1][2]=pyjets->P[4][i-1];
outlst[nout-1][3]=chrg;
outlst[nout-1][4]=pyjets->P[0][i-1];
outlst[nout-1][5]=pyjets->P[1][i-1];
outlst[nout-1][6]=pyjets->P[2][i-1];
outlst[nout-1][7]=pyjets->P[3][i-1];
outlst[nout-1][11]=0.0;
//C--
if( kh == 0 ) {
outlst[nout-1][13]=0.0;
}
else {
Int_t ipar=idhist[kh-1];
if( ipar > 0 ) {
outlst[ipar-1][11] = outlst[ipar-1][11]+1.0;
Int_t ndau = (Int_t)(outlst[ipar-1][11]+0.1);
if( ndau == 1 ) { outlst[ipar-1][12] = nout; }
outlst[nout-1][13] = ipar;
}
}
outlst[nout-1][15]=xctau;
}
//C--
//C That's it.
//C--
label900:
if( fDebug > 0 ) {
printf(" End of JSFHadronizer::Fragmentation nout=%dn",nout);
}
return ;
}
#endif
ROOT page - Home page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.