#include <Cascade2Hadronizer.h>

Inheritance diagram for gen::Cascade2Hadronizer:

Public Member Functions
	Cascade2Hadronizer (edm::ParameterSet const &ps)

void	cascadePrintParameters ()

bool	cascadeReadParameters (const std::string &ParameterString)

const char *	classname () const

bool	decay ()

bool	declareSpecialSettings (const std::vector< std::string > &)

bool	declareStableParticles (const std::vector< int > &)

void	finalizeEvent ()

bool	generatePartonsAndHadronize ()

bool	hadronize ()

bool	initializeForExternalPartons ()

bool	initializeForInternalPartons ()

void	pythia6PrintParameters ()

bool	readSettings (int)

bool	residualDecay ()

void	statistics ()

	~Cascade2Hadronizer ()

Public Member Functions inherited from gen::BaseHadronizer
	BaseHadronizer (edm::ParameterSet const &ps)

edm::Event &	getEDMEvent () const

HepMC::GenEvent *	getGenEvent ()

GenEventInfoProduct *	getGenEventInfo ()

GenRunInfoProduct &	getGenRunInfo ()

const boost::shared_ptr < lhef::LHERunInfo > &	getLHERunInfo () const

void	resetEvent (HepMC::GenEvent *event)

void	resetEventInfo (GenEventInfoProduct *eventInfo)

virtual bool	select (HepMC::GenEvent *) const

void	setEDMEvent (edm::Event &event)

void	setLHEEvent (lhef::LHEEvent *event)

void	setLHERunInfo (lhef::LHERunInfo *runInfo)

void	setRandomEngine (CLHEP::HepRandomEngine *v)

std::vector< std::string > const &	sharedResources () const

	~BaseHadronizer ()

Private Member Functions
void	fillTmpStorage ()

void	flushTmpStorage ()

void	imposeProperTime ()

Private Attributes
double	fComEnergy

bool	fConvertToPDG

bool	fDisplayPythiaBanner

bool	fDisplayPythiaCards

double	fextCrossSection

double	fextCrossSectionError

double	fFilterEfficiency

CLHEP::RandFlat *	fFlat

bool	fHepMCVerbosity

unsigned int	fMaxEventsToPrint

edm::ParameterSet	fParameters

Pythia6Service *	fPy6Service

unsigned int	fPythiaListVerbosity

Additional Inherited Members
Protected Member Functions inherited from gen::BaseHadronizer
std::auto_ptr< HepMC::GenEvent > &	event ()

std::auto_ptr < GenEventInfoProduct > &	eventInfo ()

lhef::LHEEvent *	lheEvent ()

lhef::LHERunInfo *	lheRunInfo ()

GenRunInfoProduct &	runInfo ()

Detailed Description

Definition at line 23 of file Cascade2Hadronizer.h.

Constructor & Destructor Documentation

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, edm::hlt::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;
   }

Member Function Documentation

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.

39 { return true; }

bool gen::Cascade2Hadronizer::declareStableParticles ( const std::vector< int > & _pdg )

Definition at line 529 of file Cascade2Hadronizer.cc.

References gen::call_pygive(), i, and gen::pycomp_().

                                                                         {
     vector<int> pdg = _pdg; 
     for(size_t i=0; i<pdg.size(); i++){
       int PyID = HepPID::translatePDTtoPythia(pdg[i]);
       int pyCode = pycomp_( PyID );
 
       if(pyCode > 0){
     ostringstream pyCard ;
     pyCard << "MDCY(" << pyCode << ",1)=0";
     //-- cout << "pdg= " << pdg[i] << " " << pyCard.str() << endl; 
 
     call_pygive(pyCard.str());
       }
     }
       
     return true;
   }

void gen::Cascade2Hadronizer::fillTmpStorage ( )

private

Definition at line 161 of file Cascade2Hadronizer.cc.

References i, gen::pyjets_local, and reco::return().

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 funct::abs(), gen::FortranInstance::call(), gen::BaseHadronizer::event(), fPy6Service, create_public_lumi_plots::log, gen::pycomp_(), pydat1, gen::pyr_(), mathSSE::sqrt(), lumiQTWidget::t, vbegin, vend, x, detailsBasic3DVector::y, and detailsBasic3DVector::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, edm::hlt::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;
   }