#include <PomwigHadronizer.h>

Inheritance diagram for gen::PomwigHadronizer:

Public Member Functions
const char *	classname () const
bool	decay ()
bool	declareSpecialSettings (const std::vector< std::string >)
bool	declareStableParticles (const std::vector< int > &pdgIds)
void	finalizeEvent ()
bool	generatePartonsAndHadronize ()
bool	hadronize ()
bool	initializeForExternalPartons ()
bool	initializeForInternalPartons ()
	PomwigHadronizer (const edm::ParameterSet &params)
bool	readSettings (int)
bool	residualDecay ()
void	statistics ()
	~PomwigHadronizer ()
Private Member Functions
void	clear ()
bool	initializeDPDF ()
Private Attributes
double	comEnergy
HepMC::IO_HERWIG	conv
int	diffTopology
bool	doMPInteraction
bool	doPDGConvert
int	h1fit
int	hepmcVerbosity
int	herwigVerbosity
int	maxEventsToPrint
bool	needClear
int	numTrials
gen::ParameterCollector	parameters
bool	printCards
double	survivalProbability
bool	useJimmy

Detailed Description

Definition at line 14 of file PomwigHadronizer.h.

Constructor & Destructor Documentation

gen::PomwigHadronizer::PomwigHadronizer ( const edm::ParameterSet & params )

Definition at line 89 of file PomwigHadronizer.cc.

References doPDGConvert, edm::ParameterSet::exists(), and edm::ParameterSet::getParameter().

                                                                :
        BaseHadronizer(params),
        needClear(false),
        parameters(params.getParameter<edm::ParameterSet>("HerwigParameters")),
        herwigVerbosity(params.getUntrackedParameter<int>("herwigVerbosity", 0)),
        hepmcVerbosity(params.getUntrackedParameter<int>("hepmcVerbosity", 0)),
        maxEventsToPrint(params.getUntrackedParameter<int>("maxEventsToPrint", 0)),
        printCards(params.getUntrackedParameter<bool>("printCards", false)),
        comEnergy(params.getParameter<double>("comEnergy")),
        survivalProbability(params.getParameter<double>("survivalProbability")),
        diffTopology(params.getParameter<int>("diffTopology")),
        h1fit(params.getParameter<int>("h1fit")),
        useJimmy(params.getParameter<bool>("useJimmy")),
        doMPInteraction(params.getParameter<bool>("doMPInteraction")),
        numTrials(params.getUntrackedParameter<int>("numTrialsMPI", 100)),
        doPDGConvert(false)
{
        if ( params.exists( "doPDGConvert" ) )
           doPDGConvert = params.getParameter<bool>("doPDGConvert");
}

gen::PomwigHadronizer::~PomwigHadronizer ( )

Definition at line 110 of file PomwigHadronizer.cc.

References clear().

{
        clear();
}

Member Function Documentation

const char* gen::PomwigHadronizer::classname ( ) const [inline]

Definition at line 34 of file PomwigHadronizer.h.

{ return "PomwigHadronizer"; }

void gen::PomwigHadronizer::clear ( void ) [private]

Definition at line 115 of file PomwigHadronizer.cc.

References gen::FortranInstance::call(), jmefin, needClear, and useJimmy.

Referenced by readSettings(), and ~PomwigHadronizer().

{
        if (!needClear)
                return;

        // teminate elementary process
        call(hwefin);
        if (useJimmy)
                call(jmefin);

        needClear = false;
}

bool gen::PomwigHadronizer::decay ( )

Definition at line 348 of file PomwigHadronizer.cc.

References conv, doPDGConvert, gen::BaseHadronizer::event(), and wrapper.

{
        // hadron decays

        InstanceWrapper wrapper(this);  // safe guard

        hwdhad();       // unstable particle decays
        hwdhvy();       // heavy flavour decays
        hwmevt();       // soft underlying event

        hwufne();       // finalize event

        if (hwevnt.IERROR)
                return false;

        event().reset(new HepMC::GenEvent);
        if (!conv.fill_next_event(event().get()))
                throw cms::Exception("PomwigError")
                        << "HepMC Conversion problems in event." << std::endl;

        // do particle ID conversion Herwig->PDG, if requested
        if(doPDGConvert){
           for ( HepMC::GenEvent::particle_iterator part = event()->particles_begin(); part != event()->particles_end(); ++part) {
             if ((*part)->pdg_id() != HepPID::translateHerwigtoPDT((*part)->pdg_id()))
               (*part)->set_pdg_id(HepPID::translateHerwigtoPDT((*part)->pdg_id()));
           }
        }

        return true;
}

bool gen::PomwigHadronizer::declareSpecialSettings ( const std::vector< std::string > ) [inline]

Definition at line 25 of file PomwigHadronizer.h.

{ return true; }

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

Definition at line 279 of file PomwigHadronizer.cc.

{
        for(std::vector<int>::const_iterator iter = pdgIds.begin();
            iter != pdgIds.end(); ++iter)
                if (!markStable(*iter))
                        return false;
        return true;
}

void gen::PomwigHadronizer::finalizeEvent ( )

Definition at line 339 of file PomwigHadronizer.cc.

References gen::BaseHadronizer::event(), and lhef::LHEEvent::fixHepMCEventTimeOrdering().

{
        lhef::LHEEvent::fixHepMCEventTimeOrdering(event().get());

        event()->set_signal_process_id(hwproc.IPROC);

        event()->weights().push_back(hwevnt.EVWGT);
}

bool gen::PomwigHadronizer::generatePartonsAndHadronize ( )

Definition at line 304 of file PomwigHadronizer.cc.

References gen::Herwig6Instance::callWithTimeout(), doMPInteraction, error, gen::BaseHadronizer::event(), gen::hwwarn_(), useJimmy, and wrapper.

{
        // hard process generation, parton shower, hadron formation

        InstanceWrapper wrapper(this);  // safe guard

        event().reset();

        // call herwig routines to create HEPEVT

        hwuine();       // initialize event
        
        if (callWithTimeout(10, hwepro)) { // process event and PS
          // We hung for more than 10 seconds
          int error = 199;
          hwwarn_("HWHGUP", &error);
        }
        
        hwbgen();       // parton cascades

        // call jimmy ... only if event is not killed yet by HERWIG
        if (useJimmy && doMPInteraction && !hwevnt.IERROR && call_hwmsct()) 
          return false;
        
        hwdhob();       // heavy quark decays
        hwcfor();       // cluster formation
        hwcdec();       // cluster decays
        
        // if event *not* killed by HERWIG, return true
        if (!hwevnt.IERROR) return true;

        hwufne();       // finalize event       
        return false;
}

bool gen::PomwigHadronizer::hadronize ( )

Definition at line 299 of file PomwigHadronizer.cc.

{
   return false;
}

bool gen::PomwigHadronizer::initializeDPDF ( ) [private]

Definition at line 233 of file PomwigHadronizer.cc.

References printConversionInfo::aux, edm::errors::Configuration, diffTopology, Exception, python::connectstrParser::f2, h1fit, gen::qcd_1994(), and gen::qcd_2006().

Referenced by initializeForInternalPartons().

{
   // Initialize H1 pomeron/reggeon

   if(diffTopology == 3) return true;
 
   if((diffTopology != 0)&&(diffTopology != 1)&&(diffTopology != 2)) return false;

   int nstru = hwpram.NSTRU;
   int ifit = h1fit;
   if((nstru == 9)||(nstru == 10)){
         if((ifit <= 0)||(ifit >= 7)){
              throw edm::Exception(edm::errors::Configuration,"PomwigError")
              <<" Attempted to set non existant H1 1997 fit index. Has to be 1...6";
         }
         std::string aux((nstru == 9)?"Pomeron":"Reggeon"); 
         edm::LogVerbatim("") << "   H1 1997 pdf's: " << aux << "\n"
                              << "   IFIT = " << ifit;
         double xp = 0.1;
         double Q2 = 75.0;
         double xpq[13];
         qcd_1994(xp,Q2,xpq,ifit);
    } else if((nstru >= 12)&&(nstru <= 15)){
         bool isPom = (nstru == 12)||(nstru == 14);
         bool isFitA = (nstru == 12)||(nstru == 13);
         ifit = (isFitA)?1:2;
         std::string aux_0((isFitA)?"A":"B");
         std::string aux_1((isPom)?"Pomeron":"Reggeon"); 
         edm::LogVerbatim("") << "   H1 2006 Fit " << aux_0 << " " << aux_1 << "\n"
                              << "   IFIT = "<< ifit;
         double xp = 0.1;
         double Q2 = 75.0;
         double xpq[13];
         double f2[2];
         double fl[2];
         double c2[2];
         double cl[2];
         qcd_2006(xp,Q2,ifit,xpq,f2,fl,c2,cl);
    } else{
         throw edm::Exception(edm::errors::Configuration,"PomwigError")
               <<" Only running Pomeron H1 1997 (NSTRU=9), H1 2006 fit A (NSTRU=12) and H1 2006 fit B (NSTRU=14) or Reggeon H1 1997 (NSTRU=10), H1 2006 fit A (NSTRU=13) and H1 2006 fit B (NSTRU=15)";
    }
 
    return true;
}

bool gen::PomwigHadronizer::initializeForExternalPartons ( )

Definition at line 128 of file PomwigHadronizer.cc.

{
   return false;
}

bool gen::PomwigHadronizer::initializeForInternalPartons ( )

Definition at line 219 of file PomwigHadronizer.cc.

References gen::FortranInstance::call(), and initializeDPDF().

{
 
   call(hwuinc);

   hwusta("PI0     ",1);
 
   if(!initializeDPDF()) return false;

   call(hweini);

   return true;
}

bool gen::PomwigHadronizer::readSettings ( int )

Definition at line 133 of file PomwigHadronizer.cc.

References gen::ParameterCollector::begin(), gen::FortranInstance::call(), clear(), comEnergy, edm::errors::Configuration, diffTopology, gen::ParameterCollector::end(), Exception, gen::Herwig6Instance::give(), herwigVerbosity, hwdspn, hwprch, i, geometryCSVtoXML::line, maxEventsToPrint, needClear, and parameters.

{

   clear();

   edm::LogVerbatim("") << "----------------------------------------------\n"
                        << "Initializing PomwigHadronizer\n"
                        << "----------------------------------------------\n";

   // Call hwudat to set up HERWIG block data
   //hwudat();
 
   // Setting basic parameters ...
   hwproc.PBEAM1 = comEnergy/2.;
   hwproc.PBEAM2 = comEnergy/2.;
   // Choose beam particles for POMWIG depending on topology
   switch (diffTopology){
        case 0: //DPE
                hwbmch.PART1[0]  = 'E';
                hwbmch.PART1[1]  = '-';
                hwbmch.PART2[0]  = 'E';
                hwbmch.PART2[1]  = '-';
                break;
        case 1: //SD survive PART1
                hwbmch.PART1[0]  = 'E';
                hwbmch.PART1[1]  = '-';
                hwbmch.PART2[0]  = 'P';
                hwbmch.PART2[1]  = ' ';
                break;
        case 2: //SD survive PART2
                hwbmch.PART1[0]  = 'P';
                hwbmch.PART1[1]  = ' ';
                hwbmch.PART2[0]  = 'E';
                hwbmch.PART2[1]  = '-';
                break;
        case 3: //Non diffractive
                hwbmch.PART1[0]  = 'P';
                hwbmch.PART1[1]  = ' ';
                hwbmch.PART2[0]  = 'P';
                hwbmch.PART2[1]  = ' ';
                break;
        default:
                throw edm::Exception(edm::errors::Configuration,"PomwigError")
          <<" Invalid Diff. Topology. Must be DPE(diffTopology = 0), SD particle 1 (diffTopology = 1), SD particle 2 (diffTopology = 2) and Non diffractive (diffTopology = 3)";
                break;
   }
   for(int i=2;i<8;++i){
    hwbmch.PART1[i]  = ' ';
    hwbmch.PART2[i]  = ' ';}

   // initialize other common blocks ...
   call(hwigin);

   hwevnt.MAXER = 100000000;    // O(inf)
   hwpram.LWSUD = 0;            // don't write Sudakov form factors
   hwdspn.LWDEC = 0;            // don't write three/four body decays
                                // (no fort.77 and fort.88 ...)a

   std::memset(hwprch.AUTPDF, ' ', sizeof hwprch.AUTPDF);
   for(unsigned int i = 0; i < 2; i++) {
        hwpram.MODPDF[i] = -111;
        std::memcpy(hwprch.AUTPDF[i], "HWLHAPDF", 8);
   }

   hwevnt.MAXPR = maxEventsToPrint;
   hwpram.IPRINT = herwigVerbosity;

   edm::LogVerbatim("") << "------------------------------------\n"
                        << "Reading HERWIG parameters\n"
                        << "------------------------------------\n";
    
   for(gen::ParameterCollector::const_iterator line = parameters.begin();
                                               line != parameters.end(); ++line) {
      edm::LogVerbatim("") << "   " << *line;
      if (!give(*line))
          throw edm::Exception(edm::errors::Configuration)
                << "Herwig 6 did not accept the following: \""
                << *line << "\"." << std::endl;
   }

   needClear = true;

   return true;

}

bool gen::PomwigHadronizer::residualDecay ( )

Definition at line 379 of file PomwigHadronizer.cc.

{
        return true;
}

void gen::PomwigHadronizer::statistics ( )

Definition at line 288 of file PomwigHadronizer.cc.

References gen::BaseHadronizer::runInfo(), GenRunInfoProduct::setInternalXSec(), and survivalProbability.

{
        double RNWGT = 1. / hwevnt.NWGTS;
        double AVWGT = hwevnt.WGTSUM * RNWGT;

        double xsec = 1.0e3 * AVWGT;
        xsec = survivalProbability*xsec;

        runInfo().setInternalXSec(xsec);
}