3 #include "HepMC/GenEvent.h" 4 #include "HepMC/PdfInfo.h" 5 #include "HepMC/PythiaWrapper6_4.h" 8 #include "HepMC/HEPEVT_Wrapper.h" 9 #include "HepMC/IO_HEPEVT.h" 10 #include "HepMC/IO_GenEvent.h" 19 #include "HepPID/ParticleIDTranslations.hh" 30 #include "CLHEP/Random/RandomEngine.h" 47 if(
debug && ++call < 100)
cout<<
"dcasrn from c++, call: "<<call<<
" random number: "<<rdm_nb<<endl;
63 FortranCallback::getInstance()->fillHeader();
67 FortranCallback::getInstance()->fillEvent();
80 double p[5][pyjets_maxn],
v[5][pyjets_maxn];
90 fComEnergy(pset.getParameter<double>(
"comEnergy")),
91 fextCrossSection(pset.getUntrackedParameter<double>(
"crossSection",-1.)),
92 fextCrossSectionError(pset.getUntrackedParameter<double>(
"crossSectionError",-1.)),
93 fFilterEfficiency(pset.getUntrackedParameter<double>(
"filterEfficiency",-1.)),
95 fMaxEventsToPrint(pset.getUntrackedParameter<
int>(
"maxEventsToPrint",0)),
96 fHepMCVerbosity(pset.getUntrackedParameter<
bool>(
"pythiaHepMCVerbosity",
false)),
97 fPythiaListVerbosity(pset.getUntrackedParameter<
int>(
"pythiaPylistVerbosity",0)),
99 fDisplayPythiaBanner(pset.getUntrackedParameter<
bool>(
"displayPythiaBanner",
false)),
100 fDisplayPythiaCards(pset.getUntrackedParameter<
bool>(
"displayPythiaCards",
false)){
105 if(pset.
exists(
"doPDGConvert"))
106 fConvertToPDG = pset.
getParameter<
bool>(
"doPDGConvert");
113 <<
" Pythia did not accept MSTU(12)=12345";
122 <<
" Pythia did not accept MSTU(13)=0";
144 for(
int ip=0; ip<pyjets_maxn; ip++){
145 for(
int i=0;
i<5;
i++){
159 for(
int ip=0; ip<pyjets_maxn; ip++){
160 for(
int i=0;
i<5;
i++){
175 if(
event()->signal_process_id() <= 0)
event()->set_signal_process_id(
pypars.msti[0]);
183 if(pdf.id1() <= 0) pdf.set_id1(
pyint1.mint[14] == 21 ? 0 :
pyint1.mint[14]);
184 if(pdf.id2() <= 0) pdf.set_id2(
pyint1.mint[15] == 21 ? 0 :
pyint1.mint[15]);
185 if(pdf.x1() <= 0) pdf.set_x1(
pyint1.vint[40]);
186 if(pdf.x2() <= 0) pdf.set_x2(
pyint1.vint[41]);
187 if(pdf.pdf1() <= 0) pdf.set_pdf1(
pyint1.vint[38] /
pyint1.vint[40]);
188 if(pdf.pdf2() <= 0) pdf.set_pdf2(
pyint1.vint[39] /
pyint1.vint[41]);
189 if(pdf.scalePDF() <= 0) pdf.set_scalePDF(
pyint1.vint[50]);
191 event()->set_pdf_info(pdf) ;
194 cout<<
"standard Py6 event weight: pyint1.vint[96]: "<<
pyint1.vint[96]<<endl;
195 cout<<
"event weight returned by PYEVWT: 1./(pyint1.vint[98]): "<<1./(
pyint1.vint[98])<<endl;
205 event()->weights().push_back(1.);
224 for(HepMC::GenEvent::particle_iterator
part =
event()->particles_begin();
226 (*part)->set_pdg_id(HepPID::translatePythiatoPDT((*part)->pdg_id()));
235 cout <<
"Event process = "<<
pypars.msti[0]<<endl<<
"----------------------"<<endl;
255 if(
pyint1.mint[50] != 0 ){
286 int NPartsAfterDecays =
event().get()->particles_size();
287 int barcode = NPartsAfterDecays;
297 for(
int ipart=NPartsAfterDecays; ipart>NPartsBeforeDecays; ipart-- ){
300 int status = part->status();
301 if ( status != 1 )
continue;
303 int pdgid = part->pdg_id();
306 if(pydat3.mdcy[0][py6ptr-1] != 1)
continue;
308 int py6id = HepPID::translatePDTtoPythia(pdgid);
314 if(part->momentum().t() <= part->generated_mass()){
321 for(
int i=0;
i<5;
i++){
328 pyjets.k[1][0] = py6id;
329 pyjets.p[4][0] = part->generated_mass();
330 pyjets.p[3][0] = part->momentum().t();
331 pyjets.p[0][0] = part->momentum().x();
332 pyjets.p[1][0] = part->momentum().y();
333 pyjets.p[2][0] = part->momentum().z();
334 HepMC::GenVertex* prod_vtx = part->production_vertex();
335 if ( !prod_vtx )
continue;
336 pyjets.v[0][0] = prod_vtx->position().x();
337 pyjets.v[1][0] = prod_vtx->position().y();
338 pyjets.v[2][0] = prod_vtx->position().z();
339 pyjets.v[3][0] = prod_vtx->position().t();
352 for(
int iprt1=1; iprt1<pyjets.n; iprt1++){
354 part->set_status( 2 );
356 HepMC::GenVertex* DecVtx =
new HepMC::GenVertex( HepMC::FourVector(pyjets.v[0][iprt1],
359 pyjets.v[3][iprt1]) );
360 DecVtx->add_particle_in( part );
363 HepMC::FourVector pmom(pyjets.p[0][iprt1],pyjets.p[1][iprt1],
364 pyjets.p[2][iprt1],pyjets.p[3][iprt1] );
367 if(pyjets.k[0][iprt1] >= 1 && pyjets.k[0][iprt1] <= 10) dstatus = 1;
368 else if(pyjets.k[0][iprt1] >= 11 && pyjets.k[0][iprt1] <= 20) dstatus = 2;
369 else if(pyjets.k[0][iprt1] >= 21 && pyjets.k[0][iprt1] <= 30) dstatus = 3;
370 else if(pyjets.k[0][iprt1] >= 31 && pyjets.k[0][iprt1] <= 100) dstatus = pyjets.k[0][iprt1];
373 HepPID::translatePythiatoPDT(pyjets.k[1][iprt1]),
376 daughter->suggest_barcode( barcode );
377 DecVtx->add_particle_out( daughter );
380 for(iprt2=iprt1+1; iprt2<pyjets.n; iprt2++){
382 if(pyjets.k[2][iprt2] != parent){
383 parent = pyjets.k[2][iprt2];
387 HepMC::FourVector pmomN(pyjets.p[0][iprt2],pyjets.p[1][iprt2],
388 pyjets.p[2][iprt2],pyjets.p[3][iprt2] );
391 if(pyjets.k[0][iprt2] >= 1 && pyjets.k[0][iprt2] <= 10) dstatus = 1;
392 else if(pyjets.k[0][iprt2] >= 11 && pyjets.k[0][iprt2] <= 20) dstatus = 2;
393 else if(pyjets.k[0][iprt2] >= 21 && pyjets.k[0][iprt2] <= 30) dstatus = 3;
394 else if(pyjets.k[0][iprt2] >= 31 && pyjets.k[0][iprt2] <= 100) dstatus = pyjets.k[0][iprt2];
397 HepPID::translatePythiatoPDT( pyjets.k[1][iprt2] ),
400 daughterN->suggest_barcode( barcode );
401 DecVtx->add_particle_out( daughterN );
407 event().get()->add_vertex( DecVtx );
421 for(
int i=0;
i<5;
i++){
483 for(
unsigned i=0;
i<AllSets.size();++
i) {
485 string Set = AllSets[
i];
489 for(vector<string>::const_iterator itPara = Para_Set.begin(); itPara != Para_Set.end(); ++itPara) {
493 <<
" Cascade did not accept the following parameter: \""<<*itPara<<
"\""<<endl;
523 vector<int> pdg = _pdg;
524 for(
size_t i=0;
i<pdg.size();
i++){
525 int PyID = HepPID::translatePDTtoPythia(pdg[
i]);
529 ostringstream pyCard ;
530 pyCard <<
"MDCY(" << pyCode <<
",1)=0";
546 HepMC::GenEvent::vertex_const_iterator
vbegin =
event()->vertices_begin();
547 HepMC::GenEvent::vertex_const_iterator
vend =
event()->vertices_end();
548 HepMC::GenEvent::vertex_const_iterator vitr =
vbegin;
550 for(; vitr !=
vend; ++vitr){
552 HepMC::GenVertex::particle_iterator pbegin = (*vitr)->particles_begin(
HepMC::children);
553 HepMC::GenVertex::particle_iterator pend = (*vitr)->particles_end(
HepMC::children);
554 HepMC::GenVertex::particle_iterator pitr = pbegin;
556 for(; pitr != pend; ++pitr){
558 if((*pitr)->end_vertex())
continue;
559 if((*pitr)->status()!=1)
continue;
561 int pdgcode=
abs((*pitr)->pdg_id());
563 if ( pydat3.mdcy[0][
pycomp_(pdgcode)-1] !=1 )
continue;
565 double ctau = pydat2.pmas[3][
pycomp_(pdgcode)-1];
566 HepMC::FourVector mom = (*pitr)->momentum();
567 HepMC::FourVector vin = (*vitr)->position();
572 bool decayInRange =
false;
573 while(!decayInRange){
577 double proper_length = - ctau *
log(unif_rand);
578 double factor = proper_length/mom.m();
579 x = vin.x() + factor * mom.px();
580 y = vin.y() + factor * mom.py();
581 z = vin.z() + factor * mom.pz();
582 t = vin.t() + factor * mom.e();
586 if (
sqrt(x*x+y*y)>
pydat1.parj[72] || fabs(z)>
pydat1.parj[73]) decayInRange =
true;
589 else if (
pydat1.mstj[21]==3){
590 if (
sqrt(x*x+y*y+z*z)>
pydat1.parj[71]) decayInRange =
true;
598 HepMC::GenVertex* vdec =
new HepMC::GenVertex(HepMC::FourVector(x,y,z,t));
599 event()->add_vertex(vdec);
600 vdec->add_particle_in((*pitr));
614 double sigma_CMS = 0.001 *
caeffic.avgi;
615 double error_CMS = 0.001 *
caeffic.sd;
629 return "gen::Cascade2Hadronizer";
636 bool accepted =
true;
638 if(!strncmp(ParameterString.c_str(),
"KE",2))
639 caluco.ke = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
641 else if(!strncmp(ParameterString.c_str(),
"IRES(1)",7))
642 capar6.ires[0] = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
644 else if(!strncmp(ParameterString.c_str(),
"KP",2))
645 caluco.kp = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
647 else if(!strncmp(ParameterString.c_str(),
"IRES(2)",7))
648 capar6.ires[1] = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
650 else if(!strncmp(ParameterString.c_str(),
"NFRAG",5))
651 cainpu.nfrag = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
653 else if(!strncmp(ParameterString.c_str(),
"IPST",4))
654 cashower.ipst = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
656 else if(!strncmp(ParameterString.c_str(),
"IPSIPOL",7))
657 jpsi.ipsipol = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
663 else if(!strncmp(ParameterString.c_str(),
"IFPS",4))
664 cainpu.ifps = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
666 else if(!strncmp(ParameterString.c_str(),
"ITIMSHR",7))
667 casshwr.itimshr = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
669 else if(!strncmp(ParameterString.c_str(),
"IRUNAEM",7))
670 capar1.irunaem = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
672 else if(!strncmp(ParameterString.c_str(),
"IRUNA",5))
673 capar1.iruna = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
675 else if(!strncmp(ParameterString.c_str(),
"IQ2",3))
676 capar1.iq2 = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
678 else if(!strncmp(ParameterString.c_str(),
"IPRO",4))
679 capar1.ipro = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
681 else if(!strncmp(ParameterString.c_str(),
"NFLAV",5))
682 caluco.nflav = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
684 else if(!strncmp(ParameterString.c_str(),
"INTER",5))
685 cainpu.inter = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
687 else if(!strncmp(ParameterString.c_str(),
"IHFLA",5))
688 cahflav.ihfla = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
690 else if(!strncmp(ParameterString.c_str(),
"IRPA",4))
691 cascol.irpa = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
693 else if(!strncmp(ParameterString.c_str(),
"IRPB",4))
694 cascol.irpb = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
696 else if(!strncmp(ParameterString.c_str(),
"IRPC",4))
697 cascol.irpc = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
699 else if(!strncmp(ParameterString.c_str(),
"ICCFM",5))
700 casshwr.iccfm = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
702 else if(!strncmp(ParameterString.c_str(),
"IFINAL",6))
703 cainpu.ifinal = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
705 else if(!strncmp(ParameterString.c_str(),
"IGLU",4))
706 cagluon.iglu = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
708 else if(!strncmp(ParameterString.c_str(),
"IRspl",5))
709 casprre.irspl = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
711 else if(!strncmp(ParameterString.c_str(),
"PT2CUT",6))
712 captcut.pt2cut[
capar1.ipro-1] = atof(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
714 else if(!strncmp(ParameterString.c_str(),
"NCB",3))
715 integr.ncb = atoi(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
717 else if(!strncmp(ParameterString.c_str(),
"ACC1",4))
718 integr.acc1 = atof(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
720 else if(!strncmp(ParameterString.c_str(),
"ACC2",4))
721 integr.acc2 = atof(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
723 else if(!strncmp(ParameterString.c_str(),
"SCALFAF",7))
724 scalf.scalfaf = atof(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
726 else if(!strncmp(ParameterString.c_str(),
"SCALFA",6))
727 scalf.scalfa = atof(&ParameterString[strcspn(ParameterString.c_str(),
"=")+1]);
729 else accepted =
false;
737 cout<<
"----------------------------------"<<endl;
738 cout<<
"---- Cascade Parameters ----"<<endl;
739 cout<<
"----------------------------------"<<endl;
742 cout<<
"computation switch: "<<endl;
743 cout<<
"number of calls per iteration for bases: "<<
integr.ncb<<endl;
744 cout<<
"relative precision for grid optimisation: "<<
integr.acc1<<
" %"<<endl;
745 cout<<
"relative precision for integration: "<<
integr.acc2<<
" %"<<endl;
748 cout<<
"kinematics: "<<endl;
749 cout<<
"flavour code of beam 1: "<<
caluco.ke<<endl;
750 cout<<
"direct or resolved particle 1: "<<
capar6.ires[0]<<endl;
751 cout<<
"pz of incoming beam 1: "<<
cainpu.plepin<<
" GeV"<<endl;
752 cout<<
"flavour code of beam 2: "<<
caluco.kp<<endl;
753 cout<<
"direct or resolved particle 2: "<<
capar6.ires[1]<<endl;
754 cout<<
"pz of incoming beam 2: "<<
cainpu.ppin<<
" GeV"<<endl;
755 cout<<
"number of active flavors: "<<
caluco.nflav<<endl;
758 cout<<
"hard subprocess selection: "<<endl;
759 cout<<
"hard subprocess number (IPRO): "<<
capar1.ipro<<endl;
760 cout<<
"IPRO = 10, switch to select QCD process g* g* -> q qbar: "<<
cascol.irpa<<endl;
761 cout<<
"IPRO = 10, switch to select QCD process g* g -> g g: "<<
cascol.irpb<<endl;
762 cout<<
"IPRO = 10, switch to select QCD process g* q -> g q: "<<
cascol.irpc<<endl;
763 cout<<
"flavor of heavy quark produced (IPRO = 11, 504, 514): "<<
cahflav.ihfla<<endl;
765 cout<<
"use matrix element including J/psi (Upsilon) polarisation: "<<
jpsi.ipsipol<<endl;
766 cout<<
"pt2 cut in matrix element for massless partons: "<<
captcut.pt2cut[
capar1.ipro-1]<<
" GeV^2"<<endl;
769 cout<<
"parton shower and fragmentation: "<<endl;
770 cout<<
"switch for fragmentation: "<<
cainpu.nfrag<<endl;
771 cout<<
"switch for parton shower (0 off - 1 initial state - 2 final state - 3 initial & final state): "<<
cainpu.ifps<<endl;
772 cout<<
"switch for time like parton shower in initial state: "<<
casshwr.itimshr<<endl;
773 cout<<
"select CCFM (1) or DGLAP (0) evolution: "<<
casshwr.iccfm<<endl;
774 cout<<
"scale switch for final state parton shower: "<<
cainpu.ifinal<<endl;
775 cout<<
"scale factor for final state parton shower: "<<
scalf.scalfaf<<endl;
776 cout<<
"switch for proton remnant treatment: "<<
casprre.irspl<<endl;
777 cout<<
"keep track of intermediate state in parton shower: "<<
cashower.ipst<<endl;
778 cout<<
"mode of interaction for e p: "<<
cainpu.inter<<endl;
781 cout<<
"pdfs, couplings and scales: "<<endl;
782 cout<<
"switch for running alphaem: "<<
capar1.irunaem<<endl;
783 cout<<
"switch for running alphas: "<<
capar1.iruna<<endl;
785 cout<<
"scale factor for scale in alphas: "<<
scalf.scalfa<<endl;
787 cout<<
"path where updf are stored: "<<
caspdf.pdfpath<<endl;
790 cout<<
"center of mass energy, cross section, filter efficiency: "<<endl;
800 cout<<
"----------------------------------"<<endl;
801 cout<<
"---- Pythia6 Parameters ----"<<endl;
802 cout<<
"----------------------------------"<<endl;
805 cout<<
"charm mass: "<<pydat2.pmas[0][3]<<
" GeV (default value = 1.5 GeV)"<<endl;
806 cout<<
"bottom mass: "<<pydat2.pmas[0][4]<<
" GeV (default value = 4.8 GeV)"<<endl;
807 cout<<
"Higgs mass: "<<pydat2.pmas[0][24]<<
" GeV (default value = 115 GeV)"<<endl;
809 cout<<
"lambda QCD: "<<
pydat1.paru[111]<<
" GeV (default value = 0.25 GeV)"<<endl;
811 cout<<
"alphas behaviour: "<<
pydat1.mstu[110]<<
" (default value = 1)"<<endl;
812 cout<<
"nr of flavours wrt lambda QCD: "<<
pydat1.mstu[111]<<
" (default value = 5)"<<endl;
813 cout<<
"min nr of flavours for alphas: "<<
pydat1.mstu[112]<<
" (default value = 3)"<<endl;
814 cout<<
"max nr of flavours for alphas: "<<
pydat1.mstu[113]<<
" (default value = 5)"<<endl;
816 cout<<
"maximum angle settings: "<<
pydat1.mstj[47]<<
" (default value = 0)"<<endl;
void pythia6PrintParameters()
T getParameter(std::string const &) const
double dcasrn_(int *idummy)
Pythia6Service * fPy6Service
bool call_pygive(const std::string &line)
bool exists(std::string const ¶meterName) const
checks if a parameter exists
bool generatePartonsAndHadronize()
void setInternalXSec(const XSec &xsec)
std::auto_ptr< HepMC::GenEvent > & event()
void doSetRandomEngine(CLHEP::HepRandomEngine *v) override
static const std::string kPythia6
GenRunInfoProduct & runInfo()
void call_pylist(int mode)
static struct gen::@586 pyjets_local
Pythia6ServiceWithCallback(const edm::ParameterSet &pset)
edm::ParameterSet fParameters
void call_caend(int mode)
bool initializeForInternalPartons()
unsigned int fPythiaListVerbosity
bool fDisplayPythiaBanner
Abs< T >::type abs(const T &t)
bool declareStableParticles(const std::vector< int > &)
void setPYUPDAParams(bool afterPyinit)
const char * classname() const
bool initializeForExternalPartons()
static CLHEP::HepRandomEngine * cascade2RandomEngine
std::auto_ptr< GenEventInfoProduct > & eventInfo()
bool cascadeReadParameters(const std::string &ParameterString)
double fextCrossSectionError
static const std::string kFortranInstance
static FortranCallback * getInstance()
~Cascade2Hadronizer() override
void setRandomEngine(CLHEP::HepRandomEngine *v)
return(e1-e2)*(e1-e2)+dp *dp
void resetIterationsPerEvent()
void setExternalXSecNLO(const XSec &xsec)
void cascadePrintParameters()
void setExternalXSecLO(const XSec &xsec)
HepMC::IO_HEPEVT hepevtio
unsigned int fMaxEventsToPrint