|
|
|
#include <Cascade2Hadronizer.h>
Definition at line 23 of file Cascade2Hadronizer.h.
| gen::Cascade2Hadronizer::Cascade2Hadronizer | ( | edm::ParameterSet const & | ps | ) |
Definition at line 85 of file Cascade2Hadronizer.cc.
References gen::call_pygive(), edm::errors::Configuration, gather_cfg::cout, debug, Exception, edm::ParameterSet::exists(), fConvertToPDG, fDisplayPythiaBanner, fDisplayPythiaCards, fFlat, fFlat_extern, flushTmpStorage(), fParameters, and edm::ParameterSet::getParameter().
: BaseHadronizer(pset), fPy6Service(new Pythia6ServiceWithCallback(pset)), //-- this will store py6 parameters for further settings //-- fRandomEngine(&getEngineReference()), //-- defined in GeneratorInterface/Core/src/RNDMEngineAccess.cc //-- CLHEP::HepRandomEngine& gen::getEngineReference() //-- { edm::Service<edm::RandomNumberGenerator> rng; //-- return rng->getEngine(); } fComEnergy(pset.getParameter<double>("comEnergy")), fextCrossSection(pset.getUntrackedParameter<double>("crossSection",-1.)), fextCrossSectionError(pset.getUntrackedParameter<double>("crossSectionError",-1.)), fFilterEfficiency(pset.getUntrackedParameter<double>("filterEfficiency",-1.)), fMaxEventsToPrint(pset.getUntrackedParameter<int>("maxEventsToPrint",0)), fHepMCVerbosity(pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity",false)), fPythiaListVerbosity(pset.getUntrackedParameter<int>("pythiaPylistVerbosity",0)), fDisplayPythiaBanner(pset.getUntrackedParameter<bool>("displayPythiaBanner",false)), fDisplayPythiaCards(pset.getUntrackedParameter<bool>("displayPythiaCards",false)){ fParameters = pset.getParameter<ParameterSet>("Cascade2Parameters"); edm::Service<edm::RandomNumberGenerator> rng; if(debug) cout<<"seed: "<<rng->mySeed()<<endl; CLHEP::HepRandomEngine& engine = rng->getEngine(); fFlat = new CLHEP::RandFlat(engine); fFlat_extern = fFlat; fConvertToPDG = false; if(pset.exists("doPDGConvert")) fConvertToPDG = pset.getParameter<bool>("doPDGConvert"); //-- silence Pythia6 banner printout, unless display requested if(!fDisplayPythiaBanner){ if(!call_pygive("MSTU(12)=12345")){ throw edm::Exception(edm::errors::Configuration,"PythiaError") <<" Pythia did not accept MSTU(12)=12345"; } } //-- silence printouts from PYGIVE, unless display requested if(!fDisplayPythiaCards){ if(!call_pygive("MSTU(13)=0")){ throw edm::Exception(edm::errors::Configuration,"PythiaError") <<" Pythia did not accept MSTU(13)=0"; } } //-- tmp stuff to deal with EvtGen corrupting pyjets flushTmpStorage(); }
| gen::Cascade2Hadronizer::~Cascade2Hadronizer | ( | ) |
Definition at line 142 of file Cascade2Hadronizer.cc.
References fPy6Service.
{
if(fPy6Service != 0) delete fPy6Service;
}
| void gen::Cascade2Hadronizer::cascadePrintParameters | ( | ) |
Definition at line 741 of file Cascade2Hadronizer.cc.
References cagluon, cahflav, cainpu, caluco, capar1, capar6, captcut, cascol, cashower, caspdf, casprre, casshwr, gather_cfg::cout, fComEnergy, fextCrossSection, fextCrossSectionError, fFilterEfficiency, integr, jpsi, and scalf.
Referenced by initializeForInternalPartons().
{
cout<<""<<endl;
cout<<"----------------------------------"<<endl;
cout<<"---- Cascade Parameters ----"<<endl;
cout<<"----------------------------------"<<endl;
cout<<""<<endl;
cout<<"computation switch: "<<endl;
cout<<"number of calls per iteration for bases: "<<integr.ncb<<endl;
cout<<"relative precision for grid optimisation: "<<integr.acc1<<" %"<<endl;
cout<<"relative precision for integration: "<<integr.acc2<<" %"<<endl;
cout<<""<<endl;
cout<<"kinematics: "<<endl;
cout<<"flavour code of beam 1: "<<caluco.ke<<endl;
cout<<"direct or resolved particle 1: "<<capar6.ires[0]<<endl;
cout<<"pz of incoming beam 1: "<<cainpu.plepin<<" GeV"<<endl;
cout<<"flavour code of beam 2: "<<caluco.kp<<endl;
cout<<"direct or resolved particle 2: "<<capar6.ires[1]<<endl;
cout<<"pz of incoming beam 2: "<<cainpu.ppin<<" GeV"<<endl;
cout<<"number of active flavors: "<<caluco.nflav<<endl;
cout<<""<<endl;
cout<<"hard subprocess selection: "<<endl;
cout<<"hard subprocess number (IPRO): "<<capar1.ipro<<endl;
cout<<"IPRO = 10, switch to select QCD process g* g* -> q qbar: "<<cascol.irpa<<endl;
cout<<"IPRO = 10, switch to select QCD process g* g -> g g: "<<cascol.irpb<<endl;
cout<<"IPRO = 10, switch to select QCD process g* q -> g q: "<<cascol.irpc<<endl;
cout<<"flavor of heavy quark produced (IPRO = 11, 504, 514): "<<cahflav.ihfla<<endl;
// cout<<"select vector meson state: "<<jpsi.i23s<<endl;
cout<<"use matrix element including J/psi (Upsilon) polarisation: "<<jpsi.ipsipol<<endl;
cout<<"pt2 cut in matrix element for massless partons: "<<captcut.pt2cut[capar1.ipro-1]<<" GeV^2"<<endl;
cout<<""<<endl;
cout<<"parton shower and fragmentation: "<<endl;
cout<<"switch for fragmentation: "<<cainpu.nfrag<<endl;
cout<<"switch for parton shower (0 off - 1 initial state - 2 final state - 3 initial & final state): "<<cainpu.ifps<<endl;
cout<<"switch for time like parton shower in initial state: "<<casshwr.itimshr<<endl;
cout<<"select CCFM (1) or DGLAP (0) evolution: "<<casshwr.iccfm<<endl;
cout<<"scale switch for final state parton shower: "<<cainpu.ifinal<<endl;
cout<<"scale factor for final state parton shower: "<<scalf.scalfaf<<endl;
cout<<"switch for proton remnant treatment: "<<casprre.irspl<<endl;
cout<<"keep track of intermediate state in parton shower: "<<cashower.ipst<<endl;
cout<<"mode of interaction for e p: "<<cainpu.inter<<endl;
cout<<""<<endl;
cout<<"pdfs, couplings and scales: "<<endl;
cout<<"switch for running alphaem: "<<capar1.irunaem<<endl;
cout<<"switch for running alphas: "<<capar1.iruna<<endl;
cout<<"scale in alphas: "<<capar1.iq2<<endl;
cout<<"scale factor for scale in alphas: "<<scalf.scalfa<<endl;
cout<<"unintegrated pdf: "<<cagluon.iglu<<endl;
cout<<"path where updf are stored: "<<caspdf.pdfpath<<endl;
cout<<""<<endl;
cout<<"center of mass energy, cross section, filter efficiency: "<<endl;
cout<<"center of mass energy: "<<fComEnergy<<" GeV"<<endl;
cout<<"external LO cross section: "<<fextCrossSection<<" +/- "<<fextCrossSectionError<<" pb"<<endl;
cout<<"filter efficiency: "<<fFilterEfficiency<<endl;
cout<<""<<endl;
}
| bool gen::Cascade2Hadronizer::cascadeReadParameters | ( | const std::string & | ParameterString | ) |
Definition at line 641 of file Cascade2Hadronizer.cc.
References cagluon, cahflav, cainpu, caluco, capar1, capar6, captcut, cascol, cashower, casprre, casshwr, integr, jpsi, and scalf.
Referenced by initializeForInternalPartons().
{
bool accepted = 1;
if(!strncmp(ParameterString.c_str(),"KE",2))
caluco.ke = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IRES(1)",7))
capar6.ires[0] = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"KP",2))
caluco.kp = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IRES(2)",7))
capar6.ires[1] = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"NFRAG",5))
cainpu.nfrag = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IPST",4))
cashower.ipst = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IPSIPOL",7))
jpsi.ipsipol = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
//-- from version 2.2.03 on
// else if(!strncmp(ParameterString.c_str(),"I23S",4))
// jpsi.i23s = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IFPS",4))
cainpu.ifps = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"ITIMSHR",7))
casshwr.itimshr = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IRUNAEM",7))
capar1.irunaem = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IRUNA",5))
capar1.iruna = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IQ2",3))
capar1.iq2 = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IPRO",4))
capar1.ipro = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"NFLAV",5))
caluco.nflav = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"INTER",5))
cainpu.inter = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IHFLA",5))
cahflav.ihfla = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IRPA",4))
cascol.irpa = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IRPB",4))
cascol.irpb = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IRPC",4))
cascol.irpc = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"ICCFM",5))
casshwr.iccfm = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IFINAL",6))
cainpu.ifinal = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IGLU",4))
cagluon.iglu = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"IRspl",5))
casprre.irspl = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"PT2CUT",6))
captcut.pt2cut[capar1.ipro-1] = atof(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"NCB",3))
integr.ncb = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"ACC1",4))
integr.acc1 = atof(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"ACC2",4))
integr.acc2 = atof(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"SCALFAF",7))
scalf.scalfaf = atof(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else if(!strncmp(ParameterString.c_str(),"SCALFA",6))
scalf.scalfa = atof(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
else accepted = 0;
return accepted;
}
| const char * gen::Cascade2Hadronizer::classname | ( | ) | const |
Definition at line 635 of file Cascade2Hadronizer.cc.
{
return "gen::Cascade2Hadronizer";
}
| bool gen::Cascade2Hadronizer::decay | ( | ) |
Definition at line 283 of file Cascade2Hadronizer.cc.
{
return true;
}
| bool gen::Cascade2Hadronizer::declareSpecialSettings | ( | const std::vector< std::string > | ) | [inline] |
Definition at line 39 of file Cascade2Hadronizer.h.
{ return true; }
| bool gen::Cascade2Hadronizer::declareStableParticles | ( | const std::vector< int > | ) |
| void gen::Cascade2Hadronizer::fillTmpStorage | ( | ) | [private] |
Definition at line 161 of file Cascade2Hadronizer.cc.
References i, and gen::pyjets_local.
Referenced by generatePartonsAndHadronize().
{
pyjets_local.n = pyjets.n;
pyjets_local.npad = pyjets.npad;
for(int ip=0; ip<pyjets_maxn; ip++){
for(int i=0; i<5; i++){
pyjets_local.k[i][ip] = pyjets.k[i][ip];
pyjets_local.p[i][ip] = pyjets.p[i][ip];
pyjets_local.v[i][ip] = pyjets.v[i][ip];
}
}
return ;
}
| void gen::Cascade2Hadronizer::finalizeEvent | ( | ) |
Definition at line 176 of file Cascade2Hadronizer.cc.
References gen::call_pylist(), gather_cfg::cout, debug, gen::BaseHadronizer::event(), gen::BaseHadronizer::eventInfo(), fConvertToPDG, fHepMCVerbosity, fMaxEventsToPrint, fPythiaListVerbosity, imposeProperTime(), pydat1, pyint1, and pypars.
{
HepMC::PdfInfo pdf;
//-- filling factorization "Q scale" now! pthat moved to binningValues()
if(event()->signal_process_id() <= 0) event()->set_signal_process_id(pypars.msti[0]);
if(event()->event_scale() <=0 ) event()->set_event_scale(pypars.pari[22]);
if(event()->alphaQED() <= 0) event()->set_alphaQED(pyint1.vint[56]);
if(event()->alphaQCD() <= 0) event()->set_alphaQCD(pyint1.vint[57]);
//-- get pdf info directly from Pythia6 and set it up into HepMC::GenEvent
//-- S. Mrenna: Prefer vint block
if(pdf.id1() <= 0) pdf.set_id1(pyint1.mint[14] == 21 ? 0 : pyint1.mint[14]);
if(pdf.id2() <= 0) pdf.set_id2(pyint1.mint[15] == 21 ? 0 : pyint1.mint[15]);
if(pdf.x1() <= 0) pdf.set_x1(pyint1.vint[40]);
if(pdf.x2() <= 0) pdf.set_x2(pyint1.vint[41]);
if(pdf.pdf1() <= 0) pdf.set_pdf1(pyint1.vint[38] / pyint1.vint[40]);
if(pdf.pdf2() <= 0) pdf.set_pdf2(pyint1.vint[39] / pyint1.vint[41]);
if(pdf.scalePDF() <= 0) pdf.set_scalePDF(pyint1.vint[50]);
event()->set_pdf_info(pdf) ;
if(debug) {
cout<<"standard Py6 event weight: pyint1.vint[96]: "<<pyint1.vint[96]<<endl;
cout<<"event weight returned by PYEVWT: 1./(pyint1.vint[98]): "<<1./(pyint1.vint[98])<<endl;
}
//-- this is "standard" Py6 event weight (corresponds to PYINT1/VINT(97))
// event()->weights().push_back(pyint1.vint[96]);
//-- this is event weight as 1./VINT(99) (PYINT1/VINT(99) is returned by the PYEVWT)
// event()->weights().push_back(1./(pyint1.vint[98]));
//-- all cascade events have weight = 1
event()->weights().push_back(1.);
//-- now create the GenEventInfo product from the GenEvent and fill the missing pieces
eventInfo().reset(new GenEventInfoProduct(event().get()));
//-- in Pythia6 pthat is used to subdivide samples into different bins
//-- in LHE mode the binning is done by the external ME generator
//-- which is likely not pthat, so only filling it for Py6 internal mode
eventInfo()->setBinningValues(vector<double>(1,pypars.pari[16]));
//-- here we treat long-lived particles
if (pydat1.mstj[21]==3 || pydat1.mstj[21]==4) imposeProperTime();
//-- convert particle IDs Py6 -> PDG (if requested)
if(fConvertToPDG){
for(HepMC::GenEvent::particle_iterator part = event()->particles_begin();
part != event()->particles_end(); ++part){
(*part)->set_pdg_id(HepPID::translatePythiatoPDT((*part)->pdg_id()));
}
}
//-- service printouts, if requested
if(fMaxEventsToPrint > 0){
fMaxEventsToPrint--;
if(fPythiaListVerbosity) call_pylist(fPythiaListVerbosity);
if(fHepMCVerbosity){
cout <<"Event process = "<<pypars.msti[0]<<endl<<"----------------------"<<endl;
event()->print();
}
}
return;
}
| void gen::Cascade2Hadronizer::flushTmpStorage | ( | ) | [private] |
Definition at line 146 of file Cascade2Hadronizer.cc.
References i, and gen::pyjets_local.
Referenced by Cascade2Hadronizer(), and generatePartonsAndHadronize().
{
pyjets_local.n = 0 ;
pyjets_local.npad = 0 ;
for(int ip=0; ip<pyjets_maxn; ip++){
for(int i=0; i<5; i++){
pyjets_local.k[i][ip] = 0;
pyjets_local.p[i][ip] = 0.;
pyjets_local.v[i][ip] = 0.;
}
}
return;
}
| bool gen::Cascade2Hadronizer::generatePartonsAndHadronize | ( | ) |
Definition at line 254 of file Cascade2Hadronizer.cc.
References call_event(), gen::BaseHadronizer::event(), fillTmpStorage(), flushTmpStorage(), fPy6Service, gen::FortranCallback::getInstance(), hepevtio, pyint1, and gen::FortranCallback::resetIterationsPerEvent().
{
//-- grab Py6 instance
Pythia6Service::InstanceWrapper guard(fPy6Service);
FortranCallback::getInstance()->resetIterationsPerEvent();
//-- skip event if py6 considers it bad
if(pyint1.mint[50] != 0 ){
event().reset();
return false;
}
//-- generation of the event with CASCADE
call_event();
//-- pythia pyhepc routine converts common PYJETS in common HEPEVT
call_pyhepc(1);
//-- delete the created event from memory
event().reset(hepevtio.read_next_event());
//-- this is to deal with post-gen tools & residualDecay() that may reuse PYJETS
flushTmpStorage();
fillTmpStorage();
return true;
}
| bool gen::Cascade2Hadronizer::hadronize | ( | ) |
Definition at line 250 of file Cascade2Hadronizer.cc.
{
return false;
}
| void gen::Cascade2Hadronizer::imposeProperTime | ( | ) | [private] |
Definition at line 547 of file Cascade2Hadronizer.cc.
References abs, gen::FortranInstance::call(), gen::BaseHadronizer::event(), fPy6Service, funct::log(), gen::pycomp_(), pydat1, gen::pyr_(), mathSSE::sqrt(), lumiQTWidget::t, vbegin, vend, x, detailsBasic3DVector::y, and z.
Referenced by finalizeEvent().
{
// this is practically a copy/paste of the original code by J.Alcaraz,
// taken directly from PythiaSource
int dumm=0;
HepMC::GenEvent::vertex_const_iterator vbegin = event()->vertices_begin();
HepMC::GenEvent::vertex_const_iterator vend = event()->vertices_end();
HepMC::GenEvent::vertex_const_iterator vitr = vbegin;
for(; vitr != vend; ++vitr){
HepMC::GenVertex::particle_iterator pbegin = (*vitr)->particles_begin(HepMC::children);
HepMC::GenVertex::particle_iterator pend = (*vitr)->particles_end(HepMC::children);
HepMC::GenVertex::particle_iterator pitr = pbegin;
for(; pitr != pend; ++pitr){
if((*pitr)->end_vertex()) continue;
if((*pitr)->status()!=1) continue;
int pdgcode= abs((*pitr)->pdg_id());
// Do nothing if the particle is not expected to decay
if ( pydat3.mdcy[0][pycomp_(pdgcode)-1] !=1 ) continue;
double ctau = pydat2.pmas[3][pycomp_(pdgcode)-1];
HepMC::FourVector mom = (*pitr)->momentum();
HepMC::FourVector vin = (*vitr)->position();
double x = 0.;
double y = 0.;
double z = 0.;
double t = 0.;
bool decayInRange = false;
while(!decayInRange){
double unif_rand = fPy6Service->call(pyr_, &dumm);
// Value of 0 is excluded, so following line is OK
double proper_length = - ctau * log(unif_rand);
double factor = proper_length/mom.m();
x = vin.x() + factor * mom.px();
y = vin.y() + factor * mom.py();
z = vin.z() + factor * mom.pz();
t = vin.t() + factor * mom.e();
// Decay must be happen outside a cylindrical region
if(pydat1.mstj[21]==4){
if (sqrt(x*x+y*y)>pydat1.parj[72] || fabs(z)>pydat1.parj[73]) decayInRange = true;
// Decay must be happen outside a given sphere
}
else if (pydat1.mstj[21]==3){
if (sqrt(x*x+y*y+z*z)>pydat1.parj[71]) decayInRange = true;
}
// Decay is always OK otherwise
else {
decayInRange = true;
}
}
HepMC::GenVertex* vdec = new HepMC::GenVertex(HepMC::FourVector(x,y,z,t));
event()->add_vertex(vdec);
vdec->add_particle_in((*pitr));
}
}
return;
}
| bool gen::Cascade2Hadronizer::initializeForExternalPartons | ( | ) |
Definition at line 452 of file Cascade2Hadronizer.cc.
{
return false;
}
| bool gen::Cascade2Hadronizer::initializeForInternalPartons | ( | ) |
Definition at line 456 of file Cascade2Hadronizer.cc.
References cainpu, call_caend(), call_cascade(), call_cascha(), call_casini(), cascadePrintParameters(), cascadeReadParameters(), gen::Pythia6Service::closeSLHA(), edm::errors::Configuration, Exception, fComEnergy, fParameters, fPy6Service, edm::ParameterSet::getParameter(), i, pythia6PrintParameters(), gen::Pythia6Service::setCSAParams(), gen::Pythia6Service::setGeneralParams(), gen::Pythia6Service::setPYUPDAParams(), and gen::Pythia6Service::setSLHAParams().
{
//-- grab Py6 instance
Pythia6Service::InstanceWrapper guard(fPy6Service);
//-- change standard parameters of JETSET/PYTHIA - replace call_pytcha()
fPy6Service->setGeneralParams();
fPy6Service->setCSAParams();
fPy6Service->setSLHAParams();
fPy6Service->setPYUPDAParams(false);
pythia6PrintParameters();
//-- mstu(8) is set to NMXHEP in this dummy call (version >=6.404)
call_pyhepc(1);
//-- initialise random number generator: has been changed to be CMSSW compliant
//-- dcasrn overloaded: call to rluxgo not needed anymore
//-- call_rluxgo(4,314159265,0,0);
//-- initialise CASCADE parameters (default values)
call_casini();
//-- Read the parameters and pass them to the common blocks
//-- call_steer();
//-- initialise CASCADE parameters (user values)
//-- retrieve all the different sets
vector<string> AllSets = fParameters.getParameter<vector<string> >("parameterSets");
//-- loop over the different sets
for(unsigned i=0; i<AllSets.size();++i) {
string Set = AllSets[i];
vector<string> Para_Set = fParameters.getParameter<vector<string> >(Set);
//-- loop over all the parameters and stop in case of mistake
for(vector<string>::const_iterator itPara = Para_Set.begin(); itPara != Para_Set.end(); ++itPara) {
if(!cascadeReadParameters(*itPara)) {
throw edm::Exception(edm::errors::Configuration,"CascadeError")
<<" Cascade did not accept the following parameter: \""<<*itPara<<"\""<<endl;
}
} //-- end loop over all the parameters
} //-- end loop over the different sets
cainpu.plepin = -fComEnergy/2;
cainpu.ppin = fComEnergy/2;
cascadePrintParameters();
//-- change standard parameters of CASCADE
call_cascha();
//-- change standard parameters of JETSET/PYTHIA
//-- call_pytcha();
//-- set up for running CASCADE (integration of the cross-section)
call_cascade();
//-- print cross-section result from integration
call_caend(1);
fPy6Service->setPYUPDAParams(true);
fPy6Service->closeSLHA();
return true;
}
| void gen::Cascade2Hadronizer::pythia6PrintParameters | ( | ) |
Definition at line 804 of file Cascade2Hadronizer.cc.
References gather_cfg::cout, and pydat1.
Referenced by initializeForInternalPartons().
{
cout<<""<<endl;
cout<<"----------------------------------"<<endl;
cout<<"---- Pythia6 Parameters ----"<<endl;
cout<<"----------------------------------"<<endl;
cout<<""<<endl;
cout<<"charm mass: "<<pydat2.pmas[0][3]<<" GeV (default value = 1.5 GeV)"<<endl;
cout<<"bottom mass: "<<pydat2.pmas[0][4]<<" GeV (default value = 4.8 GeV)"<<endl;
cout<<"Higgs mass: "<<pydat2.pmas[0][24]<<" GeV (default value = 115 GeV)"<<endl;
cout<<"lambda QCD: "<<pydat1.paru[111]<<" GeV (default value = 0.25 GeV)"<<endl;
cout<<"alphas behaviour: "<<pydat1.mstu[110]<<" (default value = 1)"<<endl;
cout<<"nr of flavours wrt lambda QCD: "<<pydat1.mstu[111]<<" (default value = 5)"<<endl;
cout<<"min nr of flavours for alphas: "<<pydat1.mstu[112]<<" (default value = 3)"<<endl;
cout<<"max nr of flavours for alphas: "<<pydat1.mstu[113]<<" (default value = 5)"<<endl;
cout<<"maximum angle settings: "<<pydat1.mstj[47]<<" (default value = 0)"<<endl;
cout<<""<<endl;
}
| bool gen::Cascade2Hadronizer::readSettings | ( | int | key | ) |
Definition at line 439 of file Cascade2Hadronizer.cc.
References fPy6Service, gen::Pythia6Service::setCSAParams(), gen::Pythia6Service::setGeneralParams(), and gen::Pythia6Service::setSLHAParams().
{
//-- grab Py6 instance
Pythia6Service::InstanceWrapper guard(fPy6Service);
fPy6Service->setGeneralParams();
if (key == 0) fPy6Service->setCSAParams();
fPy6Service->setSLHAParams();
return true;
}
| bool gen::Cascade2Hadronizer::residualDecay | ( | ) |
Definition at line 287 of file Cascade2Hadronizer.cc.
References gen::BaseHadronizer::event(), fPy6Service, configurableAnalysis::GenParticle, i, dbtoconf::parent, gen::pycomp_(), gen::pydecy_(), gen::pyjets_local, and ntuplemaker::status.
{
//-- grab Py6 instance
Pythia6Service::InstanceWrapper guard(fPy6Service);
int NPartsBeforeDecays = pyjets_local.n ;
int NPartsAfterDecays = event().get()->particles_size();
int barcode = NPartsAfterDecays;
// JVY: well, in principle, it's not a 100% fair to go up to NPartsBeforeDecays,
// because Photos will attach gamma's to existing vertices, i.e. in the middle
// of the event rather than at the end; but this will only shift pointers down,
// so we'll be going again over a few "original" particle...
// in the alternative, we may go all the way up to the beginning of the event
// and re-check if anything remains to decay, that's fine even if it'll take
// some extra CPU...
for(int ipart=NPartsAfterDecays; ipart>NPartsBeforeDecays; ipart-- ){
HepMC::GenParticle* part = event().get()->barcode_to_particle( ipart );
int status = part->status();
if ( status != 1 ) continue; // check only "stable" particles,
// as some undecayed may still be marked as such
int pdgid = part->pdg_id();
int py6ptr = pycomp_( pdgid );
if(pydat3.mdcy[0][py6ptr-1] != 1) continue; // particle is not expected to decay
int py6id = HepPID::translatePDTtoPythia(pdgid);
// first, will need to zero out, then fill up PYJETS
// I better do it directly (by hands) rather than via py1ent
// - to avoid (additional) mass smearing
if(part->momentum().t() <= part->generated_mass()){
continue ; // e == m -> 0-momentum, nothing to decay...
}
else{
pyjets.n = 0;
for(int i=0; i<5; i++){
pyjets.k[i][0] = 0;
pyjets.p[i][0] = 0.;
pyjets.v[i][0] = 0.;
}
pyjets.k[0][0] = 1;
pyjets.k[1][0] = py6id;
pyjets.p[4][0] = part->generated_mass();
pyjets.p[3][0] = part->momentum().t();
pyjets.p[0][0] = part->momentum().x();
pyjets.p[1][0] = part->momentum().y();
pyjets.p[2][0] = part->momentum().z();
HepMC::GenVertex* prod_vtx = part->production_vertex();
if ( !prod_vtx ) continue; // in principle, should never happen but...
pyjets.v[0][0] = prod_vtx->position().x();
pyjets.v[1][0] = prod_vtx->position().y();
pyjets.v[2][0] = prod_vtx->position().z();
pyjets.v[3][0] = prod_vtx->position().t();
pyjets.v[4][0] = 0.;
pyjets.n = 1;
pyjets.npad = pyjets_local.npad;
}
//-- now call Py6 decay routine
int parent = 1; // since we always pass to Py6 a single particle
pydecy_(parent);
//-- now attach decay products to mother
for(int iprt1=1; iprt1<pyjets.n; iprt1++){
part->set_status( 2 );
HepMC::GenVertex* DecVtx = new HepMC::GenVertex( HepMC::FourVector(pyjets.v[0][iprt1],
pyjets.v[1][iprt1],
pyjets.v[2][iprt1],
pyjets.v[3][iprt1]) );
DecVtx->add_particle_in( part ); // this will cleanup end_vertex if exists, replace with the new one
// I presume (vtx) barcode will be given automatically
HepMC::FourVector pmom(pyjets.p[0][iprt1],pyjets.p[1][iprt1],
pyjets.p[2][iprt1],pyjets.p[3][iprt1] );
int dstatus = 0;
if(pyjets.k[0][iprt1] >= 1 && pyjets.k[0][iprt1] <= 10) dstatus = 1;
else if(pyjets.k[0][iprt1] >= 11 && pyjets.k[0][iprt1] <= 20) dstatus = 2;
else if(pyjets.k[0][iprt1] >= 21 && pyjets.k[0][iprt1] <= 30) dstatus = 3;
else if(pyjets.k[0][iprt1] >= 31 && pyjets.k[0][iprt1] <= 100) dstatus = pyjets.k[0][iprt1];
HepMC::GenParticle* daughter = new HepMC::GenParticle(pmom,
HepPID::translatePythiatoPDT(pyjets.k[1][iprt1]),
dstatus);
barcode++;
daughter->suggest_barcode( barcode );
DecVtx->add_particle_out( daughter );
int iprt2;
for(iprt2=iprt1+1; iprt2<pyjets.n; iprt2++){ // the pointer is shifted by -1, c++ style
if(pyjets.k[2][iprt2] != parent){
parent = pyjets.k[2][iprt2];
break; // another parent particle; reset & break the loop
}
HepMC::FourVector pmomN(pyjets.p[0][iprt2],pyjets.p[1][iprt2],
pyjets.p[2][iprt2],pyjets.p[3][iprt2] );
dstatus = 0;
if(pyjets.k[0][iprt2] >= 1 && pyjets.k[0][iprt2] <= 10) dstatus = 1;
else if(pyjets.k[0][iprt2] >= 11 && pyjets.k[0][iprt2] <= 20) dstatus = 2;
else if(pyjets.k[0][iprt2] >= 21 && pyjets.k[0][iprt2] <= 30) dstatus = 3;
else if(pyjets.k[0][iprt2] >= 31 && pyjets.k[0][iprt2] <= 100) dstatus = pyjets.k[0][iprt2];
HepMC::GenParticle* daughterN = new HepMC::GenParticle(pmomN,
HepPID::translatePythiatoPDT( pyjets.k[1][iprt2] ),
dstatus);
barcode++;
daughterN->suggest_barcode( barcode );
DecVtx->add_particle_out( daughterN );
} //-- end iprt2 loop
iprt1 = iprt2-1; // reset counter such that it doesn't go over the same child more than once
// don't forget to offset back into c++ counting, as it's already +1 forward
event().get()->add_vertex( DecVtx );
} //-- end iprt1 loop
} //-- end loop over decay products
// now restore the very original Py6 event record
if(pyjets_local.n != pyjets.n){
// restore pyjets to its state as it was before external decays -
// might have been jammed by action above or by py1ent calls in EvtGen
pyjets.n = pyjets_local.n;
pyjets.npad = pyjets_local.npad;
for(int ip=0; ip<pyjets_local.n; ip++){
for(int i=0; i<5; i++){
pyjets.k[i][ip] = pyjets_local.k[i][ip];
pyjets.p[i][ip] = pyjets_local.p[i][ip];
pyjets.v[i][ip] = pyjets_local.v[i][ip];
}
}
}
return true;
}
| void gen::Cascade2Hadronizer::statistics | ( | ) |
Definition at line 614 of file Cascade2Hadronizer.cc.
References caeffic, call_caend(), fextCrossSection, fextCrossSectionError, gen::BaseHadronizer::runInfo(), GenRunInfoProduct::setExternalXSecLO(), GenRunInfoProduct::setExternalXSecNLO(), and GenRunInfoProduct::setInternalXSec().
{
runInfo().setExternalXSecLO(GenRunInfoProduct::XSec(fextCrossSection,fextCrossSectionError));
runInfo().setExternalXSecNLO(GenRunInfoProduct::XSec(-1.,-1.));
if(!runInfo().internalXSec()){
double sigma_CMS = 0.001 * caeffic.avgi;
double error_CMS = 0.001 * caeffic.sd;
runInfo().setInternalXSec(GenRunInfoProduct::XSec(sigma_CMS,error_CMS));
}
//-- print out generated event summary
call_caend(2);
//-- write out some information from Pythia to the screen
call_pystat(1);
return;
}
double gen::Cascade2Hadronizer::fComEnergy [private] |
Definition at line 66 of file Cascade2Hadronizer.h.
Referenced by cascadePrintParameters(), and initializeForInternalPartons().
bool gen::Cascade2Hadronizer::fConvertToPDG [private] |
Definition at line 78 of file Cascade2Hadronizer.h.
Referenced by Cascade2Hadronizer(), and finalizeEvent().
bool gen::Cascade2Hadronizer::fDisplayPythiaBanner [private] |
Definition at line 75 of file Cascade2Hadronizer.h.
Referenced by Cascade2Hadronizer().
bool gen::Cascade2Hadronizer::fDisplayPythiaCards [private] |
Definition at line 76 of file Cascade2Hadronizer.h.
Referenced by Cascade2Hadronizer().
double gen::Cascade2Hadronizer::fextCrossSection [private] |
Definition at line 67 of file Cascade2Hadronizer.h.
Referenced by cascadePrintParameters(), and statistics().
double gen::Cascade2Hadronizer::fextCrossSectionError [private] |
Definition at line 68 of file Cascade2Hadronizer.h.
Referenced by cascadePrintParameters(), and statistics().
double gen::Cascade2Hadronizer::fFilterEfficiency [private] |
Definition at line 69 of file Cascade2Hadronizer.h.
Referenced by cascadePrintParameters().
CLHEP::RandFlat* gen::Cascade2Hadronizer::fFlat [private] |
Definition at line 64 of file Cascade2Hadronizer.h.
Referenced by Cascade2Hadronizer().
bool gen::Cascade2Hadronizer::fHepMCVerbosity [private] |
Definition at line 72 of file Cascade2Hadronizer.h.
Referenced by finalizeEvent().
unsigned int gen::Cascade2Hadronizer::fMaxEventsToPrint [private] |
Definition at line 71 of file Cascade2Hadronizer.h.
Referenced by finalizeEvent().
Definition at line 60 of file Cascade2Hadronizer.h.
Referenced by Cascade2Hadronizer(), and initializeForInternalPartons().
Definition at line 62 of file Cascade2Hadronizer.h.
Referenced by generatePartonsAndHadronize(), imposeProperTime(), initializeForInternalPartons(), readSettings(), residualDecay(), and ~Cascade2Hadronizer().
unsigned int gen::Cascade2Hadronizer::fPythiaListVerbosity [private] |
Definition at line 73 of file Cascade2Hadronizer.h.
Referenced by finalizeEvent().
1.7.3