diff --git a/Configuration/Applications/python/ConfigBuilder.py b/Configuration/Applications/python/ConfigBuilder.py index da6611baf976b..cf9e8535eff95 100644 --- a/Configuration/Applications/python/ConfigBuilder.py +++ b/Configuration/Applications/python/ConfigBuilder.py @@ -1379,6 +1379,8 @@ def prepare_GEN(self, sequence = None): if self._options.scenario == 'HeavyIons': if self._options.pileup=='HiMixGEN': self.loadAndRemember("Configuration/StandardSequences/GeneratorMix_cff") + elif self._options.pileup=='HiMixEmbGEN': + self.loadAndRemember("Configuration/StandardSequences/GeneratorEmbMix_cff") else: self.loadAndRemember("Configuration/StandardSequences/GeneratorHI_cff") @@ -2175,6 +2177,9 @@ def prepare(self, doChecking = False): self.pythonCfgCode +='\n' self.pythonCfgCode +=dumpPython(self.process,object) + if self._options.pileup=='HiMixEmbGEN': + self.pythonCfgCode += "\nprocess.generator.embeddingMode=cms.bool(True)\n" + # dump all paths self.pythonCfgCode += "\n# Path and EndPath definitions\n" for path in self.process.paths: @@ -2239,7 +2244,7 @@ def prepare(self, doChecking = False): MassReplaceInputTag(self.process, new="rawDataMapperByLabel", old="rawDataCollector") # special treatment in case of production filter sequence 2/2 - if self.productionFilterSequence: + if self.productionFilterSequence and not (self._options.pileup=='HiMixEmbGEN'): self.pythonCfgCode +='# filter all path with the production filter sequence\n' self.pythonCfgCode +='for path in process.paths:\n' if len(self.conditionalPaths): diff --git a/Configuration/Generator/python/Pyquen2013Settings_cff.py b/Configuration/Generator/python/Pyquen2013Settings_cff.py index 8ff1ceabeafba..21d44f40a5e7c 100644 --- a/Configuration/Generator/python/Pyquen2013Settings_cff.py +++ b/Configuration/Generator/python/Pyquen2013Settings_cff.py @@ -28,7 +28,6 @@ pyquenParameters = cms.PSet(doIsospin = cms.bool(True), angularSpectrumSelector = cms.int32(0), ## angular emitted gluon spectrum : embeddingMode = cms.bool(False), - backgroundLabel = cms.InputTag("generator","unsmeared") ## ineffective in no mixing ) hydjetParameters = cms.PSet(sigmaInelNN = cms.double(64), diff --git a/Configuration/Generator/python/Pyquen2015Settings_cff.py b/Configuration/Generator/python/Pyquen2015Settings_cff.py index 87c4f6c9684ba..230c19ab5f98e 100644 --- a/Configuration/Generator/python/Pyquen2015Settings_cff.py +++ b/Configuration/Generator/python/Pyquen2015Settings_cff.py @@ -28,7 +28,6 @@ pyquenParameters = cms.PSet(doIsospin = cms.bool(True), angularSpectrumSelector = cms.int32(0), ## angular emitted gluon spectrum embeddingMode = cms.bool(False), - backgroundLabel = cms.InputTag("generator","unsmeared") ## ineffective in no mixing ) hydjetParameters = cms.PSet(sigmaInelNN = cms.double(70), diff --git a/Configuration/Generator/python/PyquenDefaultSettings_cff.py b/Configuration/Generator/python/PyquenDefaultSettings_cff.py index 96e1e389cee93..5e8d72738a521 100644 --- a/Configuration/Generator/python/PyquenDefaultSettings_cff.py +++ b/Configuration/Generator/python/PyquenDefaultSettings_cff.py @@ -26,7 +26,6 @@ pyquenParameters = cms.PSet(doIsospin = cms.bool(True), angularSpectrumSelector = cms.int32(0), ## angular emitted gluon spectrum : embeddingMode = cms.bool(False), - backgroundLabel = cms.InputTag("generator","unsmeared") ## ineffective in no mixing ) hydjetParameters = cms.PSet(sigmaInelNN = cms.double(58), diff --git a/Configuration/StandardSequences/python/GeneratorEmbMix_cff.py b/Configuration/StandardSequences/python/GeneratorEmbMix_cff.py new file mode 100644 index 0000000000000..40227e2c86f7c --- /dev/null +++ b/Configuration/StandardSequences/python/GeneratorEmbMix_cff.py @@ -0,0 +1,4 @@ +import FWCore.ParameterSet.Config as cms +from Configuration.StandardSequences.GeneratorMix_cff import * + +pgen = cms.Sequence(cms.SequencePlaceholder("mix")+cms.SequencePlaceholder("generator")+cms.SequencePlaceholder("randomEngineStateProducer")+VertexSmearing+GenSmeared+genParticles+hiGenJets) diff --git a/Configuration/StandardSequences/python/Mixing.py b/Configuration/StandardSequences/python/Mixing.py index 20f4b07f4ffda..47bafd2d16e68 100644 --- a/Configuration/StandardSequences/python/Mixing.py +++ b/Configuration/StandardSequences/python/Mixing.py @@ -53,6 +53,7 @@ def addMixingScenario(label,dict): addMixingScenario("E8TeV_AVE_10_BX_25ns_300ns_spread",{'file':'SimGeneral.MixingModule.mix_E8TeV_AVE_10_BX_25ns_300ns_spread_cfi'}) addMixingScenario("HiMix",{'file': 'SimGeneral.MixingModule.HiMix_cff'}) addMixingScenario("HiMixGEN",{'file': 'SimGeneral.MixingModule.HiMixGEN_cff'}) +addMixingScenario("HiMixEmbGEN",{'file': 'SimGeneral.MixingModule.HiMixEmbGEN_cff'}) addMixingScenario("HiMixNoPU",{'file': 'SimGeneral.MixingModule.HiMixNoPU_cff'}) addMixingScenario("HighLumiPileUp",{'file': 'SimGeneral.MixingModule.mixHighLumPU_cfi'}) addMixingScenario("InitialLumiPileUp",{'file': 'SimGeneral.MixingModule.mixInitialLumPU_cfi'}) diff --git a/GeneratorInterface/Core/interface/GeneratorFilter.h b/GeneratorInterface/Core/interface/GeneratorFilter.h index 499023e829847..de8816869f0cd 100644 --- a/GeneratorInterface/Core/interface/GeneratorFilter.h +++ b/GeneratorInterface/Core/interface/GeneratorFilter.h @@ -61,10 +61,13 @@ namespace edm { void preallocThreads(unsigned int iThreads) override; private: + // two-phase construction to allow specialization of the constructor + void init(ParameterSet const& ps); + Hadronizer hadronizer_; //gen::ExternalDecayDriver* decayer_; - Decayer* decayer_; - unsigned int nEventsInLumiBlock_; + Decayer* decayer_ = nullptr; + unsigned int nEventsInLumiBlock_ = 0; unsigned int nThreads_{1}; bool initialized_ = false; }; @@ -74,8 +77,12 @@ namespace edm { // Implementation template - GeneratorFilter::GeneratorFilter(ParameterSet const& ps) - : EDFilter(), hadronizer_(ps), decayer_(nullptr), nEventsInLumiBlock_(0) { + GeneratorFilter::GeneratorFilter(ParameterSet const& ps) : hadronizer_(ps) { + init(ps); + } + + template + void GeneratorFilter::init(ParameterSet const& ps) { // TODO: // Put the list of types produced by the filters here. // The current design calls for: diff --git a/GeneratorInterface/HydjetInterface/BuildFile.xml b/GeneratorInterface/HydjetInterface/BuildFile.xml index 0fae129434faf..2726fc4683476 100644 --- a/GeneratorInterface/HydjetInterface/BuildFile.xml +++ b/GeneratorInterface/HydjetInterface/BuildFile.xml @@ -1,12 +1,12 @@ - - - - - - - - - + + + + + + + + + - + diff --git a/GeneratorInterface/HydjetInterface/interface/HydjetGeneratorFilter.h b/GeneratorInterface/HydjetInterface/interface/HydjetGeneratorFilter.h index 0840330161612..f09eea7889dfc 100644 --- a/GeneratorInterface/HydjetInterface/interface/HydjetGeneratorFilter.h +++ b/GeneratorInterface/HydjetInterface/interface/HydjetGeneratorFilter.h @@ -5,6 +5,14 @@ #include "GeneratorInterface/Core/interface/GeneratorFilter.h" #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h" +namespace edm { + template <> + GeneratorFilter::GeneratorFilter(ParameterSet const& ps) + : hadronizer_(ps, consumesCollector()) { + init(ps); + } +} // namespace edm + namespace gen { typedef edm::GeneratorFilter HydjetGeneratorFilter; } diff --git a/GeneratorInterface/HydjetInterface/interface/HydjetHadronizer.h b/GeneratorInterface/HydjetInterface/interface/HydjetHadronizer.h index be93ac47c81ec..0e6b8bde80ce7 100644 --- a/GeneratorInterface/HydjetInterface/interface/HydjetHadronizer.h +++ b/GeneratorInterface/HydjetInterface/interface/HydjetHadronizer.h @@ -1,19 +1,24 @@ #ifndef HydjetHadronizer_h #define HydjetHadronizer_h -/** \class HydjetHadronizer -* -* Generates HYDJET ==> HepMC events -* -* Yetkin Yilmaz -* for the Generator Interface. May 2009 -*********************************************/ +/** + \class HydjetHadronizer + \brief Interface to the HYDJET generator (since core v. 1.9.1), produces HepMC events + \version 2.0 + \authors Camelia Mironov, Yetkin Yilmaz, Andrey Belyaev +*/ #define PYCOMP pycomp_ #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "GeneratorInterface/Core/interface/BaseHadronizer.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "SimDataFormats/CrossingFrame/interface/MixCollection.h" +#include "SimDataFormats/PileupSummaryInfo/interface/PileupMixingContent.h" + #include #include #include @@ -27,6 +32,7 @@ namespace HepMC { class GenEvent; class GenParticle; class GenVertex; + class FourVector; } // namespace HepMC namespace gen { @@ -34,7 +40,7 @@ namespace gen { class HydjetHadronizer : public BaseHadronizer { public: - HydjetHadronizer(const edm::ParameterSet&); + HydjetHadronizer(const edm::ParameterSet&, edm::ConsumesCollector&&); ~HydjetHadronizer() override; bool generatePartonsAndHadronize(); @@ -68,68 +74,69 @@ namespace gen { HepMC::GenEvent* evt; edm::ParameterSet pset_; - double abeamtarget_; // beam/target atomic mass number - int angularspecselector_; // angular emitted gluon spectrum selection - double bfixed_; // fixed impact param (fm); valid only if cflag_=0 - double bmax_; // max impact param; - // units of nucl radius - double bmin_; // min impact param; - // units of nucl radius - int cflag_; // centrality flag - // = 0 fixed impact param, - // <> 0 between bmin and bmax - bool embedding_; // Switch for embedding mode - double comenergy; // collision energy - bool doradiativeenloss_; //! DEFAULT = true - bool docollisionalenloss_; //! DEFAULT = true - double fracsoftmult_; // fraction of soft hydro induced hadronic multiplicity - // proportional to no of nucleon participants - // (1-fracsoftmult_)--- fraction of soft - // multiplicity proportional to the numebr - // of nucleon-nucleon binary collisions - // DEFAULT=1., allowed range [0.01,1] - double hadfreeztemp_; // hadron freez-out temperature - // DEFAULT=0.14MeV, allowed ranges [0.08,0.2]MeV - std::string hymode_; // Hydjet running mode - unsigned int maxEventsToPrint_; // Events to print if verbosity - double maxlongy_; // max longitudinal collective rapidity: - // controls width of eta-spectra - // DEFAULT=4, allowed range [0.01,7.0] - double maxtrany_; // max transverse collective rapidity: - // controls slope of low-pt spectra - // DEFAULT=1.5, allowed range [0.01,3.0] - int nsub_; // number of sub-events - int nhard_; // multiplicity of PYTHIA(+PYQUEN)-induced particles in event - int nmultiplicity_; // mean soft multiplicity in central PbPb - // automatically calculated for other centralitie and beams - int nsoft_; // multiplicity of HYDRO-induced particles in event - unsigned int nquarkflavor_; //! number of active quark flavors in qgp - //! DEFAULT=0; allowed values: 0,1,2,3. - unsigned int pythiaPylistVerbosity_; // pythia verbosity; def=1 - double qgpt0_; // initial temperature of QGP - // DEFAULT = 1GeV; allowed range [0.2,2.0]GeV; - double qgptau0_; // proper time of QGP formation - // DEFAULT = 0.1 fm/c; allowed range [0.01,10.0]fm/ - - double phi0_; // Event plane angle + double abeamtarget_; ///< beam/target atomic mass number + int angularspecselector_; ///< angular emitted gluon spectrum selection + double bfixed_; ///< fixed impact param (fm); valid only if cflag_=0 + double bmax_; ///< max impact param; + ///< units of nucl radius + double bmin_; ///< min impact param; + ///< units of nucl radius + int cflag_; ///< centrality flag + ///< = 0 fixed impact param, + ///< <> 0 between bmin and bmax + bool embedding_; ///< Switch for embedding mode + double comenergy; ///< collision energy + bool doradiativeenloss_; ///< DEFAULT = true + bool docollisionalenloss_; ///< DEFAULT = true + double fracsoftmult_; ///< fraction of soft hydro induced hadronic multiplicity + ///< proportional to no of nucleon participants + ///< (1-fracsoftmult_)--- fraction of soft + ///< multiplicity proportional to the numebr + ///< of nucleon-nucleon binary collisions + ///< DEFAULT=1., allowed range [0.01,1] + double hadfreeztemp_; ///< hadron freez-out temperature + ///< DEFAULT=0.14MeV, allowed ranges [0.08,0.2]MeV + std::string hymode_; ///< Hydjet running mode + unsigned int maxEventsToPrint_; ///< Events to print if verbosity + double maxlongy_; ///< max longitudinal collective rapidity: + ///< controls width of eta-spectra + ///< DEFAULT=4, allowed range [0.01,7.0] + double maxtrany_; ///< max transverse collective rapidity: + ///< controls slope of low-pt spectra + ///< DEFAULT=1.5, allowed range [0.01,3.0] + int nsub_; ///< number of sub-events + int nhard_; ///< multiplicity of PYTHIA(+PYQUEN)-induced particles in event + int nmultiplicity_; ///< mean soft multiplicity in central PbPb + ///< automatically calculated for other centralitie and beams + int nsoft_; ///< multiplicity of HYDRO-induced particles in event + unsigned int nquarkflavor_; ///< number of active quark flavors in qgp + ///< DEFAULT=0; allowed values: 0,1,2,3. + unsigned int pythiaPylistVerbosity_; ///< pythia verbosity; def=1 + double qgpt0_; ///< initial temperature of QGP + ///< DEFAULT = 1GeV; allowed range [0.2,2.0]GeV; + double qgptau0_; ///< proper time of QGP formation + ///< DEFAULT = 0.1 fm/c; allowed range [0.01,10.0]fm/ + + double phi0_; ///< Event plane angle double sinphi0_; double cosphi0_; - bool rotate_; // Switch to rotate event plane + bool rotate_; ///< Switch to rotate event plane - unsigned int shadowingswitch_; // shadowing switcher - // 1-ON, 0-OFF - double signn_; // inelastic nucleon nucleon cross section [mb] - // DEFAULT= 58 mb + unsigned int shadowingswitch_; ///< shadowing switcher + ///< 1-ON, 0-OFF + double signn_; ///< inelastic nucleon nucleon cross section [mb] + ///< DEFAULT= 58 mb + HepMC::FourVector* fVertex_; ///< Event signal vertex + std::vector signalVtx_; ///< Pset double vector to set event signal vertex Pythia6Service* pythia6Service_; - edm::InputTag src_; + edm::EDGetTokenT > src_; }; - double HydjetHadronizer::nuclear_radius() const { - // Return the nuclear radius derived from the - // beam/target atomic mass number. - - return 1.15 * pow((double)abeamtarget_, 1. / 3.); - } + /** + Return the nuclear radius derived from the + beam/target atomic mass number. + */ + double HydjetHadronizer::nuclear_radius() const { return 1.15 * pow((double)abeamtarget_, 1. / 3.); } } // namespace gen #endif diff --git a/GeneratorInterface/HydjetInterface/interface/HydjetWrapper.h b/GeneratorInterface/HydjetInterface/interface/HydjetWrapper.h index 11138e8968594..96cfa94b72e65 100644 --- a/GeneratorInterface/HydjetInterface/interface/HydjetWrapper.h +++ b/GeneratorInterface/HydjetInterface/interface/HydjetWrapper.h @@ -1,16 +1,11 @@ #ifndef GeneratorInterface_HydjetInterface_HydjetWrapper #define GeneratorInterface_HydjetInterface_HydjetWrapper -// -// - -/* - * - * Wrapper for FORTRAN version of HYDJET - * - * Camelia Mironov - * - */ +/** + \brief Wrapper for FORTRAN version of HYDJET + \version 2.0 + \authors Camelia Mironov +*/ extern "C" { void hyinit_(double& energy, double& a, int& ifb1, double& bmin, double& bmax, double& bfix1, int& nh1); @@ -20,10 +15,15 @@ void hyinit_(double& energy, double& a, int& ifb1, double& bmin, double& bmax, d #define _MAXMULsize_ 150000 extern "C" { -void hyevnt_(); +void hyevnt_(double& bfix1); } #define HYEVNT hyevnt_ +extern "C" { +void hyjver_(int&, int&, int&, int&); +} +#define HYJVER hyjver_ + extern "C" { extern struct { double psi3; } hypsi3_; } diff --git a/GeneratorInterface/HydjetInterface/python/hydjetDefaultParameters_cff.py b/GeneratorInterface/HydjetInterface/python/hydjetDefaultParameters_cff.py index 6f97f67eb6c16..d1f3b31cc3972 100644 --- a/GeneratorInterface/HydjetInterface/python/hydjetDefaultParameters_cff.py +++ b/GeneratorInterface/HydjetInterface/python/hydjetDefaultParameters_cff.py @@ -16,7 +16,7 @@ comEnergy = cms.double(5020.0) ) -qgpParameters = cms.PSet(qgpInitialTemperature = cms.double(1.1), ## initial temperature of QGP; allowed range [0.2,2.0]GeV; +qgpParameters = cms.PSet(qgpInitialTemperature = cms.double(1.), ## initial temperature of QGP; allowed range [0.2,2.0]GeV; qgpProperTimeFormation = cms.double(0.1), ## proper time of QGP formation; allowed range [0.01,10.0]fm/c; hadronFreezoutTemperature = cms.double(0.125), doRadiativeEnLoss = cms.bool(True), ## if true, perform partonic radiative en loss @@ -27,20 +27,20 @@ hydjetParameters = cms.PSet(sigmaInelNN = cms.double(70), shadowingSwitch = cms.int32(1), - nMultiplicity = cms.int32(25000), + nMultiplicity = cms.int32(18545), fracSoftMultiplicity = cms.double(1.), - maxLongitudinalRapidity = cms.double(4.5), - maxTransverseRapidity = cms.double(1.25), + maxLongitudinalRapidity = cms.double(3.75), + maxTransverseRapidity = cms.double(1.160), rotateEventPlane = cms.bool(True), allowEmptyEvents = cms.bool(False), - angularSpectrumSelector = cms.int32(1), ## angular emitted gluon spectrum : + angularSpectrumSelector = cms.int32(0), ## angular emitted gluon spectrum : embeddingMode = cms.bool(False) ) pyquenPythiaDefaultBlock = cms.PSet( pythiaUESettingsBlock, hydjetPythiaDefault = cms.vstring('MSEL=0 ! user processes', - 'CKIN(3)=10.',# ! ptMin + 'CKIN(3)=9.2',# ! ptMin 'MSTP(81)=1' ), ppDefault = cms.vstring('MSEL=1 ! QCD hight pT processes (only jets)', diff --git a/GeneratorInterface/HydjetInterface/src/HydjetHadronizer.cc b/GeneratorInterface/HydjetInterface/src/HydjetHadronizer.cc index 679142892249c..9f6fd7d565109 100644 --- a/GeneratorInterface/HydjetInterface/src/HydjetHadronizer.cc +++ b/GeneratorInterface/HydjetInterface/src/HydjetHadronizer.cc @@ -1,14 +1,8 @@ -/* - - ######################## - # Hydjet1 # - # version: 1.9 patch1 # - ######################## - - * Interface to the HYDJET generator, produces HepMC events - * - * Original Author: Camelia Mironov - */ +/** + \brief Interface to the HYDJET generator (since core v. 1.9.1), produces HepMC events + \version 2.0 + \authors Camelia Mironov, Andrey Belyaev +*/ #include #include @@ -28,6 +22,7 @@ #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Declarations.h" #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h" +#include "HepMC/IO_HEPEVT.h" #include "HepMC/PythiaWrapper6_4.h" #include "HepMC/GenEvent.h" #include "HepMC/HeavyIon.h" @@ -42,6 +37,8 @@ using namespace edm; using namespace std; using namespace gen; +HepMC::IO_HEPEVT hepevtio; + namespace { int convertStatus(int st) { if (st <= 0) @@ -61,7 +58,7 @@ const std::vector HydjetHadronizer::theSharedResources = {edm::Shar gen::FortranInstance::kFortranInstance}; //_____________________________________________________________________ -HydjetHadronizer::HydjetHadronizer(const ParameterSet& pset) +HydjetHadronizer::HydjetHadronizer(const ParameterSet& pset, edm::ConsumesCollector&& iC) : BaseHadronizer(pset), evt(nullptr), pset_(pset), @@ -95,9 +92,22 @@ HydjetHadronizer::HydjetHadronizer(const ParameterSet& pset) rotate_(pset.getParameter("rotateEventPlane")), shadowingswitch_(pset.getParameter("shadowingSwitch")), signn_(pset.getParameter("sigmaInelNN")), + fVertex_(nullptr), pythia6Service_(new Pythia6Service(pset)) { // Default constructor + if (pset.exists("signalVtx")) + signalVtx_ = pset.getUntrackedParameter >("signalVtx"); + + if (signalVtx_.size() == 4) { + if (!fVertex_) + fVertex_ = new HepMC::FourVector(); + LogDebug("EventSignalVertex") << "Setting event signal vertex " + << " x = " << signalVtx_.at(0) << " y = " << signalVtx_.at(1) + << " z= " << signalVtx_.at(2) << " t = " << signalVtx_.at(3) << endl; + fVertex_->set(signalVtx_.at(0), signalVtx_.at(1), signalVtx_.at(2), signalVtx_.at(3)); + } + // PYLIST Verbosity Level // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13 pythiaPylistVerbosity_ = pset.getUntrackedParameter("pythiaPylistVerbosity", 0); @@ -107,8 +117,15 @@ HydjetHadronizer::HydjetHadronizer(const ParameterSet& pset) maxEventsToPrint_ = pset.getUntrackedParameter("maxEventsToPrint", 0); LogDebug("Events2Print") << "Number of events to be printed = " << maxEventsToPrint_; - if (embedding_) - src_ = pset.getParameter("backgroundLabel"); + if (embedding_) { + cflag_ = 0; + src_ = iC.consumes >( + pset.getUntrackedParameter("backgroundLabel", edm::InputTag("mix", "generatorSmeared"))); + } + + int cm = 1, va, vb, vc; + HYJVER(cm, va, vb, vc); + HepMC::HEPEVT_Wrapper::set_max_number_entries(4000); } //_____________________________________________________________________ @@ -150,7 +167,6 @@ void HydjetHadronizer::add_heavy_ion_rec(HepMC::GenEvent* evt) { //___________________________________________________________________ HepMC::GenParticle* HydjetHadronizer::build_hyjet(int index, int barcode) { // Build particle object corresponding to index in hyjets (soft+hard) - double x0 = hyjets.phj[0][index]; double y0 = hyjets.phj[1][index]; @@ -172,7 +188,6 @@ HepMC::GenParticle* HydjetHadronizer::build_hyjet(int index, int barcode) { //___________________________________________________________________ HepMC::GenVertex* HydjetHadronizer::build_hyjet_vertex(int i, int id) { // build verteces for the hyjets stored events - double x0 = hyjets.vhj[0][i]; double y0 = hyjets.vhj[1][i]; double x = x0 * cosphi0_ - y0 * sinphi0_; @@ -191,20 +206,53 @@ bool HydjetHadronizer::generatePartonsAndHadronize() { // generate single event if (embedding_) { - cflag_ = 0; const edm::Event& e = getEDMEvent(); - Handle input; - e.getByLabel(src_, input); - const HepMC::GenEvent* inev = input->GetEvent(); + HepMC::GenVertex* genvtx = nullptr; + const HepMC::GenEvent* inev = nullptr; + Handle > cf; + e.getByToken(src_, cf); + MixCollection mix(cf.product()); + if (mix.size() < 1) { + throw cms::Exception("MatchVtx") << "Mixing has " << mix.size() << " sub-events, should have been at least 1" + << endl; + } + const HepMCProduct& bkg = mix.getObject(0); + if (!(bkg.isVtxGenApplied())) { + throw cms::Exception("MatchVtx") << "Input background does not have smeared vertex!" << endl; + } else { + inev = bkg.GetEvent(); + } + + genvtx = inev->signal_process_vertex(); + + if (!genvtx) + throw cms::Exception("MatchVtx") << "Input background does not have signal process vertex!" << endl; + + double aX, aY, aZ, aT; + + aX = genvtx->position().x(); + aY = genvtx->position().y(); + aZ = genvtx->position().z(); + aT = genvtx->position().t(); + + if (!fVertex_) { + fVertex_ = new HepMC::FourVector(); + } + LogInfo("MatchVtx") << " setting vertex " + << " aX " << aX << " aY " << aY << " aZ " << aZ << " aT " << aT << endl; + fVertex_->set(aX, aY, aZ, aT); + const HepMC::HeavyIon* hi = inev->heavy_ion(); + if (hi) { - bfixed_ = hi->impact_parameter(); + bfixed_ = (hi->impact_parameter()) / nuclear_radius(); phi0_ = hi->event_plane_angle(); sinphi0_ = sin(phi0_); cosphi0_ = cos(phi0_); } else { LogWarning("EventEmbedding") << "Background event does not have heavy ion record!"; } + } else if (rotate_) rotateEvtPlane(); @@ -217,7 +265,6 @@ bool HydjetHadronizer::generatePartonsAndHadronize() { edm::LogInfo("HYDJETinTemp") << "##### HYDJET: QGP init temperature, T0 =" << pyqpar.T0u; edm::LogInfo("HYDJETinTau") << "##### HYDJET: QGP formation time,tau0 =" << pyqpar.tau0u; - // generate a HYDJET event int ntry = 0; while (nsoft_ == 0 && nhard_ == 0) { if (ntry > 100) { @@ -231,7 +278,7 @@ bool HydjetHadronizer::generatePartonsAndHadronize() { edm::Exception except(edm::errors::EventCorruption, sstr.str()); throw except; } else { - HYEVNT(); + HYEVNT(bfixed_); nsoft_ = hyfpar.nhyd; nsub_ = hyjpar.njet; nhard_ = hyfpar.npyt; @@ -252,6 +299,21 @@ bool HydjetHadronizer::generatePartonsAndHadronize() { evt->set_event_scale(pypars.pari[16]); // Q^2 add_heavy_ion_rec(evt); + if (fVertex_) { + // generate new vertex & apply the shift + + // Copy the HepMC::GenEvent + std::unique_ptr HepMCEvt(new edm::HepMCProduct(evt)); + + HepMCEvt->applyVtxGen(fVertex_); + evt = new HepMC::GenEvent((*HepMCEvt->GetEvent())); + } + + HepMC::HEPEVT_Wrapper::check_hepevt_consistency(); + LogDebug("HEPEVT_info") << "Ev numb: " << HepMC::HEPEVT_Wrapper::event_number() + << " Entries number: " << HepMC::HEPEVT_Wrapper::number_entries() << " Max. entries " + << HepMC::HEPEVT_Wrapper::max_number_entries() << std::endl; + event().reset(evt); return true; } @@ -266,76 +328,90 @@ bool HydjetHadronizer::get_particles(HepMC::GenEvent* evt) { // The SubEvent information is kept by storing indeces of main vertices // of subevents as a vector in GenHIEvent. - LogDebug("SubEvent") << "Number of sub events " << nsub_; - LogDebug("Hydjet") << "Number of hard events " << hyjpar.njet; - LogDebug("Hydjet") << "Number of hard particles " << nhard_; - LogDebug("Hydjet") << "Number of soft particles " << nsoft_; - - vector sub_vertices(nsub_); + LogDebug("SubEvent") << " Number of sub events " << nsub_; + LogDebug("Hydjet") << " Number of hard events " << hyjpar.njet; + LogDebug("Hydjet") << " Number of hard particles " << nhard_; + LogDebug("Hydjet") << " Number of soft particles " << nsoft_; + LogDebug("Hydjet") << " nhard_ + nsoft_ = " << nhard_ + nsoft_ << " hyjets.nhj = " << hyjets.nhj << endl; int ihy = 0; - for (int isub = 0; isub < nsub_; isub++) { - LogDebug("SubEvent") << "Sub Event ID : " << isub; - - int sub_up = (isub + 1) * 50000; // Upper limit in mother index, determining the range of Sub-Event - vector particles; - vector mother_ids; - vector prods; - - sub_vertices[isub] = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), isub); - evt->add_vertex(sub_vertices[isub]); - if (!evt->signal_process_vertex()) - evt->set_signal_process_vertex(sub_vertices[isub]); - - while (ihy < nhard_ + nsoft_ && (hyjets.khj[2][ihy] < sub_up || ihy > nhard_)) { - particles.push_back(build_hyjet(ihy, ihy + 1)); - prods.push_back(build_hyjet_vertex(ihy, isub)); - mother_ids.push_back(hyjets.khj[2][ihy]); - LogDebug("DecayChain") << "Mother index : " << hyjets.khj[2][ihy]; - - ihy++; - } + int isub = -1; + int isub_l = -1; + int stab = 0; - //Produce Vertices and add them to the GenEvent. Remember that GenParticles are adopted by - //GenVertex and GenVertex is adopted by GenEvent. + vector primary_particle(hyjets.nhj); + vector particle(hyjets.nhj); - LogDebug("Hydjet") << "Number of particles in vector " << particles.size(); + HepMC::GenVertex* sub_vertices = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), 0); // just initialization - for (unsigned int i = 0; i < particles.size(); i++) { - HepMC::GenParticle* part = particles[i]; + // contain the last index in for each subevent + vector index(nsub_); - //The Fortran code is modified to preserve mother id info, by seperating the beginning - //mother indices of successive subevents by 5000 - int mid = mother_ids[i] - isub * 50000 - 1; - LogDebug("DecayChain") << "Particle " << i; - LogDebug("DecayChain") << "Mother's ID " << mid; - LogDebug("DecayChain") << "Particle's PDG ID " << part->pdg_id(); + while (ihy < hyjets.nhj) { + isub = std::floor((hyjets.khj[2][ihy] / 50000)); + int hjoffset = isub * 50000; - if (mid <= 0) { - sub_vertices[isub]->add_particle_out(part); - continue; + if (isub != isub_l) { + sub_vertices = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), isub); + evt->add_vertex(sub_vertices); + if (!evt->signal_process_vertex()) + evt->set_signal_process_vertex(sub_vertices); + index[isub] = ihy - 1; + isub_l = isub; + } + + if (convertStatus(hyjets.khj[0][ihy]) == 1) + stab++; + LogDebug("Hydjet_array") << ihy << " MULTin ev.:" << hyjets.nhj << " SubEv.#" << isub << " Part #" << ihy + 1 + << ", PDG: " << hyjets.khj[1][ihy] << " (st. " << convertStatus(hyjets.khj[0][ihy]) + << ") mother=" << hyjets.khj[2][ihy] - (isub * 50000) + index[isub] + 1 << " (" + << hyjets.khj[2][ihy] << "), childs (" + << hyjets.khj[3][ihy] - (isub * 50000) + index[isub] + 1 << "-" + << hyjets.khj[4][ihy] - (isub * 50000) + index[isub] + 1 << "), vtx (" + << hyjets.vhj[0][ihy] << "," << hyjets.vhj[1][ihy] << "," << hyjets.vhj[2][ihy] << ") " + << std::endl; + + if (hyjets.khj[2][ihy] == 0) { + primary_particle[ihy] = build_hyjet(ihy, ihy + 1); + sub_vertices->add_particle_out(primary_particle[ihy]); + LogDebug("Hydjet_array") << " ---> " << ihy + 1 << std::endl; + } else { + particle[ihy] = build_hyjet(ihy, ihy + 1); + int mid = hyjets.khj[2][ihy] - hjoffset + index[isub]; + int mid_t = mid; + while ((mid < ihy) && (hyjets.khj[1][mid] < 100) && (hyjets.khj[3][mid + 1] - hjoffset + index[isub] == ihy)) + mid++; + if (hyjets.khj[1][mid] < 100) + mid = mid_t; + + HepMC::GenParticle* mother = primary_particle.at(mid); + HepMC::GenVertex* prods = build_hyjet_vertex(ihy, isub); + + if (!mother) { + mother = particle[mid]; + primary_particle[mid] = mother; } - if (mid > 0) { - HepMC::GenParticle* mother = particles[mid]; - LogDebug("DecayChain") << "Mother's PDG ID " << mother->pdg_id(); - - HepMC::GenVertex* prod_vertex = mother->end_vertex(); - if (!prod_vertex) { - prod_vertex = prods[i]; - prod_vertex->add_particle_in(mother); - evt->add_vertex(prod_vertex); - prods[i] = nullptr; // mark to protect deletion - } - prod_vertex->add_particle_out(part); + HepMC::GenVertex* prod_vertex = mother->end_vertex(); + if (!prod_vertex) { + prod_vertex = prods; + prod_vertex->add_particle_in(mother); + LogDebug("Hydjet_array") << " <--- " << mid + 1 << std::endl; + evt->add_vertex(prod_vertex); + prods = nullptr; } + + prod_vertex->add_particle_out(particle[ihy]); + LogDebug("Hydjet_array") << " ---" << mid + 1 << "---> " << ihy + 1 << std::endl; + + if (prods) + delete prods; } - // cleanup vertices not assigned to evt - for (unsigned int i = 0; i < prods.size(); i++) { - if (prods[i]) - delete prods[i]; - } + ihy++; } + LogDebug("Hydjet_array") << " MULTin ev.:" << hyjets.nhj << ", last index: " << ihy - 1 + << ", Sub events: " << isub + 1 << ", stable particles: " << stab << std::endl; + return true; } diff --git a/GeneratorInterface/HydjetInterface/src/hydjet.F b/GeneratorInterface/HydjetInterface/src/hydjet.F deleted file mode 100644 index af6508825bb21..0000000000000 --- a/GeneratorInterface/HydjetInterface/src/hydjet.F +++ /dev/null @@ -1,1676 +0,0 @@ -*---------------------------------------------------------------------- -* -* Filename : HYDJET1_9.F -* -* Author : Igor Lokhtin (Igor.Lokhtin@cern.ch) -* Version : HYDJET1_9.f -* Last revision : 1-JUL-2013 -* -*====================================================================== -* -* Description : Event generator for simulation of jet production, jet -* quenching and flow effects in ultrarelativistic heavy -* ion AA collisions -* -* Method: I.P. Lokhtin and A.M. Snigirev, Eur. Phys. J C 45 (2006) 211 -* -*====================================================================== - -********************************* HYINIT *************************** - SUBROUTINE HYINIT(energy,A,ifb1,bmin,bmax,bfix1,nh1) -* PYTHIA inizialization, calculation of total and hard cross sections and -* # of participants and binary sub-collisions at "reference point" (Pb,b=0), -* tabulation of nuclear thickness function and nuclear overlap function, -* test of input parametes - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - double precision numpar,npar0,nbco0 - external numpar,rhoaa,hfunc3 - common /pyint7/ sigt(0:6,0:6,0:5) - common /pypars/ mstp(200),parp(200),msti(200),pari(200) - common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200) - common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet - common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr - common /hypyin/ ene,rzta,rnta,bfix,ifb,nh - common /hyflow/ ytfl,ylfl,Tf,fpart - common /hygeom/ BC - common /hythic/ BAB(110),TAB(110),TAAB(110) - common /hynup1/ bp,x - save /pyint7/,/pypars/,/pysubs/,/hyjpar/,/hyipar/,/hypyin/, - > /hyflow/,/hygeom/,/hythic/ - -* start HYDJET initialization - init=1 - ipr=1 - -* set beam paramters - ene=energy ! c.m.s. energy per nucleon - AW=A ! atomic weight - RA=1.15d0*AW**0.333333d0 ! nuclear radius - rzta=1.d0/(1.98d0+0.015d0*A**0.666667d0) ! fraction of protons in nucleus - rnta=1.d0-rzta ! fraction of neutrons in nucleus - ifb=ifb1 ! centrality flag - bfix=bfix1 ! fixed impact parameter - nh=nh1 ! mean soft mult. in central PbPb - - ptmin=ckin(3) ! minimum pt of hard scattering - -* printout of HYDJET versioning - write(6,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - >%%%%%%%%%%%%%%%%%%%%%%%%' - write(6,*) '%% This is HYDJET version 1.9 - > %%' - write(6,*) '%% Corresponding author: Igor Lokhtin (Igor.Lo - >khtin@cern.ch) %%' - write(6,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - >%%%%%%%%%%%%%%%%%%%%%%%%' - -* Pythia inizialization for pp collisions - call pyinit('cms','p','p',ene) - -* calculation of hard scattering cross section - if(nhsel.ne.0) then - - mstp(111)=0 -* no printout of Pythia initialization information hereinafter - mstp(122)=0 -* Pythia test pp event run - do i=1,1000 - call pyevnt - end do -* hard scattering pp cross section - sjpp=pari(1) - -* Pythia inizialization for pn collisions - call pyinit('cms','p','n',ene) -* Pythia test pn event run - do i=1,1000 - call pyevnt - end do -* hard scattering pn cross section - sjpn=pari(1) - -* Pythia inizialization for nn collisions - call pyinit('cms','n','n',ene) -* Pythia test nn event run - do i=1,1000 - call pyevnt - end do -* hard scattering nn cross section - sjnn=pari(1) - - sigjet=rzta*rzta*sjpp+rnta*rnta*sjnn+2.d0*rzta*rnta*sjpn - - end if - -* total inelastic cross section - if(sigin.lt.10.d0.or.sigin.gt.200.d0) - > sigin=sigt(0,0,0)-sigt(0,0,1) - -* # of nucelons-participants and NN sub-collisions at "reference point" (Pb,b=0) - Apb=207.d0 - Rpb=1.15d0*Apb**0.333333d0 - EPS=0.005d0 - Z2=4.d0*Rpb - Z1=-1.d0*Z2 - H=0.01d0*(Z2-Z1) - do ib=1,110 - BC=3.d0*Rpb*(ib-1)/109.d0 - CALL SIMPA(Z1,Z2,H,EPS,1.d-8,rhoaa,Z,RES,AIH,AIABS) - BAB(ib)=BC - TAB(ib)=Apb*RES - end do - Z1=0.d0 - Z2=6.28318d0 - H=0.01d0*(Z2-Z1) - bp=0.d0 - CALL SIMPA(Z1,Z2,H,EPS,1.d-8,HFUNC3,X,TAAPB0,AIH,AIABS) - - npar0=numpar(0.d0) ! Npart(Pb,b=0) - nbco0=0.1d0*sigin*TAAPB0 ! Nsub(Pb,b=0) - - init=0 - -* creation of arrays for tabulation of beam/target nuclear thickness function - Z2=4.d0*RA - Z1=-1.d0*Z2 - H=0.01d0*(Z2-Z1) - do ib=1,110 - BC=3.d0*RA*(ib-1)/109.d0 - CALL SIMPA(Z1,Z2,H,EPS,1.d-8,rhoaa,Z,RES,AIH,AIABS) - BAB(ib)=BC - TAB(ib)=AW*RES - end do - -* creation of arrays for tabulation of nuclear overlap function - Z1=0.d0 - Z2=6.28318d0 - H=0.01d0*(Z2-Z1) - do ib=1,110 - bp=BAB(ib) - CALL SIMPA(Z1,Z2,H,EPS,1.d-8,HFUNC3,X,RES,AIH,AIABS) - TAAB(ib)=RES - end do - -* test of centrality selection - if(ifb.eq.0) then - if(bfix.lt.0.d0) then - write(6,*) 'Impact parameter less than zero!' - bfix=0.d0 - end if - if (bfix.gt.3.d0) then - write(6,*) 'Impact parameter larger than three nuclear radius!' - bfix=3.d0 - end if - else - if(bmin.lt.0.d0) then - write(6,*) 'Impact parameter less than zero!' - bmin=0.d0 - end if - if(bmax.gt.3.d0) then - write(6,*) 'Impact parameter larger than three nuclear radius!' - bmax=3.d0 - end if - bminh=bmin - bmaxh=bmax - end if - -* test of flow parameter selection - if (Tf.lt.0.08d0.or.Tf.gt.0.2d0) Tf=0.1d0 ! freeze-out temperature - if (ylfl.lt.0.01d0.or.ylfl.gt.7.d0) ylfl=4.d0 ! longitudinal flow rapidity - if (ytfl.lt.0.01d0.or.ytfl.gt.3.d0) ytfl=1.5d0 ! transverse flow rapidity - if (fpart.le.0.d0.or.fpart.gt.1.d0) fpart=1.d0 ! fraction of soft multiplicity - ! proport. to # of participants -* test of 'nhsel' selection - if(nhsel.ne.1.and.nhsel.ne.2.and.nhsel.ne.3.and.nhsel.ne.4) - > nhsel=0 - - return - end -****************************** END HYINIT *************************** - -********************************* HYEVNT **************************** - SUBROUTINE HYEVNT -* generation of single HYDJET event (soft+hard parts) at given parameters - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - double precision numpar,npar0,nbco0,npart,nbcol - integer pycomp - external hsin,gauss,hftaa,numpar,hyhard,hipsear,pyr,pymass - common /pyjets/ nhi,npadhi,khi(4000,5),phi(4000,5),vhi(4000,5) - COMMON /PYDAT3/ MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5) - common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5) - common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr - common /hypyin/ ene,rzta,rnta,bfix,ifb,nh - common /hyfpar/ bgen,nbcol,npart,npyt,nhyd - common /hypsi3/ psi3 - common /hyflow/ ytfl,ylfl,Tf,fpart - common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet - common /pypars/ mstp(200),parp(200),msti(200),pari(200) - save /pypars/,/pyjets/,/pydat3/,/hyjets/,/hyipar/,/hyfpar/, - + /hyflow/,/hyjpar/,/hypyin/,/hypsi3/ - -* reset lujets and hyjets arrays before event generation -c n=0 - nhj=0 - do ncl=1,150000 - do j=1,5 -c p(ncl,j)=0. - phj(ncl,j)=0.d0 -c v(ncl,j)=0. - vhj(ncl,j)=0.d0 -c k(ncl,j)=0 - khj(ncl,j)=0 - enddo - end do -* - pi=3.14159d0 - -* generate impact parameter of A-A collision - if(ifb.eq.0) then - b1=bfix*RA - bgen=bfix - else - call hipsear(fmax1,xmin1) - fmax=fmax1 - xmin=xmin1 - 3 bb1=xmin*pyr(0)+bminh*RA - ff1=fmax*pyr(0) - fb=hsin(bb1) - if(ff1.gt.fb) goto 3 - b1=bb1 - bgen=bb1/RA - end if - -* calculate # of nucelons-participants and binary NN sub-collisions - npart=numpar(b1) ! Npart(b) - nbcol=0.1d0*sigin*hftaa(b1) ! Nsub(b) - -* generate soft multiplicity in event, np, -* fpart - fraction of soft multiplicity proportional to # of participants, -* fbcol=1-fpart - fraction of multiplicity proprtional to # of NN subcollisions - fbcol=1.d0-fpart - rnp=nh*(fpart*npart+fbcol*nbcol)/(fpart*npar0+fbcol*nbco0) - np=int(rnp) - sign=dsqrt(rnp) - 1 if(nhsel.lt.3.and.np.gt.0) then - np=max(0,int(gauss(rnp,sign))) - else - np=0 - end if - if(np.gt.150000) then - write(6,*) 'Warning, soft multiplicity too large!' - goto 1 - end if - -* generate hard parton-parton scatterings (Q>ptmin) 'njet' times - njet=0 - if(nhsel.ne.0) then - pjet=sigjet/sigin - do i=1,int(nbcol) - if(pyr(0).lt.pjet) njet=njet+1 - end do - call hyhard - end if - -* reset pyjets array before soft component generation - nhi=0 - do ncl=1,4000 - do j=1,5 - phi(ncl,j)=0.d0 - vhi(ncl,j)=0.d0 - khi(ncl,j)=0 - enddo - end do - - npyt=nhj - nhyd=np - -* Offset for separating Hydro as a sub-event. YY - ihydoff = njet * 50000 - - if(nhyd.eq.0) goto 4 - if(bgen.gt.1.999d0) b1=1.999d0*RA - -* set of parameters for azimuthal anisotropy in hadron emission region - anc=2.8d0 - ank=0.25d0 - ane=0.5d0*ank*b1/RA - anb=max(0.000000000001d0,anc*(1.-ane*ane)*ane) - anl=max(0.000000000001d0,0.5d0*(dsqrt(1.d0+4.d0*anb* - > (ane+anb))-1.d0)/anb) - anr=RA*dsqrt(abs(1.d0-ane)) - eps3=0.06*(b1/RA)**0.33333d0 - psi3=2.d0*pi*(-0.5d0+pyr(0))/3.d0 - -* generate sort of hadrons in jetset/pythia format - 5 yy=9999.d0*pyr(0) - if(yy.lt.1376.d0) then - kf=211 - elseif(yy.lt.2752.d0) then - kf=-211 - elseif(yy.lt.4142.d0) then - kf=111 - elseif(yy.lt.4590.d0) then - kf=321 - elseif(yy.lt.5038.d0) then - kf=-321 - elseif(yy.lt.5479.5d0) then - kf=311 - elseif(yy.lt.5921.d0) then - kf=-311 - elseif(yy.lt.6063.5d0) then - kf=2212 - elseif(yy.lt.6206.d0) then - kf=-2212 - elseif(yy.lt.6347.7d0) then - kf=2112 - elseif(yy.lt.6489.4d0) then - kf=-2112 - elseif(yy.lt.6935.4d0) then - kf=213 - elseif(yy.lt.7381.4d0) then - kf=-213 - elseif(yy.lt.7827.4d0) then - kf=113 - elseif(yy.lt.8248.8d0) then - kf=223 - elseif(yy.lt.8614.5d0) then - kf=221 - elseif(yy.lt.8680.d0) then - kf=331 - elseif(yy.lt.8829.d0) then - kf=333 - elseif(yy.lt.9092.d0) then - kf=323 - elseif(yy.lt.9355.d0) then - kf=-323 - elseif(yy.lt.9613.d0) then - kf=313 - elseif(yy.lt.9871.d0) then - kf=-313 - elseif(yy.lt.9935.d0) then - kf=3122 - else - kf=-3122 - end if - nhj=nhj+1 - khj(nhj,1)=1 - khj(nhj,2)=kf - do j=3,5 - khj(nhj,j)=ihydoff - end do - do j=1,5 - vhj(nhj,j)=0.d0 - enddo - kfa=iabs(kf) - phj(nhj,5)=pymass(kfa) - -* generate uniform distribution in system of a fluid element - 2 ep0=-1.d0*Tf*(dlog(max(1.d-10,pyr(0)))+dlog(max(1.d-10,pyr(0))) - > +dlog(max(1.d-10,pyr(0)))) - if(ep0.le.phj(nhj,5)) go to 2 - pp0=dsqrt(abs(ep0**2-abs(phj(nhj,5)**2))) - probt=pp0/ep0 - if(pyr(0).gt.probt) go to 2 - ctp0=2.d0*pyr(0)-1.d0 - stp0=dsqrt(abs(1.d0-ctp0**2)) - php0=2.d0*pi*pyr(0) - -* generate coordinates of a fluid element -c etaf=ylfl*(2.d0*pyr(0)-1.d0) ! flat initial eta-spectrum - etaf=gauss(0.d0,ylfl) ! gaussian initial eta-spectrum - phif=2.d0*pi*pyr(0) - RF=anr*dsqrt(abs(1.d0-ane*ane))/ - > (dsqrt(abs(1.d0+ane*dcos(2.d0*phif)))) - RF=RF*(1.d0+eps3*dcos(3.d0*(phif-psi3))) - rrf=RF*RF*pyr(0) - -* generate four-velocity of a fluid element - vradf=dsinh(ytfl*dsqrt(rrf)/anr) - rn1=dcosh(etaf) - rn2=dsinh(etaf) - uxf=vradf*dsqrt(abs(1.d0+anl))*dcos(phif) - uyf=vradf*dsqrt(abs(1.d0-anl))*dsin(phif) - urf=dsqrt(uxf*uxf+uyf*uyf) - rn0=dsqrt(abs(1.d0+urf*urf)) - utf=rn1*rn0 - uzf=rn2*rn0 - -* Von Neumann rejection/acceptance procedure for Cooper-Frye prescription - rfac=(rn1+rn0)/(1.d0+utf) - rnx=-1.d0*uxf*rfac - rny=-1.d0*uyf*rfac - rnz=rn2-uzf*rfac - pnx=pp0*stp0*dcos(php0) - pny=pp0*stp0*dsin(php0) - pnz=pp0*ctp0 - pn0=dsqrt(pnx**2+pny**2+pnz**2+phj(nhj,5)**2) - wesn=rn0-(rnx*pnx+rny*pny+rnz*pnz)/pn0 - wksi=2.d0*dcosh(ytfl)*pyr(0) - if(wksi.gt.wesn) goto 2 - -* boost in the laboratory system - uipi=pp0*(urf*stp0*dcos(phif-php0)+uzf*ctp0) - bfac=uipi/(utf+1.d0)+ep0 - phj(nhj,1)=pnx+urf*dcos(phif)*bfac - phj(nhj,2)=pny+urf*dsin(phif)*bfac - phj(nhj,3)=pnz+uzf*bfac - phj(nhj,4)=dsqrt(phj(nhj,1)**2+phj(nhj,2)**2+phj(nhj,3)**2 - > +phj(nhj,5)**2) - -* skip if too energetic forward particle - YY - if(abs(phj(nhj,3)).gt.500) then - nhj = nhj-1 - np = np-1 - goto 5 - end if - -* treat hadron decay (if it switched on) with pyexec routine - if(mdcy(pycomp(kf),1).eq.1) then - nhi=1 - do j=1,5 - phi(1,j)=phj(nhj,j) - khi(1,j)=khj(nhj,j) - vhi(1,j)=vhj(nhj,j) - end do - call pyexec - do i=1,nhi - do j=1,5 - icor=0 - if(j.gt.2.and.khi(i,j).ne.ihydoff) - > icor=ihydoff + nhj - npyt -1 - phj(nhj+i-1,j)=phi(i,j) - khj(nhj+i-1,j)=khi(i,j)+icor - vhj(nhj+i-1,j)=vhi(i,j) - end do - end do - nhj=nhj+nhi-1 - np=np+nhi-1 - end if - if(nhj.lt.npyt+np) goto 5 - nhyd=np - - 4 continue - -* fill array 'lujets' -c n=nhyd+npyt -c do ih=1,nhj -c do jh=1,5 -c p(ih,jh)=phj(ih,jh) -c k(ih,jh)=khj(ih,jh) -c v(ih,jh)=vhj(ih,jh) -c end do -c end do - - return - end -****************************** END HYEVNT ************************** - -********************************* HYHARD *************************** - SUBROUTINE HYHARD -* generate 'njet' number of hard parton-parton scatterings with PYTHIA - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - double precision npar0,nbco0,npart,nbcol - CHARACTER beam*2,targ*2 - external pydata - external pyp,pyr,pyk,pyquen,shad1 - common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5) - common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5) - COMMON /PYDAT1/ MSTU(200),PARU(200),MSTJ(200),PARJ(200) - COMMON /PYDAT2/ KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) - COMMON /PYDAT3/ MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5) - COMMON /PYSUBS/ MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200) - common /pypars/ mstp(200),parp(200),msti(200),pari(200) - common /parimp/ b1, psib1, rb1, rb2, noquen - common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet - common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr - common /hyfpar/ bgen,nbcol,npart,npyt,nhyd - common /hypyin/ ene,rzta,rnta,bfix,ifb,nh - save /pyjets/,/pypars/,/pydat1/,/pydat2/,/pydat3/,/pysubs/, - + /hyjets/,/parimp/,/hyjpar/,/hyipar/,/hyfpar/,/hypyin/ - - -* generate 'njet' PYTHIA events and fill arrays for partons and hadrons - xhardestPtHat=-0.1 - nshad=0 - noquen=0 - if(nhsel.eq.1.or.nhsel.eq.3) noquen=1 - if(njet.ge.1) then - ifbp=0 ! fix impact parameter - bfixp=RA*bgen - Ap=AW ! atomic weight - - do ihard=1,njet - mstp(111)=0 - -* generate type of NN sub-collision (pp, pn or nn) - rand1=pyr(0) - if(rand1.lt.rzta) then - beam='p' - else - beam='n' - end if - rand2=pyr(0) - if(rand2.lt.rzta) then - targ='p' - else - targ='n' - end if - call pyinit('cms',beam,targ,ene) -c mstj(41)=0 ! vacuum showering off - call pyevnt ! generate hard scattering - -* PYQUEN: quenched jets if noquen=0 or non-quenched jets if noquen=1 and ishad=1 - if(ishad.eq.1.or.nhsel.eq.2.or.nhsel.eq.4) - > call pyquen(Ap,ifbp,bfixp) - -* treatment of "nuclear shadowing" (for Pb, Au, Pd or Ca beams only) - if(ishad.eq.1) then - kfh1=abs(k(3,2)) - kfh2=abs(k(4,2)) - xh1=pari(33) - xh2=pari(34) - Q2=pari(22) - shad=shad1(kfh1,xh1,Q2,rb1)*shad1(kfh2,xh2,Q2,rb2) - if(pyr(0).gt.shad) then - nshad=nshad+1 - goto 53 - end if - end if - - call pyexec ! hadronization done -c call pyedit(2) ! remove partons & leave hadrons - if(pari(17).gt.xhardestPtHat) xhardestPtHat = pari(17) - -* fill array of final particles -#ifndef __HIROOT_FORTRAN_ONLY__ - ievoff=(ihard-nshad-1)*50000 ! event offset -#endif - nu=nhj+n - if(nu.gt.150000-np) then - write(6,*) 'Warning, multiplicity too large! Cut hard part.' - goto 52 - end if - nhj=nu - do i=nhj-n+1,nhj - ip=i+n-nhj - do j=1,5 - phj(i,j)=p(ip,j) - end do - do j=1,5 - vhj(i,j)=v(ip,j) - end do - do j=1,5 - khj(i,j)=k(ip,j) - end do - do j=3,5 - kk=k(ip,j) -#ifndef __HIROOT_FORTRAN_ONLY__ - khj(i,j)=kk+ievoff -#else - if(kk.gt.0) khj(i,j)=kk+nhj-n -#endif - end do - end do - 53 continue - end do - 52 njet=ihard-1 - end if - njet=njet-nshad - - pari(17)=xhardestPtHat - - return - end -****************************** END HYHARD ************************** - -********************************* HIPSEAR *************************** - SUBROUTINE HIPSEAR (fmax,xmin) -* find maximum and 'sufficient minimum' of differential inelasic AA cross -* section as a function of impact paramater (xm, fm are outputs) - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - double precision npar0,nbco0 - external hsin - common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr - save /hyipar/ - xmin=(bmaxh-bminh)*RA - fmax=0.d0 - do j=1,1000 - x=bminh*RA+xmin*(j-1)/999.d0 - f=hsin(x) - if(f.gt.fmax) then - fmax=f - endif - end do - return - end -****************************** END HIPSEAR ************************** - -* differential inelastic AA cross section - double precision function hsin(x) - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - external hftaa - common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet - save /hyjpar/ - br=x - hsin=br*(1.d0-dexp(-0.1d0*hftaa(br)*sigin)) - return - end - -* number of nucleons-participants at impact parameter b - double precision function numpar(c) - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - external HFUNC1 - common /hynup1/ bp,x - EPS=0.005d0 - A=0.d0 - B=6.28318d0 - H=0.01d0*(B-A) - bp=c - CALL SIMPA(A,B,H,EPS,1.d-8,HFUNC1,X,RES,AIH,AIABS) - numpar=2.d0*RES - return - end -* - double precision function HFUNC1(x) - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - double precision npar0,nbco0 - external HFUNC2 - common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr - common /hynup1/ bp,xx - save /hyipar/ - if(init.eq.1) then - Rl=Rpb - else - Rl=RA - end if - xx=x - EPS=0.005d0 - A=0.d0 - B=3.d0*Rl - H=0.01d0*(B-A) - CALL SIMPB(A,B,H,EPS,1.d-8,HFUNC2,Y,RES,AIH,AIABS) - HFUNC1=RES - return - end -* - double precision function HFUNC2(y) - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - double precision npar0,nbco0 - external hythik - common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr - common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet - common /hynup1/ bp,x - save /hyipar/,/hyjpar/ - r1=dsqrt(abs(y*y+bp*bp/4.d0+y*bp*dcos(x))) - r2=dsqrt(abs(y*y+bp*bp/4.d0-y*bp*dcos(x))) - s=1.d0-dexp(-0.1d0*sigin*hythik(r2)) - HFUNC2=y*hythik(r1)*s - return - end - -* nuclear overlap function at impact parameter b - double precision function hftaa(c) - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - common /hythic/ BAB(110),TAB(110),TAAB(110) - save /hythic/ - call parinv(c,BAB,TAAB,110,RES) - hftaa=RES - return - end -* - double precision function HFUNC3(x) - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - double precision npar0,nbco0 - external HFUNC4 - common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr - common /hynup1/ bp,xx - save /hyipar/ - if(init.eq.1) then - Rl=Rpb - else - Rl=RA - end if - xx=x - EPS=0.005d0 - A=0.d0 - B=3.d0*Rl - H=0.01d0*(B-A) - CALL SIMPB(A,B,H,EPS,1.d-8,HFUNC4,Y,RES,AIH,AIABS) - HFUNC3=RES - return - end -* - double precision function HFUNC4(y) - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - double precision npar0,nbco0 - external hythik - common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr - common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet - common /hynup1/ bp,x - save /hyipar/,/hyjpar/ - r1=dsqrt(abs(y*y+bp*bp/4.d0+y*bp*dcos(x))) - r2=dsqrt(abs(y*y+bp*bp/4.d0-y*bp*dcos(x))) - HFUNC4=y*hythik(r1)*hythik(r2) - return - end - -* nuclear thickness function - double precision function hythik(r) - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - common /hythic/ BAB(110),TAB(110),TAAB(110) - save /hythic/ - call parinv(r,BAB,TAB,110,RES) - hythik=RES - return - end - -* Wood-Saxon nucleon distrubution - double precision function rhoaa(z) - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - double precision npar0,nbco0 - common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr - common /hygeom/ BC - save /hyipar/,/hygeom/ - if(init.eq.1) then - Rl=Rpb - else - Rl=RA - end if - pi=3.14159d0 - df=0.54d0 - r=sqrt(bc*bc+z*z) - rho0=3.d0/(4.d0*pi*Rl**3)/(1.d0+(pi*df/Rl)**2) - rhoaa=rho0/(1.d0+dexp((r-Rl)/df)) - return - end - -* function to calculate nuclear shadowing factor (for Pb, Au, Pd or Ca beams) - double precision function shad1(kfh,xbjh,Q2h,r) - IMPLICIT DOUBLE PRECISION(A-H, O-Z) - double precision npar0,nbco0,nbcol,npart - external ggshad - common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr - common /hyfpar/ bgen,nbcol,npart,npyt,nhyd - common /hyshad/ bbmin,bbmax,inuc - save /hyipar/,/hyfpar/,/hyshad/ - dimension res(2) - kf=kfh - xbj=xbjh - Q2=Q2h - bb=r - inuc=0 - if(AW.gt.205.d0.and.AW.lt.209.d0) inuc=4 ! Pb-206, 207 or 208 - if(AW.gt.196.d0.and.AW.lt.198.d0) inuc=3 ! Au-197 - if(AW.gt.109.d0.and.AW.lt.111.d0) inuc=2 ! Pd-110 - if(AW.gt.39.d0.and.AW.lt.41.d0) inuc=1 ! Ca-40 - if(inuc.eq.0.and.ipr.eq.1) then - write(6,*) - > 'Warning! Shadowing is not foreseen for atomic weigth A =' - > ,AW - write(6,*)'****************************************************** - >************************' - ipr=0 - end if - if(inuc.ne.0) then - xbj=max(5.d-5,xbj) - xbj=min(0.95d0,xbj) - Q2=max(4.d0,Q2) - Q2=min(520.d0,Q2) - bb=max(0.d0,bb) - call ggshad(inuc,xbj,Q2,bb,res,ta) - if(kf.eq.21) then - shad1=res(1) - elseif(kf.eq.1.or.kf.eq.2.or.kf.eq.3) then - shad1=res(2) - else - shad1=1.d0 - end if - else - shad1=1.d0 - end if - return - end - -****************************************************************************** -* The part of the code which follows below, includes nuclear shadowing model * -* and has been written by Konrad Tywoniuk (Oslo University, Norway) * -****************************************************************************** -c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$ -c$$$ -c$$$ -c$$$ Shadowing from Glauber-Gribov theory for -c$$$ gluons and light quarks (u,d,s and their antiquarks). -c$$$ All Pomeron tree diagrams have been summed (Schwimmer model). -c$$$ More details about the model in -c$$$ K. Tywoniuk, I.C. Arsene, L. Bravina, A. Kaidalov and E. -c$$$ Zabrodin, Phys. Lett. B 657 (2007) 170 -c$$$ -c$$$ We use FIT B from the H1 parameterization published in -c$$$ A. Aktas et al. (H1 Collaboration), Eur. Phys. J C 48 (2006) 715 -c$$$ A. Aktas et al. (H1 Collaboration), Eur. Phys. J C 48 (2006) 749 -c$$$ -c$$$ Main routine is GGSHAD which provides the user with the -c$$$ ratio of nuclear and nucleon parton distribution function -c$$$ normalized by the atomic number for a given -c$$$ INUCL: nucleus (1=Ca,2=Pd,3=Au,4=Pb) -c$$$ X: Bjorken x -c$$$ Q2: q**2, scale squared -c$$$ B: impact parameter -c$$$ -c$$$ Then RES(1): gluon shadowing -c$$$ RES(2): sea quark shadowing -c$$$ TA: nuclear profile function -c$$$ -c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$ - SUBROUTINE GGSHAD(INUCL,X,Q2,B,RES,TAF) - IMPLICIT NONE - DOUBLE PRECISION TA(31,4),IMPAR(31),ANUCL(4) - DOUBLE PRECISION XB(36),Q2V(13) - DOUBLE PRECISION XMAX,XMAXX,Q2MIN,Q2MAX,BMAX - DOUBLE PRECISION G(36,13,4),LQ(36,13,4) - DOUBLE PRECISION TATMP(31),SHAD(2) - DOUBLE PRECISION C(100),D(100),E(100) - DOUBLE PRECISION X,Q2,B - DOUBLE PRECISION RES(2) - DOUBLE PRECISION TAF - DOUBLE PRECISION SEVAL - INTEGER INUCL,TMAX - INTEGER I,IK - - PARAMETER(TMAX=31) - PARAMETER(XMAX=0.1) - PARAMETER(XMAXX=0.95) - PARAMETER(Q2MIN=4.) - PARAMETER(Q2MAX=520.) - PARAMETER(BMAX=30.) - - COMMON/GG07/XB,Q2V,G,LQ - - DATA ANUCL/40.,110.,197.,206./ - DATA IMPAR - > /0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,6.5,7.,7.5,8., - > 8.5,9.,9.5,10.,10.5,11.,11.5,12.,12.5,13.,13.5,14.,14.5,15./ - DATA TA - > /0.0291662,0.0288636,0.0279354,0.0263044,0.0238341,0.0203565, - > 0.0158333,0.0107343,0.00617758,0.00307859,0.00139738, - > 0.000603839,0.000254843,0.000106329,4.40965e-05,1.82206e-05, - > 7.50949e-06,3.08881e-06,1.26837e-06,5.20082e-07,2.12978e-07, - > 8.71154e-08,3.55956e-08,1.45304e-08,5.92622e-09,2.41504e-09, - > 9.83429e-10,4.00184e-10,1.62741e-10,6.61407e-11,2.68657e-11, - > 0.0153487,0.0152775,0.0150614,0.0146926,0.0141561,0.0134268, - > 0.0124651,0.0112115,0.00959258,0.00756907,0.00527514, - > 0.00313482,0.00159852,0.000732479,0.000316628,0.000133116, - > 5.52508e-05,2.27904e-05,9.3692e-06,3.84349e-06,1.57423e-06, - > 6.43962e-07,2.63132e-07,1.07414e-07,4.38089e-08,1.78529e-08, - > 7.2699e-09,2.95832e-09,1.20304e-09,4.8894e-10,1.98603e-10, - > 0.0105299,0.010498,0.0104017,0.010239,0.010006,0.00969677, - > 0.00930205,0.00880781,0.00819264,0.00742455,0.00646058, - > 0.00526252,0.00385835,0.00243988,0.00131252,0.000621079, - > 0.000272322,0.000114993,4.77312e-05,1.96569e-05,8.06377e-06, - > 3.30066e-06,1.34904e-06,5.50751e-07,2.24633e-07,9.15434e-08, - > 3.72777e-08,1.51694e-08,6.16885e-09,2.50714e-09,1.01838e-09, - > 0.0102179,0.0101879,0.0100975,0.00994467,0.00972606, - > 0.00943628,0.0090671,0.00860605,0.00803416,0.00732289, - > 0.00643266,0.00532296,0.00400025,0.00261264,0.00144993, - > 0.00070098,0.000310867,0.000131949,5.48887e-05,2.26247e-05, - > 9.28457e-06,3.80091e-06,1.55358e-06,6.34273e-07,2.58701e-07, - > 1.05427e-07,4.29316e-08,1.74701e-08,7.10447e-09,2.88739e-09, - > 1.17283e-09/ - DATA Q2V - > /4.000, 6.000, 9.000, 13.500, 20.250, 30.375, 45.562, - > 68.344, 102.516, 153.773, 230.660, 345.990, 518.985/ - DATA XB - > /0.99999998E-05, 0.13000000E-04, 0.16899999E-04, - > 0.21970000E-04, 0.28561000E-04, 0.37129299E-04, - > 0.48268099E-04, 0.62748499E-04, 0.81573096E-04, - > 0.10604500E-03, 0.13785800E-03, 0.17921600E-03, - > 0.23298099E-03, 0.30287501E-03, 0.39373801E-03, - > 0.51185902E-03, 0.66541700E-03, 0.86504198E-03, - > 0.11245500E-02, 0.14619200E-02, 0.19005000E-02, - > 0.24706500E-02, 0.32118401E-02, 0.41753901E-02, - > 0.54280101E-02, 0.70564100E-02, 0.91733299E-02, - > 0.11925300E-01, 0.15502900E-01, 0.20153800E-01, - > 0.26200000E-01, 0.34059901E-01, 0.44277899E-01, - > 0.57561301E-01, 0.74829698E-01, 0.97278602E-01/ - DATA G - > /2.542140,2.457280,2.372200,2.283920,2.195320,2.105540, - > 2.014040,1.921480,1.827600,1.733930,1.639990,1.546170, - > 1.453220,1.361040,1.268750,1.178500,1.088730,1.000740, - > 0.914606,0.830842,0.749315,0.671419,0.596883,0.526173, - > 0.460067,0.398489,0.341895,0.289998,0.242628,0.199405, - > 0.159625,0.122521,0.087252,0.053699,0.022804,0.000414, - > 1.825600,1.768440,1.712550,1.654660,1.595830,1.533980, - > 1.472860,1.409910,1.346060,1.282030,1.217520,1.153050, - > 1.089030,1.025930,0.962709,0.900094,0.837787,0.776733, - > 0.716072,0.656570,0.598295,0.541607,0.486661,0.433846, - > 0.383421,0.335585,0.290685,0.248542,0.209305,0.172614, - > 0.138166,0.105552,0.074246,0.044592,0.017935,0.000291, - > 1.445600,1.405920,1.364150,1.322050,1.276900,1.232020, - > 1.184670,1.136620,1.086950,1.037370,0.987225,0.937488, - > 0.888244,0.839061,0.790370,0.741798,0.693426,0.645511, - > 0.598298,0.551516,0.505500,0.460400,0.416326,0.373488, - > 0.332281,0.292614,0.254859,0.218944,0.184976,0.152759, - > 0.122077,0.092747,0.064477,0.037873,0.014527,0.000216, - > 1.209700,1.179300,1.146560,1.112700,1.077790,1.041400, - > 1.003150,0.964598,0.923825,0.882423,0.841304,0.800381, - > 0.759486,0.718947,0.678587,0.638579,0.598671,0.558978, - > 0.519524,0.480617,0.442116,0.404231,0.367057,0.330599, - > 0.295291,0.261100,0.228220,0.196605,0.166367,0.137331, - > 0.109496,0.082677,0.056835,0.032696,0.012013,0.000165, - > 1.049390,1.024810,0.998839,0.971705,0.942633,0.912223, - > 0.880504,0.847236,0.812716,0.777477,0.741670,0.706578, - > 0.671021,0.636183,0.601384,0.566923,0.532196,0.497956, - > 0.463884,0.430069,0.396612,0.363536,0.330915,0.298900, - > 0.267627,0.237155,0.207697,0.179169,0.151566,0.124961, - > 0.099266,0.074436,0.050544,0.028462,0.010013,0.000127, - > 0.943086,0.922388,0.900700,0.877793,0.853134,0.826808, - > 0.798937,0.769284,0.738709,0.706829,0.674818,0.643041, - > 0.611250,0.580168,0.548982,0.517908,0.486720,0.456021, - > 0.425312,0.394899,0.364616,0.334682,0.305173,0.276042, - > 0.247465,0.219584,0.192474,0.166033,0.140391,0.115481, - > 0.091347,0.068001,0.045636,0.025180,0.008510,0.000101, - > 0.858923,0.841883,0.823827,0.803952,0.782124,0.759447, - > 0.734659,0.708141,0.680281,0.651690,0.622372,0.593451, - > 0.564610,0.536045,0.507521,0.479248,0.450845,0.422731, - > 0.394687,0.366774,0.339095,0.311580,0.284350,0.257527, - > 0.231101,0.205195,0.179929,0.155178,0.131077,0.107555, - > 0.084712,0.062629,0.041557,0.022497,0.007326,0.000082, - > 0.790086,0.776082,0.760408,0.743282,0.724441,0.703812, - > 0.681673,0.657893,0.632354,0.605918,0.578976,0.552364, - > 0.525820,0.499665,0.473306,0.447125,0.421027,0.394926, - > 0.369017,0.343250,0.317596,0.292061,0.266805,0.241770, - > 0.217121,0.192830,0.169123,0.145809,0.122984,0.100685, - > 0.078967,0.057983,0.038065,0.020236,0.006364,0.000067, - > 0.732957,0.720881,0.707633,0.692613,0.676200,0.657914, - > 0.637626,0.615846,0.592373,0.567813,0.542877,0.518036, - > 0.493386,0.468899,0.444596,0.420291,0.395908,0.371641, - > 0.347423,0.323358,0.299396,0.275511,0.251830,0.228313, - > 0.205146,0.182274,0.159813,0.137693,0.115983,0.094703, - > 0.073974,0.053979,0.035068,0.018327,0.005580,0.000056, - > 0.684126,0.674389,0.662432,0.649473,0.634783,0.618458, - > 0.600173,0.580108,0.558464,0.535421,0.512051,0.488760, - > 0.465717,0.442837,0.419997,0.397309,0.374393,0.351598, - > 0.328943,0.306305,0.283685,0.261210,0.238875,0.216696, - > 0.194737,0.173035,0.151692,0.130607,0.109864,0.089483, - > 0.069628,0.050489,0.032492,0.016714,0.004940,0.000047, - > 0.642684,0.634012,0.624311,0.612687,0.599751,0.584885, - > 0.568165,0.549716,0.529352,0.507582,0.485673,0.463844, - > 0.442031,0.420480,0.398892,0.377446,0.355834,0.334383, - > 0.312903,0.291520,0.270128,0.248790,0.227632,0.206518, - > 0.185624,0.164939,0.144527,0.124348,0.104415,0.084843, - > 0.065754,0.047403,0.030217,0.015305,0.004396,0.000040, - > 0.606784,0.599501,0.591216,0.581053,0.569262,0.555644, - > 0.540298,0.522979,0.503888,0.483449,0.462725,0.441941, - > 0.421308,0.400965,0.380515,0.360166,0.339667,0.319306, - > 0.298896,0.278531,0.258224,0.237932,0.217694,0.197574, - > 0.177567,0.157731,0.138168,0.118766,0.099587,0.080708, - > 0.062310,0.044656,0.028210,0.014081,0.003935,0.000035, - > 0.575250,0.569055,0.561834,0.552756,0.542358,0.530173, - > 0.515880,0.499697,0.481592,0.462197,0.442485,0.422845, - > 0.403137,0.383704,0.364332,0.344943,0.325467,0.306020, - > 0.286496,0.267106,0.247644,0.228251,0.208881,0.189557, - > 0.170377,0.151316,0.132475,0.113748,0.095226,0.076992, - > 0.059214,0.042195,0.026424,0.013004,0.003541,0.000030, - > 2.510860,2.424890,2.338800,2.251430,2.161030,2.070080, - > 1.976500,1.884380,1.789520,1.695430,1.600900,1.507110, - > 1.413980,1.321430,1.229430,1.138950,1.049380,0.961914, - > 0.875863,0.792403,0.711728,0.634220,0.560088,0.490193, - > 0.424574,0.363647,0.307681,0.256400,0.209654,0.167105, - > 0.128273,0.092911,0.060892,0.033178,0.011804,0.000172, - > 1.797890,1.740890,1.683730,1.625330,1.564220,1.502940, - > 1.439440,1.375560,1.311540,1.247060,1.182260,1.117820, - > 1.053940,0.990598,0.927193,0.864982,0.802868,0.741596, - > 0.681232,0.621931,0.563990,0.507817,0.452982,0.400714, - > 0.350809,0.303514,0.259139,0.217505,0.178849,0.142987, - > 0.109689,0.079022,0.051121,0.027212,0.009208,0.000121, - > 1.421860,1.380980,1.338550,1.294150,1.249120,1.201850, - > 1.153940,1.105280,1.054910,1.004680,0.954677,0.905086, - > 0.855141,0.806467,0.757283,0.708817,0.660443,0.613233, - > 0.565935,0.519307,0.473770,0.428741,0.385002,0.342676, - > 0.301797,0.262591,0.225442,0.189974,0.156683,0.125290, - > 0.095921,0.068630,0.043899,0.022874,0.007411,0.000090, - > 1.187690,1.155530,1.122280,1.087030,1.051020,1.013270, - > 0.973852,0.933953,0.892905,0.851713,0.810455,0.769192, - > 0.728379,0.687798,0.647423,0.607517,0.567475,0.528048, - > 0.488850,0.450119,0.411897,0.374353,0.337466,0.301396, - > 0.266522,0.232713,0.200387,0.169339,0.139725,0.111661, - > 0.085284,0.060612,0.038324,0.019572,0.006096,0.000069, - > 1.028240,1.002610,0.975714,0.947100,0.916605,0.885287, - > 0.852214,0.817988,0.782800,0.747390,0.711939,0.676325, - > 0.641102,0.606200,0.571466,0.536994,0.502546,0.468438, - > 0.434547,0.400992,0.367663,0.334910,0.302578,0.270890, - > 0.240102,0.210041,0.181153,0.153175,0.126406,0.100826, - > 0.076605,0.054074,0.033774,0.016896,0.005055,0.000053, - > 0.923239,0.901357,0.878307,0.853928,0.827994,0.800448, - > 0.770833,0.740864,0.709208,0.677436,0.645505,0.613815, - > 0.582211,0.551121,0.519863,0.488847,0.457919,0.427291, - > 0.396728,0.366593,0.336576,0.306991,0.277689,0.248974, - > 0.220914,0.193485,0.166892,0.141077,0.116297,0.092513, - > 0.069960,0.049024,0.030239,0.014840,0.004277,0.000042, - > 0.840109,0.821572,0.802221,0.780737,0.758007,0.733570, - > 0.707391,0.680060,0.651741,0.622598,0.593525,0.564790, - > 0.535923,0.507607,0.479180,0.450770,0.422703,0.394808, - > 0.366922,0.339201,0.311718,0.284555,0.257727,0.231240, - > 0.205291,0.179821,0.155134,0.131145,0.107877,0.085551, - > 0.064423,0.044815,0.027334,0.013173,0.003668,0.000034, - > 0.771898,0.756403,0.739454,0.720917,0.700462,0.678667, - > 0.655360,0.630295,0.604073,0.577412,0.550894,0.524237, - > 0.497613,0.471552,0.445325,0.419311,0.393354,0.367582, - > 0.341816,0.316292,0.290960,0.265746,0.240859,0.216190, - > 0.191974,0.168255,0.145100,0.122534,0.100715,0.079595, - > 0.059666,0.041209,0.024869,0.011778,0.003176,0.000028, - > 0.715432,0.701697,0.687335,0.671043,0.652774,0.633196, - > 0.611618,0.588703,0.564555,0.539941,0.515046,0.490273, - > 0.465666,0.441417,0.417189,0.393083,0.368842,0.344918, - > 0.320853,0.296990,0.273346,0.249853,0.226554,0.203403, - > 0.180727,0.158303,0.136484,0.115160,0.094448,0.074509, - > 0.055551,0.038122,0.022769,0.010610,0.002775,0.000023, - > 0.666960,0.655627,0.643169,0.628497,0.612317,0.594375, - > 0.574670,0.553372,0.531024,0.507805,0.484742,0.461630, - > 0.438595,0.415762,0.393171,0.370479,0.347803,0.325380, - > 0.302957,0.280541,0.258199,0.236074,0.214141,0.192375, - > 0.170880,0.149710,0.129011,0.108739,0.089084,0.069999, - > 0.052020,0.035461,0.020980,0.009629,0.002450,0.000020, - > 0.626089,0.616427,0.605000,0.592191,0.577470,0.561221, - > 0.543098,0.523316,0.502350,0.480560,0.458749,0.436948, - > 0.415413,0.393933,0.372648,0.351193,0.329853,0.308661, - > 0.287323,0.266214,0.245090,0.224191,0.203371,0.182736, - > 0.162324,0.142202,0.122478,0.103078,0.084237,0.066064, - > 0.048902,0.033114,0.019409,0.008778,0.002175,0.000017, - > 0.590664,0.582220,0.572479,0.560768,0.547507,0.532490, - > 0.515657,0.497052,0.477369,0.456893,0.436199,0.415686, - > 0.395192,0.374819,0.354542,0.334339,0.314120,0.293919, - > 0.273897,0.253708,0.233747,0.213772,0.193977,0.174220, - > 0.154737,0.135524,0.116625,0.098102,0.080014,0.062548, - > 0.046083,0.031041,0.018032,0.008042,0.001942,0.000014, - > 0.559453,0.552473,0.543647,0.533124,0.521070,0.507229, - > 0.491452,0.474143,0.455506,0.435945,0.416256,0.396977, - > 0.377326,0.358052,0.338792,0.319598,0.300262,0.281122, - > 0.261883,0.242706,0.223579,0.204566,0.185577,0.166710, - > 0.148088,0.129570,0.111434,0.093593,0.076187,0.059415, - > 0.043600,0.029214,0.016816,0.007398,0.001744,0.000012, - > 2.494430,2.407750,2.320680,2.230570,2.139410,2.045820, - > 1.954300,1.859160,1.764370,1.669990,1.574490,1.481120, - > 1.386830,1.294720,1.202200,1.111670,1.021960,0.934479, - > 0.848525,0.765125,0.684539,0.607178,0.533502,0.463759, - > 0.398700,0.338188,0.282593,0.231728,0.185473,0.143600, - > 0.105860,0.072262,0.043334,0.020543,0.005838,0.000062, - > 1.781830,1.724140,1.666660,1.605560,1.544230,1.481280, - > 1.418440,1.353130,1.288210,1.223120,1.158330,1.093950, - > 1.029490,0.965787,0.902413,0.840091,0.777784,0.716620, - > 0.656180,0.597156,0.539338,0.483075,0.428748,0.376710, - > 0.327154,0.280125,0.236105,0.194961,0.156820,0.121629, - > 0.089495,0.060685,0.035892,0.016616,0.004508,0.000044, - > 1.407480,1.365110,1.320730,1.276090,1.228930,1.181910, - > 1.132930,1.082570,1.032770,0.982143,0.932084,0.882066, - > 0.832465,0.783250,0.734272,0.685821,0.637505,0.589736, - > 0.542652,0.496248,0.450545,0.405906,0.362451,0.320187, - > 0.279650,0.240757,0.203901,0.169026,0.136216,0.105650, - > 0.077490,0.052176,0.030469,0.013808,0.003599,0.000032, - > 1.173640,1.140320,1.105880,1.069540,1.032320,0.993897, - > 0.953913,0.913021,0.871617,0.829958,0.788517,0.747238, - > 0.706258,0.666094,0.625534,0.585422,0.545512,0.506014, - > 0.466972,0.428374,0.390208,0.352773,0.316057,0.280282, - > 0.245648,0.212195,0.180160,0.149604,0.120645,0.093445, - > 0.068267,0.045631,0.026330,0.011698,0.002940,0.000025, - > 1.014900,0.987942,0.959684,0.930270,0.899044,0.866678, - > 0.832642,0.797620,0.762091,0.726232,0.690987,0.655077, - > 0.620091,0.585191,0.550543,0.515986,0.481432,0.447342, - > 0.413584,0.380034,0.346933,0.314295,0.282265,0.250834, - > 0.220232,0.190535,0.162003,0.134530,0.108400,0.083726, - > 0.060878,0.040604,0.022978,0.010002,0.002422,0.000019, - > 0.909560,0.887017,0.862965,0.837572,0.810211,0.781875, - > 0.751746,0.720666,0.688703,0.656776,0.624847,0.593122, - > 0.561398,0.530391,0.499138,0.468280,0.437229,0.406831, - > 0.376391,0.346350,0.316458,0.286984,0.257979,0.229524, - > 0.201720,0.174606,0.148473,0.123248,0.099111,0.076308, - > 0.055184,0.036277,0.020402,0.008711,0.002038,0.000015, - > 0.827323,0.808240,0.787194,0.764785,0.741102,0.715710, - > 0.688724,0.660472,0.631563,0.602307,0.573344,0.544275, - > 0.515643,0.487203,0.458901,0.430660,0.402505,0.374606, - > 0.346855,0.319408,0.292040,0.265141,0.238553,0.212326, - > 0.186666,0.161641,0.137398,0.113940,0.091458,0.070180, - > 0.050481,0.032923,0.018295,0.007672,0.001739,0.000012, - > 0.759389,0.742705,0.724897,0.705170,0.684041,0.661178, - > 0.636646,0.610897,0.584460,0.557556,0.530715,0.504308, - > 0.477779,0.451618,0.425597,0.399640,0.373702,0.348041, - > 0.322356,0.297020,0.271779,0.246872,0.222152,0.197859, - > 0.173964,0.150642,0.127976,0.105994,0.084903,0.064934, - > 0.046463,0.030071,0.016527,0.006814,0.001499,0.000010, - > 0.703019,0.688814,0.673059,0.655722,0.636535,0.615667, - > 0.593245,0.569605,0.545086,0.520194,0.495411,0.470863, - > 0.446201,0.421921,0.397793,0.373652,0.349528,0.325620, - > 0.301845,0.278188,0.254649,0.231349,0.208318,0.185538, - > 0.163153,0.141239,0.119901,0.099172,0.079281,0.060422, - > 0.043027,0.027647,0.015032,0.006099,0.001305,0.000008, - > 0.655549,0.643053,0.628966,0.613794,0.596440,0.577313, - > 0.556691,0.534653,0.511805,0.488542,0.465448,0.442353, - > 0.419418,0.396807,0.374153,0.351526,0.328993,0.306561, - > 0.284259,0.262079,0.240000,0.218086,0.196361,0.174927, - > 0.153790,0.133102,0.112916,0.093252,0.074414,0.056530, - > 0.040069,0.025572,0.013769,0.005504,0.001148,0.000007, - > 0.614886,0.603957,0.591709,0.577655,0.561834,0.544400, - > 0.525193,0.504891,0.483435,0.461525,0.439729,0.418166, - > 0.396478,0.375064,0.353725,0.332565,0.311308,0.290188, - > 0.269125,0.248171,0.227323,0.206600,0.186070,0.165736, - > 0.145686,0.125993,0.106807,0.088112,0.070133,0.053112, - > 0.037464,0.023751,0.012665,0.004991,0.001015,0.000006, - > 0.579569,0.570010,0.559085,0.546465,0.532125,0.516026, - > 0.498089,0.478853,0.458681,0.438173,0.417404,0.396958, - > 0.376499,0.356301,0.336115,0.316014,0.295930,0.275886, - > 0.255871,0.236010,0.216223,0.196524,0.176991,0.157626, - > 0.138515,0.119742,0.101404,0.083545,0.066339,0.050078, - > 0.035167,0.022152,0.011703,0.004550,0.000904,0.000005, - > 0.548753,0.540580,0.530520,0.519276,0.506004,0.491061, - > 0.474436,0.456132,0.436996,0.417488,0.397755,0.378421, - > 0.359085,0.339851,0.320700,0.301518,0.282357,0.263296, - > 0.244273,0.225311,0.206421,0.187624,0.168959,0.150424, - > 0.132145,0.114172,0.096557,0.079446,0.062967,0.047378, - > 0.033126,0.020736,0.010987,0.004166,0.000809,0.000004, - > 2.486670,2.401390,2.313150,2.223100,2.132210,2.040620, - > 1.946240,1.852960,1.758660,1.663790,1.569300,1.475400, - > 1.381630,1.288980,1.197380,1.106900,1.017420,0.930068, - > 0.844430,0.760955,0.680549,0.603582,0.530075,0.460570, - > 0.395641,0.335338,0.279919,0.229212,0.183101,0.141385, - > 0.103826,0.070754,0.041898,0.019580,0.005430,0.000056, - > 1.776150,1.719900,1.661320,1.600440,1.539110,1.477010, - > 1.412790,1.348550,1.283220,1.218900,1.153830,1.088770, - > 1.025080,0.961812,0.898582,0.836002,0.773940,0.712963, - > 0.652650,0.593610,0.536099,0.480001,0.425916,0.373940, - > 0.324454,0.277574,0.233716,0.192654,0.154624,0.119622, - > 0.087676,0.059112,0.034659,0.015816,0.004189,0.000039, - > 1.403340,1.361000,1.316930,1.271610,1.225380,1.177160, - > 1.128240,1.078580,1.028640,0.977947,0.928063,0.878521, - > 0.828550,0.779706,0.730569,0.682199,0.633948,0.586522, - > 0.539401,0.493171,0.447639,0.403145,0.359673,0.317729, - > 0.277222,0.238454,0.201716,0.166913,0.134251,0.103840, - > 0.075864,0.050772,0.029388,0.013129,0.003341,0.000029, - > 1.169950,1.136380,1.102100,1.065820,1.028190,0.989956, - > 0.950189,0.908974,0.867893,0.826368,0.784940,0.743815, - > 0.702872,0.662550,0.622187,0.582284,0.542431,0.503183, - > 0.464112,0.425583,0.387521,0.350185,0.313576,0.277943, - > 0.243398,0.210036,0.178125,0.147636,0.118821,0.091774, - > 0.066794,0.044372,0.025371,0.011112,0.002727,0.000022, - > 1.011370,0.984724,0.956435,0.926814,0.895443,0.862964, - > 0.829274,0.794479,0.758875,0.722999,0.687335,0.652160, - > 0.616882,0.582106,0.547373,0.512881,0.478596,0.444627, - > 0.410813,0.377446,0.344502,0.311951,0.279953,0.248594, - > 0.218124,0.188529,0.160091,0.132732,0.106696,0.082174, - > 0.059499,0.039204,0.022122,0.009492,0.002245,0.000017, - > 0.906689,0.883943,0.860162,0.834079,0.807067,0.778589, - > 0.748446,0.717460,0.685716,0.653595,0.621726,0.590095, - > 0.558631,0.527321,0.496368,0.465511,0.434680,0.404249, - > 0.373821,0.343908,0.314123,0.284707,0.255819,0.227414, - > 0.199673,0.172714,0.146637,0.121510,0.097519,0.074849, - > 0.053902,0.035215,0.019623,0.008259,0.001888,0.000014, - > 0.824084,0.805081,0.784198,0.761980,0.737877,0.712401, - > 0.685185,0.657419,0.628573,0.599438,0.570456,0.541554, - > 0.512854,0.484478,0.456174,0.428019,0.399930,0.372215, - > 0.344523,0.317092,0.289826,0.262998,0.236396,0.210320, - > 0.184715,0.159801,0.135667,0.112306,0.089938,0.068791, - > 0.049274,0.031936,0.017585,0.007270,0.001610,0.000011, - > 0.756740,0.739909,0.722048,0.702184,0.681164,0.658293, - > 0.633526,0.608140,0.581469,0.554572,0.527904,0.501517, - > 0.475145,0.448970,0.422970,0.397121,0.371248,0.345564, - > 0.320059,0.294788,0.269623,0.244755,0.220148,0.195923, - > 0.172134,0.148871,0.126312,0.104417,0.083450,0.063624, - > 0.045337,0.029156,0.015873,0.006450,0.001387,0.000009, - > 0.700680,0.686067,0.670562,0.652740,0.633579,0.612848, - > 0.590476,0.566706,0.542207,0.517556,0.492696,0.468105, - > 0.443562,0.419350,0.395237,0.371204,0.347185,0.323339, - > 0.299543,0.276011,0.252614,0.229355,0.206336,0.183677, - > 0.161386,0.139527,0.118283,0.097686,0.077903,0.059183, - > 0.041956,0.026789,0.014429,0.005771,0.001207,0.000008, - > 0.652991,0.640812,0.626563,0.610980,0.593505,0.574586, - > 0.553939,0.531965,0.509196,0.485938,0.462796,0.439793, - > 0.416876,0.394269,0.371716,0.349126,0.326720,0.304355, - > 0.282067,0.259966,0.237971,0.216140,0.194490,0.173138, - > 0.152089,0.131483,0.111379,0.091839,0.073081,0.055347, - > 0.039055,0.024764,0.013208,0.005204,0.001061,0.000006, - > 0.612493,0.601527,0.589260,0.574987,0.559293,0.541793, - > 0.522714,0.502071,0.480580,0.459021,0.437202,0.415483, - > 0.394052,0.372703,0.351511,0.330271,0.309080,0.288012, - > 0.267062,0.246156,0.225353,0.204703,0.184235,0.163954, - > 0.144013,0.124421,0.105317,0.086709,0.068861,0.051973, - > 0.036498,0.022991,0.012143,0.004716,0.000938,0.000005, - > 0.577224,0.567718,0.556768,0.543992,0.529614,0.513264, - > 0.495593,0.476332,0.456118,0.435643,0.414958,0.394475, - > 0.374159,0.353948,0.333917,0.313860,0.293761,0.273842, - > 0.253865,0.234066,0.214285,0.194711,0.175202,0.155932, - > 0.136903,0.118198,0.099962,0.082192,0.065107,0.048988, - > 0.034247,0.021430,0.011215,0.004297,0.000835,0.000005, - > 0.546580,0.538342,0.528410,0.516932,0.503593,0.488594, - > 0.471844,0.453700,0.434588,0.415081,0.395496,0.376054, - > 0.356702,0.337597,0.318488,0.299383,0.280225,0.261271, - > 0.242274,0.223404,0.204576,0.185839,0.167243,0.148787, - > 0.130586,0.112676,0.095184,0.078140,0.061781,0.046329, - > 0.032244,0.020050,0.010400,0.003932,0.000747,0.000004/ - DATA LQ - > /0.514104,0.515832,0.517091,0.517536,0.517832,0.517221, - > 0.516574,0.515224,0.513287,0.510858,0.508000,0.504917, - > 0.501246,0.497403,0.492983,0.488152,0.482358,0.476122, - > 0.468649,0.460339,0.450645,0.439812,0.427166,0.412836, - > 0.396439,0.377753,0.356388,0.331566,0.302870,0.269780, - > 0.231439,0.187668,0.138736,0.086850,0.036211,0.000604, - > 0.579741,0.577476,0.575013,0.571386,0.567464,0.562828, - > 0.557414,0.551410,0.544332,0.536891,0.528953,0.521049, - > 0.512612,0.503967,0.494945,0.485619,0.475531,0.464958, - > 0.453985,0.442176,0.429817,0.416392,0.402026,0.386267, - > 0.369251,0.350413,0.329631,0.306106,0.279316,0.248539, - > 0.213030,0.172492,0.127054,0.078865,0.032177,0.000513, - > 0.617205,0.612889,0.607448,0.601280,0.593987,0.585725, - > 0.577239,0.567645,0.557436,0.546468,0.535124,0.523690, - > 0.511913,0.500192,0.488271,0.476132,0.463482,0.450653, - > 0.437460,0.423920,0.409767,0.395178,0.379804,0.363620, - > 0.346492,0.327958,0.307882,0.285606,0.260283,0.231441, - > 0.198228,0.160224,0.117620,0.072457,0.029003,0.000445, - > 0.635918,0.629589,0.622092,0.614001,0.604751,0.594717, - > 0.583683,0.572041,0.559188,0.546115,0.532523,0.518899, - > 0.505201,0.491782,0.477675,0.463752,0.449738,0.435623, - > 0.421062,0.406484,0.391587,0.376279,0.360631,0.344349, - > 0.327299,0.309280,0.289970,0.268596,0.244665,0.217365, - > 0.185991,0.150104,0.109826,0.067166,0.026423,0.000392, - > 0.642211,0.634584,0.625818,0.616705,0.605951,0.594718, - > 0.582351,0.569037,0.554986,0.540276,0.525303,0.510521, - > 0.495596,0.480637,0.465827,0.451058,0.435872,0.421066, - > 0.405880,0.390796,0.375397,0.359903,0.344185,0.327987, - > 0.311213,0.293614,0.274949,0.254454,0.231575,0.205516, - > 0.175652,0.141477,0.103138,0.062618,0.024219,0.000348, - > 0.639180,0.631092,0.621893,0.612137,0.600935,0.589006, - > 0.576122,0.561754,0.547048,0.531226,0.515512,0.499878, - > 0.484089,0.468662,0.453258,0.437735,0.422292,0.407026, - > 0.391544,0.376258,0.360735,0.345208,0.329470,0.313579, - > 0.297158,0.280066,0.261996,0.242200,0.220305,0.195335, - > 0.166765,0.134023,0.097367,0.058717,0.022356,0.000312, - > 0.631156,0.622788,0.613580,0.603527,0.592277,0.579779, - > 0.566484,0.551813,0.536423,0.520508,0.504251,0.488021, - > 0.472093,0.456329,0.440551,0.424908,0.409169,0.393745, - > 0.378229,0.362803,0.347473,0.331992,0.316442,0.300771, - > 0.284762,0.268174,0.250682,0.231645,0.210500,0.186498, - > 0.158976,0.127569,0.092372,0.055348,0.020770,0.000282, - > 0.619698,0.611519,0.602336,0.591943,0.580814,0.568514, - > 0.554982,0.540509,0.524970,0.508623,0.492412,0.476046, - > 0.459935,0.444036,0.428166,0.412462,0.396815,0.381305, - > 0.365892,0.350546,0.335395,0.320121,0.304877,0.289437, - > 0.273841,0.257707,0.240774,0.222305,0.201894,0.178726, - > 0.152209,0.121898,0.087984,0.052399,0.019398,0.000257, - > 0.606343,0.598499,0.589365,0.579821,0.568434,0.556327, - > 0.542872,0.528454,0.512900,0.496651,0.480390,0.464149, - > 0.447972,0.432310,0.416261,0.400748,0.385116,0.369727, - > 0.354542,0.339434,0.324280,0.309328,0.294433,0.279362, - > 0.264086,0.248363,0.231853,0.214025,0.194261,0.171833, - > 0.146197,0.116876,0.084113,0.049805,0.018206,0.000236, - > 0.591896,0.584256,0.575797,0.566248,0.555494,0.543797, - > 0.530664,0.516396,0.501105,0.484951,0.468541,0.452368, - > 0.436358,0.420764,0.405075,0.389617,0.374238,0.359020, - > 0.343972,0.329043,0.314268,0.299607,0.284884,0.270198, - > 0.255337,0.239949,0.223985,0.206644,0.187466,0.165701, - > 0.140843,0.112412,0.080664,0.047512,0.017165,0.000218, - > 0.576697,0.569585,0.561998,0.552972,0.542790,0.531174, - > 0.518335,0.504287,0.489160,0.473285,0.457139,0.441230, - > 0.425520,0.409903,0.394534,0.379177,0.364141,0.349046, - > 0.334214,0.319602,0.305086,0.290690,0.276241,0.261890, - > 0.247381,0.232411,0.216817,0.199969,0.181340,0.160182, - > 0.135997,0.108375,0.077551,0.045446,0.016236,0.000202, - > 0.561779,0.555491,0.548244,0.539719,0.529917,0.518731, - > 0.506334,0.492791,0.477798,0.462049,0.446208,0.430534, - > 0.414938,0.399619,0.384500,0.369572,0.354611,0.339912, - > 0.325296,0.310898,0.296649,0.282465,0.268439,0.254375, - > 0.240158,0.225578,0.210376,0.193910,0.175777,0.155185, - > 0.131618,0.104727,0.074745,0.043586,0.015406,0.000188, - > 0.547406,0.541522,0.534772,0.526637,0.517320,0.506665, - > 0.494929,0.481541,0.466958,0.451581,0.435977,0.420382, - > 0.405163,0.390039,0.375242,0.360562,0.345699,0.331352, - > 0.316934,0.302875,0.288897,0.275075,0.261211,0.247494, - > 0.233546,0.219335,0.204477,0.188465,0.170694,0.150613, - > 0.127623,0.101395,0.072188,0.041898,0.014659,0.000176, - > 0.510949,0.512461,0.513125,0.513775,0.513627,0.512953, - > 0.511676,0.510153,0.508026,0.505586,0.502488,0.499514, - > 0.495918,0.491777,0.487221,0.482318,0.476322,0.469865, - > 0.462254,0.453539,0.443579,0.431957,0.418405,0.402714, - > 0.384513,0.363211,0.338441,0.309525,0.276023,0.237418, - > 0.194095,0.146991,0.098692,0.053921,0.018675,0.000251, - > 0.573840,0.571539,0.568660,0.564636,0.560125,0.555429, - > 0.549250,0.542654,0.535827,0.528339,0.520207,0.512007, - > 0.503505,0.494854,0.485721,0.476289,0.466253,0.455649, - > 0.444534,0.432471,0.419685,0.405709,0.390526,0.373656, - > 0.355038,0.334208,0.310578,0.283443,0.252376,0.217005, - > 0.177144,0.133976,0.089700,0.048645,0.016531,0.000213, - > 0.610274,0.604837,0.598760,0.591967,0.584049,0.575825, - > 0.566634,0.556589,0.545957,0.534847,0.523533,0.511998, - > 0.500328,0.488702,0.476646,0.464532,0.451851,0.438900, - > 0.425698,0.412022,0.397498,0.382650,0.366763,0.349574, - > 0.330976,0.310680,0.288132,0.262614,0.233680,0.200673, - > 0.163728,0.123628,0.082475,0.044455,0.014854,0.000185, - > 0.627192,0.620322,0.611968,0.603362,0.593426,0.582683, - > 0.571079,0.558798,0.545824,0.532347,0.519069,0.505318, - > 0.491597,0.478201,0.464023,0.450314,0.436183,0.422097, - > 0.407576,0.392963,0.377856,0.362363,0.346096,0.329071, - > 0.310942,0.291304,0.269699,0.245618,0.218210,0.187279, - > 0.152695,0.115069,0.076563,0.041003,0.013496,0.000163, - > 0.631625,0.623836,0.614536,0.604467,0.593362,0.581189, - > 0.568126,0.554473,0.540005,0.525342,0.510320,0.495556, - > 0.480543,0.465856,0.450901,0.436047,0.421171,0.406166, - > 0.391154,0.375993,0.360616,0.344868,0.328695,0.311826, - > 0.294125,0.275057,0.254485,0.231491,0.205428,0.176164, - > 0.143362,0.107930,0.071534,0.038049,0.012338,0.000145, - > 0.628674,0.619289,0.609887,0.598892,0.587603,0.574461, - > 0.560634,0.545957,0.530648,0.515019,0.499229,0.483499, - > 0.467857,0.452587,0.437187,0.421817,0.406489,0.391076, - > 0.375987,0.360550,0.345003,0.329399,0.313389,0.296812, - > 0.279488,0.261191,0.241276,0.219209,0.194576,0.166551, - > 0.135508,0.101696,0.067217,0.035540,0.011363,0.000130, - > 0.619655,0.610842,0.600741,0.589506,0.577444,0.564098, - > 0.550163,0.534861,0.519180,0.503427,0.486955,0.471011, - > 0.455040,0.439305,0.423732,0.408132,0.392629,0.377169, - > 0.361831,0.346498,0.331073,0.315724,0.299859,0.283638, - > 0.266810,0.249063,0.229985,0.208716,0.185019,0.158309, - > 0.128580,0.096355,0.063475,0.033376,0.010535,0.000117, - > 0.607474,0.598863,0.589070,0.578071,0.565959,0.552353, - > 0.538238,0.523141,0.507028,0.491211,0.474478,0.458296, - > 0.442260,0.426619,0.410731,0.395250,0.379629,0.364356, - > 0.348971,0.333841,0.318608,0.303265,0.287813,0.272025, - > 0.255702,0.238406,0.219910,0.199549,0.176755,0.151092, - > 0.122529,0.091666,0.060228,0.031496,0.009822,0.000107, - > 0.593760,0.585285,0.576038,0.564716,0.553007,0.539720, - > 0.525715,0.510544,0.494716,0.478059,0.461888,0.445882, - > 0.429808,0.413922,0.398545,0.382937,0.367557,0.352300, - > 0.337238,0.322197,0.307168,0.292182,0.277095,0.261674, - > 0.245761,0.228939,0.211113,0.191442,0.169413,0.144749, - > 0.117225,0.087581,0.057361,0.029844,0.009202,0.000098, - > 0.579324,0.570883,0.561883,0.551297,0.539642,0.526840, - > 0.512914,0.497903,0.481998,0.465838,0.449788,0.433693, - > 0.417703,0.402286,0.386792,0.371437,0.356215,0.341283, - > 0.326313,0.311637,0.296947,0.282140,0.267321,0.252270, - > 0.236761,0.220580,0.203236,0.184230,0.162954,0.139137, - > 0.112557,0.083990,0.054824,0.028391,0.008663,0.000090, - > 0.564315,0.556562,0.547782,0.537905,0.526572,0.514011, - > 0.500033,0.485378,0.470013,0.453915,0.438022,0.422153, - > 0.406530,0.391112,0.375843,0.360792,0.345939,0.331089, - > 0.316383,0.301913,0.287442,0.273075,0.258564,0.243883, - > 0.228827,0.213013,0.196181,0.177712,0.157167,0.134036, - > 0.108376,0.080670,0.052555,0.027085,0.008183,0.000084, - > 0.549474,0.542540,0.534001,0.524222,0.513538,0.501429, - > 0.487998,0.473665,0.458530,0.442760,0.426794,0.411187, - > 0.395839,0.380807,0.365652,0.350917,0.336154,0.321565, - > 0.307328,0.293072,0.278886,0.264801,0.250668,0.236339, - > 0.221517,0.206147,0.189835,0.171902,0.151860,0.129450, - > 0.104545,0.077707,0.050505,0.025913,0.007753,0.000078, - > 0.535066,0.528539,0.520546,0.511350,0.500893,0.489256, - > 0.476407,0.462092,0.447135,0.431648,0.416302,0.400932, - > 0.385885,0.370871,0.356128,0.341562,0.327354,0.312933, - > 0.298890,0.284864,0.270948,0.257177,0.243326,0.229321, - > 0.214967,0.199986,0.183994,0.166574,0.147078,0.125313, - > 0.101089,0.075065,0.048640,0.024855,0.007368,0.000073, - > 0.509409,0.510722,0.511674,0.511738,0.511286,0.510766, - > 0.509198,0.507564,0.505152,0.502704,0.499650,0.496452, - > 0.492756,0.488747,0.483986,0.478817,0.472835,0.466164, - > 0.458412,0.449362,0.438892,0.426796,0.412430,0.395391, - > 0.375497,0.352110,0.324651,0.292275,0.254919,0.212461, - > 0.165886,0.117521,0.071431,0.033526,0.009193,0.000091, - > 0.570888,0.568118,0.564921,0.561074,0.555970,0.550478, - > 0.544438,0.537826,0.530394,0.523005,0.514685,0.506342, - > 0.497778,0.489102,0.479852,0.470348,0.460202,0.449558, - > 0.438121,0.426046,0.412891,0.398444,0.382493,0.364765, - > 0.344689,0.321925,0.295886,0.265915,0.231610,0.192817, - > 0.150406,0.106352,0.064426,0.030028,0.008099,0.000077, - > 0.605452,0.599834,0.593296,0.586627,0.578039,0.569200, - > 0.559754,0.549624,0.538672,0.527535,0.516054,0.504396, - > 0.492657,0.481104,0.468974,0.456682,0.443949,0.431100, - > 0.417683,0.403750,0.389161,0.373852,0.357151,0.339250, - > 0.319633,0.297684,0.273111,0.245042,0.213167,0.177272, - > 0.138115,0.097498,0.058870,0.027271,0.007249,0.000067, - > 0.621534,0.614157,0.605835,0.596272,0.585946,0.574335, - > 0.562822,0.550233,0.536841,0.523598,0.509920,0.496293, - > 0.482476,0.468781,0.455150,0.441275,0.427053,0.412946, - > 0.398426,0.383702,0.368436,0.352438,0.335827,0.317985, - > 0.298878,0.277758,0.254411,0.228040,0.198146,0.164659, - > 0.128112,0.090261,0.054331,0.025017,0.006563,0.000059, - > 0.625646,0.616575,0.607120,0.596395,0.584780,0.572447, - > 0.558731,0.544391,0.529948,0.515041,0.499832,0.485344, - > 0.470145,0.455555,0.440783,0.425851,0.410922,0.396061, - > 0.381038,0.365825,0.350234,0.334195,0.317648,0.300202, - > 0.281570,0.261329,0.238985,0.213898,0.185675,0.154079, - > 0.119726,0.084159,0.050488,0.023108,0.005981,0.000052, - > 0.621816,0.611958,0.601503,0.590515,0.577713,0.564667, - > 0.550311,0.535181,0.519529,0.503933,0.488235,0.472641, - > 0.456837,0.441572,0.426210,0.410935,0.395680,0.380387, - > 0.365098,0.349710,0.334114,0.318243,0.301798,0.284744, - > 0.266642,0.247149,0.225778,0.201848,0.175065,0.145114, - > 0.112575,0.078976,0.047163,0.021485,0.005492,0.000047, - > 0.612456,0.602615,0.592112,0.580629,0.567547,0.554000, - > 0.539266,0.523505,0.507494,0.491317,0.475170,0.459335, - > 0.443356,0.427588,0.412089,0.396613,0.381018,0.365836, - > 0.350418,0.335168,0.319712,0.303995,0.287942,0.271269, - > 0.253737,0.234846,0.214374,0.191514,0.165911,0.137420, - > 0.106436,0.074522,0.044401,0.020096,0.005079,0.000042, - > 0.599906,0.590658,0.579871,0.568198,0.555552,0.541752, - > 0.526573,0.511088,0.494690,0.478321,0.462067,0.446064, - > 0.430078,0.414174,0.398700,0.383180,0.367731,0.352392, - > 0.337125,0.322012,0.306729,0.291406,0.275582,0.259369, - > 0.242384,0.224158,0.204420,0.182500,0.157984,0.130669, - > 0.101087,0.070603,0.041960,0.018893,0.004724,0.000039, - > 0.586385,0.576832,0.566175,0.554843,0.542371,0.528570, - > 0.513481,0.497995,0.481584,0.465454,0.449166,0.432955, - > 0.417165,0.401569,0.386012,0.370517,0.355274,0.340119, - > 0.324992,0.310107,0.295066,0.280058,0.264679,0.248876, - > 0.232336,0.214712,0.195681,0.174542,0.150960,0.124776, - > 0.096383,0.067222,0.039821,0.017842,0.004416,0.000035, - > 0.571562,0.562230,0.553483,0.541106,0.528820,0.515154, - > 0.500724,0.484959,0.468895,0.452796,0.436478,0.420645, - > 0.404877,0.389321,0.373971,0.358762,0.343687,0.328774, - > 0.313968,0.299216,0.284619,0.269828,0.254780,0.239398, - > 0.223345,0.206321,0.187886,0.167537,0.144752,0.119523, - > 0.092229,0.064212,0.037929,0.016920,0.004149,0.000033, - > 0.556430,0.547716,0.538118,0.527507,0.515330,0.502167, - > 0.487523,0.472378,0.456637,0.440513,0.424536,0.408900, - > 0.393308,0.378059,0.362874,0.347899,0.333001,0.318366, - > 0.303866,0.289409,0.275001,0.260570,0.245889,0.230899, - > 0.215327,0.198784,0.180938,0.161191,0.139240,0.114867, - > 0.088515,0.061527,0.036242,0.016094,0.003912,0.000030, - > 0.541088,0.533114,0.524472,0.513843,0.502197,0.489243, - > 0.475392,0.460194,0.444773,0.428968,0.413337,0.397800, - > 0.382488,0.367387,0.352568,0.337799,0.323189,0.308771, - > 0.294545,0.280383,0.266247,0.252107,0.237863,0.223231, - > 0.208050,0.191998,0.174706,0.155516,0.134237,0.110630, - > 0.085174,0.059101,0.034727,0.015356,0.003700,0.000028, - > 0.526746,0.519306,0.510794,0.500680,0.489449,0.477026, - > 0.463177,0.448622,0.433350,0.417867,0.402498,0.387171, - > 0.372235,0.357466,0.342866,0.328463,0.314082,0.299986, - > 0.285991,0.272106,0.258305,0.244422,0.230511,0.216246, - > 0.201466,0.185822,0.169011,0.150402,0.129717,0.106819, - > 0.082146,0.056897,0.033348,0.014690,0.003511,0.000026, - > 0.508294,0.509537,0.510232,0.510364,0.510028,0.509066, - > 0.507574,0.506096,0.503764,0.501269,0.498261,0.494940, - > 0.491357,0.487139,0.482585,0.477518,0.471487,0.464685, - > 0.456954,0.448138,0.437547,0.425359,0.410953,0.393893, - > 0.373877,0.350282,0.322653,0.290145,0.252556,0.209889, - > 0.163221,0.114885,0.069163,0.031966,0.008546,0.000081, - > 0.569358,0.566803,0.562983,0.559001,0.554723,0.548955, - > 0.542742,0.536150,0.528626,0.520981,0.513074,0.504898, - > 0.496364,0.487653,0.478428,0.468815,0.458593,0.447996, - > 0.436563,0.424615,0.411428,0.396895,0.381014,0.363131, - > 0.342953,0.320173,0.294051,0.263852,0.229272,0.190367, - > 0.147875,0.103902,0.062322,0.028615,0.007526,0.000069, - > 0.603670,0.598438,0.591834,0.584727,0.576222,0.567680, - > 0.558196,0.547591,0.536787,0.525654,0.514148,0.502844, - > 0.491102,0.479283,0.467390,0.455046,0.442422,0.429358, - > 0.416062,0.402212,0.387705,0.372222,0.355668,0.337678, - > 0.318005,0.295922,0.271145,0.242997,0.210936,0.174968, - > 0.135724,0.095201,0.057278,0.025971,0.006733,0.000060, - > 0.619938,0.612660,0.603237,0.594399,0.584220,0.572945, - > 0.560906,0.547974,0.535000,0.521640,0.508074,0.494503, - > 0.480788,0.467052,0.453309,0.439430,0.425326,0.411154, - > 0.396834,0.381957,0.366731,0.350903,0.334155,0.316435, - > 0.297197,0.275976,0.252539,0.225982,0.196020,0.162415, - > 0.125838,0.088092,0.052512,0.023816,0.006094,0.000053, - > 0.623752,0.615162,0.604853,0.594199,0.582617,0.570319, - > 0.556672,0.542512,0.527993,0.513108,0.498185,0.483234, - > 0.468378,0.453629,0.438899,0.424139,0.409158,0.394331, - > 0.379329,0.364110,0.348592,0.332662,0.316045,0.298521, - > 0.279860,0.259548,0.237136,0.211942,0.183600,0.151947, - > 0.117534,0.082104,0.048767,0.021982,0.005551,0.000047, - > 0.619727,0.610281,0.600069,0.588100,0.575830,0.562394, - > 0.548038,0.533038,0.517667,0.501730,0.486231,0.470564, - > 0.454906,0.439606,0.424254,0.409150,0.393803,0.378587, - > 0.363406,0.347957,0.332413,0.316594,0.300201,0.283135, - > 0.264987,0.245411,0.223977,0.199987,0.173072,0.143020, - > 0.110486,0.077018,0.045586,0.020431,0.005096,0.000042, - > 0.610495,0.600588,0.589924,0.578181,0.565270,0.551795, - > 0.537196,0.521624,0.505437,0.489150,0.472983,0.457321, - > 0.441422,0.425796,0.410213,0.394710,0.379352,0.363971, - > 0.348802,0.333422,0.317927,0.302300,0.286286,0.269619, - > 0.252021,0.233168,0.212605,0.189631,0.163957,0.135368, - > 0.104416,0.072632,0.042854,0.019101,0.004711,0.000038, - > 0.597856,0.588522,0.577800,0.566289,0.553520,0.539436, - > 0.524435,0.508812,0.492453,0.476231,0.459950,0.443865, - > 0.427942,0.412263,0.396763,0.381179,0.365851,0.350734, - > 0.335459,0.320164,0.305048,0.289712,0.274016,0.257736, - > 0.240735,0.222511,0.202688,0.180671,0.156046,0.128701, - > 0.099140,0.068812,0.040482,0.017948,0.004381,0.000035, - > 0.584084,0.574752,0.564446,0.552705,0.540225,0.526223, - > 0.511434,0.495863,0.479694,0.463424,0.447226,0.431058, - > 0.415167,0.399560,0.384026,0.368648,0.353514,0.338375, - > 0.323406,0.308441,0.293400,0.278400,0.263020,0.247257, - > 0.230713,0.213051,0.193943,0.172756,0.149112,0.122832, - > 0.094502,0.065475,0.038404,0.016946,0.004094,0.000032, - > 0.569517,0.560385,0.550236,0.539113,0.526421,0.513222, - > 0.498378,0.482790,0.466696,0.450744,0.434441,0.418603, - > 0.402901,0.387328,0.372129,0.356950,0.341855,0.326972, - > 0.312231,0.297539,0.282901,0.268148,0.253147,0.237801, - > 0.221735,0.204645,0.186195,0.165769,0.142932,0.117672, - > 0.090399,0.062530,0.036573,0.016065,0.003846,0.000029, - > 0.554361,0.545724,0.536138,0.525127,0.513348,0.500011, - > 0.485538,0.470383,0.454424,0.438520,0.422677,0.406888, - > 0.391373,0.376154,0.360890,0.345988,0.331296,0.316666, - > 0.302083,0.287727,0.273311,0.258904,0.244293,0.229265, - > 0.213676,0.197142,0.179286,0.159481,0.137449,0.113024, - > 0.086731,0.059884,0.034929,0.015279,0.003625,0.000027, - > 0.539566,0.531375,0.522021,0.511690,0.499982,0.487119, - > 0.474212,0.458150,0.442576,0.426822,0.411177,0.395600, - > 0.380520,0.365481,0.350710,0.336033,0.321423,0.307074, - > 0.292885,0.278746,0.264660,0.250555,0.236263,0.221649, - > 0.206449,0.190350,0.173011,0.154255,0.132811,0.108865, - > 0.083444,0.057506,0.033461,0.014573,0.003428,0.000025, - > 0.524865,0.517319,0.508721,0.498692,0.487348,0.474698, - > 0.461245,0.446504,0.431299,0.416072,0.400507,0.385278, - > 0.370252,0.355605,0.341042,0.326661,0.312354,0.298301, - > 0.284289,0.270473,0.256711,0.242903,0.228985,0.214684, - > 0.199885,0.184233,0.167367,0.148740,0.128003,0.105093, - > 0.080442,0.055358,0.032128,0.013936,0.003252,0.000024/ - -C IF (X.GE.XMAXX.OR.Q2.LT.Q2MIN.OR.Q2.GT.Q2MAX) THEN -C WRITE(*,*) 'X or Q2 out of range!!' -C WRITE(*,*) 'Valid range is:' -C WRITE(*,*) '1E-05 < X < 0.95 and' -C WRITE(*,*) '4 < Q2 < 520 GeV**2' -C RETURN -C END IF - - IF (X.GE.XMAX.OR.B.GT.BMAX) THEN - RES(1) = 1.0 - RES(2) = 1.0 - RETURN - END IF - IF(Q2.LT.Q2MIN) THEN - Q2 = Q2MIN - ELSE IF(Q2.GT.Q2MAX) THEN - Q2 = Q2MAX - ENDIF - - CALL GGINTER(INUCL,X,Q2,SHAD) - - DO I=1,31 - TATMP(I) = TA(I,INUCL) - ENDDO - CALL SPLINE(TMAX,IMPAR,TATMP,C,D,E) - TAF = SEVAL(TMAX,B,IMPAR,TATMP,C,D,E) - - DO IK=1,2 - RES(IK) = 1.0/(1.0 + (ANUCL(INUCL)-1.0)*SHAD(IK)*TAF) - ENDDO - - RETURN - END - - SUBROUTINE GGINTER(INUCL,X,Q2,SHAD) - IMPLICIT DOUBLE PRECISION(A-H,L-Z) - IMPLICIT INTEGER(I-K) - - DIMENSION SHAD(2) - DIMENSION G(36,13,4),LQ(36,13,4) - DIMENSION XB(36),Q2V(13) - DIMENSION GQTMP(13),LQQTMP(13) - DIMENSION GXTMP(36),LQXTMP(36) - DIMENSION C(100),D(100),E(100) - - PARAMETER(INMAX=13,IMMAX=36) - - COMMON/GG07/ XB,Q2V,G,LQ - - DO K=1,36 - DO J=1,13 - GQTMP(J) = G(K,J,INUCL) - LQQTMP(J) = LQ(K,J,INUCL) - ENDDO - CALL SPLINE(INMAX,Q2V,GQTMP,C,D,E) - GXTMP(K) = SEVAL(INMAX,Q2,Q2V,GQTMP,C,D,E) - CALL SPLINE(INMAX,Q2V,LQQTMP,C,D,E) - LQXTMP(K) = SEVAL(INMAX,Q2,Q2V,LQQTMP,C,D,E) - ENDDO - - CALL SPLINE(IMMAX,XB,GXTMP,C,D,E) - SHAD(1) = SEVAL(IMMAX,X,XB,GXTMP,C,D,E) - CALL SPLINE(IMMAX,XB,LQXTMP,C,D,E) - SHAD(2) = SEVAL(IMMAX,X,XB,LQXTMP,C,D,E) - - RETURN - END - -C --------------------------------------------------------------------- - SUBROUTINE SPLINE(N,X,Y,B,C,D) -C --------------------------------------------------------------------- -c*************************************************************************** -C CALCULATE THE COEFFICIENTS B,C,D IN A CUBIC SPLINE INTERPOLATION. -C INTERPOLATION SUBROUTINES ARE TAKEN FROM -C G.E. FORSYTHE, M.A. MALCOLM AND C.B. MOLER, -C COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS (PRENTICE-HALL, 1977). - IMPLICIT double precision (A-H,O-Z) - DIMENSION X(N),Y(N),B(N),C(N),D(N) - NM1=N-1 - IF(N.LT.2) RETURN - IF(N.LT.3) GO TO 250 - D(1)=X(2)-X(1) - C(2)=(Y(2)-Y(1))/D(1) - DO 210 I=2,NM1 - D(I)=X(I+1)-X(I) - B(I)=2.0D0*(D(I-1)+D(I)) - C(I+1)=(Y(I+1)-Y(I))/D(I) - C(I)=C(I+1)-C(I) - 210 CONTINUE - B(1)=-D(1) - B(N)=-D(N-1) - C(1)=0.0D0 - C(N)=0.0D0 - IF(N.EQ.3) GO TO 215 - C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1)) - C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3)) - C(1)=C(1)*D(1)**2.0D0/(X(4)-X(1)) - C(N)=-C(N)*D(N-1)**2.0D0/(X(N)-X(N-3)) - 215 CONTINUE - DO 220 I=2,N - T=D(I-1)/B(I-1) - B(I)=B(I)-T*D(I-1) - C(I)=C(I)-T*C(I-1) - 220 CONTINUE - C(N)=C(N)/B(N) - DO 230 IB=1,NM1 - I=N-IB - C(I)=(C(I)-D(I)*C(I+1))/B(I) - 230 CONTINUE - B(N)=(Y(N)-Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2.0D0*C(N)) - DO 240 I=1,NM1 - B(I)=(Y(I+1)-Y(I))/D(I)-D(I)*(C(I+1)+2.0D0*C(I)) - D(I)=(C(I+1)-C(I))/D(I) - C(I)=3.0D0*C(I) - 240 CONTINUE - C(N)=3.0D0*C(N) - D(N)=D(N-1) - RETURN - 250 CONTINUE - B(1)=(Y(2)-Y(1))/(X(2)-X(1)) - C(1)=0.0D0 - D(1)=0.0D0 - B(2)=B(1) - C(2)=0.0D0 - D(2)=0.0D0 - RETURN - END -c -c*************************************************************************** -C --------------------------------------------------------------------- - double precision FUNCTION SEVAL(N,XX,X,Y,B,C,D) -C --------------------------------------------------------------------- -c*************************************************************************** -C CALCULATE THE DISTRIBUTION AT XX BY CUBIC SPLINE INTERPOLATION. - implicit double precision(A-H,O-Z) - DIMENSION X(N),Y(N),B(N),C(N),D(N) - DATA I/1/ - IF(I.GE.N) I=1 - IF(XX.LT.X(I)) GO TO 310 - IF(XX.LE.X(I+1)) GO TO 330 - 310 CONTINUE - I=1 - J=N+1 - 320 CONTINUE - K=(I+J)/2 - IF(XX.LT.X(K)) J=K - IF(XX.GE.X(K)) I=K - IF(J.GT.I+1) GO TO 320 - 330 CONTINUE - DX=XX-X(I) - SEVAL=Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I))) - RETURN - END -c -c*************************************************************************** diff --git a/GeneratorInterface/HydjetInterface/test/HydjetAnalyzer.cc b/GeneratorInterface/HydjetInterface/test/HydjetAnalyzer.cc index ba6ff09d1be02..e59d53065a13b 100644 --- a/GeneratorInterface/HydjetInterface/test/HydjetAnalyzer.cc +++ b/GeneratorInterface/HydjetInterface/test/HydjetAnalyzer.cc @@ -1,20 +1,9 @@ -// -*- C++ -*- -// -// Package: HydjetAnalyzer -// Class: HydjetAnalyzer -// -/**\class HydjetAnalyzer HydjetAnalyzer.cc yetkin/HydjetAnalyzer/src/HydjetAnalyzer.cc - - Description: - - Implementation: - +/** + \class HydjetAnalyzer + \brief HepMC events analyzer + \version 2.0 + \authors Yetkin Yilmaz, Andrey Belyaev */ -// -// Original Author: Yetkin Yilmaz -// Created: Tue Dec 18 09:44:41 EST 2007 -// -// // system include files #include @@ -54,11 +43,14 @@ #include "TFile.h" #include "TNtuple.h" -using namespace std; +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +using namespace edm; +using namespace std; static const int MAXPARTICLES = 5000000; static const int ETABINS = 3; // Fix also in branch string - // // class decleration // @@ -71,11 +63,14 @@ struct HydjetEvent { float nhard; float phi0; float scale; - int n[ETABINS]; float ptav[ETABINS]; - int mult; + float px[MAXPARTICLES]; + float py[MAXPARTICLES]; + float pz[MAXPARTICLES]; + float e[MAXPARTICLES]; + float pseudoRapidity[MAXPARTICLES]; float pt[MAXPARTICLES]; float eta[MAXPARTICLES]; float phi[MAXPARTICLES]; @@ -87,36 +82,27 @@ struct HydjetEvent { float vz; float vr; }; - class HydjetAnalyzer : public edm::EDAnalyzer { public: - explicit HydjetAnalyzer(const edm::ParameterSet&); + explicit HydjetAnalyzer(const edm::ParameterSet &); ~HydjetAnalyzer(); private: - virtual void beginRun(const edm::Run&, const edm::EventSetup&); + virtual void beginRun(const edm::Run &, const edm::EventSetup &); virtual void beginJob(); - virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void analyze(const edm::Event &, const edm::EventSetup &); virtual void endJob(); - // ----------member data --------------------------- - std::ofstream out_b; std::string fBFileName; - std::ofstream out_n; std::string fNFileName; - std::ofstream out_m; std::string fMFileName; - - TTree* hydjetTree_; + TTree *hydjetTree_; HydjetEvent hev_; - - TNtuple* nt; - + TNtuple *nt; std::string output; // Output filename - bool doAnalysis_; bool printLists_; bool doCF_; @@ -124,18 +110,169 @@ class HydjetAnalyzer : public edm::EDAnalyzer { bool useHepMCProduct_; bool doHI_; bool doParticles_; - double etaMax_; double ptMin_; - edm::InputTag src_; + bool doHistos_, userHistos_; + + float *ptBins; + float *etaBins; + float *phiBins; + float *v2ptBins; + float *v2etaBins; + + vector uPtBins_; + vector uEtaBins_; + vector uPhiBins_; + vector uV2ptBins_; + vector uV2etaBins_; + + int uPDG_1; + int uPDG_2; + int uPDG_3; + int uStatus_; + int nintPt = 0; + int nintEta = 0; + int nintPhi = 0; + int nintV2pt = 0; + int nintV2eta = 0; + + double minPt = 0.; + double minEta = 0.; + double minPhi = 0.; + double minV2pt = 0.; + double minV2eta = 0.; + + double maxPt = 0.; + double maxEta = 0.; + double maxPhi = 0.; + double maxV2pt = 0.; + double maxV2eta = 0.; + + double upTetaCut_ = 0.; + double downTetaCut_ = -1.; + const double pi = 3.14159265358979; + + edm::EDGetTokenT srcT_; + edm::EDGetTokenT > srcTmix_; + + // edm::InputTag src_; edm::InputTag genParticleSrc_; edm::InputTag genHIsrc_; edm::InputTag simVerticesTag_; - edm::ESHandle pdt; edm::Service f; -}; + //common + TH1D *dhpt; + TH1D *dhpt_ch; + + TH1D *dhphi; + TH1D *dhphi_ch; + + TH1D *dhv2pt; + TH1D *dhv0pt; + + TH1D *dhv2eta; + TH1D *dhv0eta; + + TH1D *dheta; + TH1D *dheta_ch; + + TH1D *hNev; + + //Users + TH1D *uhpt; + TH1D *uhpth; + TH1D *uhptj; + + TH1D *uhpt_db; + TH1D *uhpth_db; + TH1D *uhptj_db; + + TH1D *uhNpart; + TH1D *uhNparth; + TH1D *uhNpartj; + + TH1D *uhNpart_db; + TH1D *uhNparth_db; + TH1D *uhNpartj_db; + + TH1D *uhPtNpart; + TH1D *uhPtNparth; + TH1D *uhPtNpartj; + + TH1D *uhPtNpart_db; + TH1D *uhPtNparth_db; + TH1D *uhPtNpartj_db; + + TH1D *uhv2Npart; + TH1D *uhv2Nparth; + TH1D *uhv2Npartj; + + TH1D *uhv2Npart_db; + TH1D *uhv2Nparth_db; + TH1D *uhv2Npartj_db; + + TH1D *uheta; + TH1D *uhetah; + TH1D *uhetaj; + + TH1D *uhphi; + TH1D *uhphih; + TH1D *uhphij; + + TH1D *uhv2pt; + TH1D *uhv2pth; + TH1D *uhv2ptj; + + TH1D *uhv3pt; + TH1D *uhv4pt; + TH1D *uhv5pt; + TH1D *uhv6pt; + + TH1D *uhv0pt; + TH1D *uhv0pth; + TH1D *uhv0ptj; + + TH1D *uhv2pt_db; + TH1D *uhv2pth_db; + TH1D *uhv2ptj_db; + + TH1D *uhv3pt_db; + TH1D *uhv4pt_db; + TH1D *uhv5pt_db; + TH1D *uhv6pt_db; + + TH1D *uhv0pt_db; + TH1D *uhv0pth_db; + TH1D *uhv0ptj_db; + + TH1D *uhv2eta; + TH1D *uhv2etah; + TH1D *uhv2etaj; + + TH1D *uhv3eta; + TH1D *uhv4eta; + TH1D *uhv5eta; + TH1D *uhv6eta; + + TH1D *uhv0eta; + TH1D *uhv0etah; + TH1D *uhv0etaj; + + TH1D *uhv2eta_db; + TH1D *uhv2etah_db; + TH1D *uhv2etaj_db; + + TH1D *uhv3eta_db; + TH1D *uhv4eta_db; + TH1D *uhv5eta_db; + TH1D *uhv6eta_db; + + TH1D *uhv0eta_db; + TH1D *uhv0etah_db; + TH1D *uhv0etaj_db; +}; // // constants, enums and typedefs // @@ -147,12 +284,12 @@ class HydjetAnalyzer : public edm::EDAnalyzer { // // constructors and destructor // -HydjetAnalyzer::HydjetAnalyzer(const edm::ParameterSet& iConfig) { +HydjetAnalyzer::HydjetAnalyzer(const edm::ParameterSet &iConfig) { //now do what ever initialization is needed fBFileName = iConfig.getUntrackedParameter("output_b", "b_values.txt"); fNFileName = iConfig.getUntrackedParameter("output_n", "n_values.txt"); fMFileName = iConfig.getUntrackedParameter("output_m", "m_values.txt"); - doAnalysis_ = iConfig.getUntrackedParameter("doAnalysis", true); + doAnalysis_ = iConfig.getUntrackedParameter("doAnalysis", false); useHepMCProduct_ = iConfig.getUntrackedParameter("useHepMCProduct", true); printLists_ = iConfig.getUntrackedParameter("printLists", false); doCF_ = iConfig.getUntrackedParameter("doMixed", false); @@ -160,14 +297,128 @@ HydjetAnalyzer::HydjetAnalyzer(const edm::ParameterSet& iConfig) { if (doVertex_) { simVerticesTag_ = iConfig.getParameter("simVerticesTag"); } - etaMax_ = iConfig.getUntrackedParameter("etaMax", 2); + etaMax_ = iConfig.getUntrackedParameter("etaMax", 2.); ptMin_ = iConfig.getUntrackedParameter("ptMin", 0); - src_ = iConfig.getUntrackedParameter("src", edm::InputTag("VtxSmeared")); + srcT_ = mayConsume( + iConfig.getUntrackedParameter("src", edm::InputTag("generator", "unsmeared"))); + srcTmix_ = consumes >( + iConfig.getUntrackedParameter("srcMix", edm::InputTag("mix", "generatorSmeared"))); + + // src_ = iConfig.getUntrackedParameter("src", edm::InputTag("VtxSmeared")); genParticleSrc_ = iConfig.getUntrackedParameter("src", edm::InputTag("hiGenParticles")); genHIsrc_ = iConfig.getUntrackedParameter("src", edm::InputTag("heavyIon")); - doParticles_ = iConfig.getUntrackedParameter("doParticles", true); -} + doParticles_ = iConfig.getUntrackedParameter("doParticles", false); + + doHistos_ = iConfig.getUntrackedParameter("doHistos", false); + + if (doHistos_) { + userHistos_ = iConfig.getUntrackedParameter("userHistos", false); + if (userHistos_) { + uStatus_ = iConfig.getUntrackedParameter("uStatus"); + uPDG_1 = iConfig.getUntrackedParameter("uPDG_1"); + uPDG_2 = iConfig.getUntrackedParameter("uPDG_2", uPDG_1); + uPDG_3 = iConfig.getUntrackedParameter("uPDG_3", uPDG_1); + upTetaCut_ = iConfig.getUntrackedParameter("uPTetaCut", 0.8); + downTetaCut_ = iConfig.getUntrackedParameter("dPTetaCut", -1.); + uPtBins_ = iConfig.getUntrackedParameter >("PtBins"); + uEtaBins_ = iConfig.getUntrackedParameter >("EtaBins"); + uPhiBins_ = iConfig.getUntrackedParameter >("PhiBins"); + uV2ptBins_ = iConfig.getUntrackedParameter >("v2PtBins"); + uV2etaBins_ = iConfig.getUntrackedParameter >("v2EtaBins"); + + //Pt + int PtSize = uPtBins_.size(); + if (PtSize > 1) { + ptBins = new float[PtSize]; + nintPt = PtSize - 1; + for (int k = 0; k < PtSize; k++) { + ptBins[k] = uPtBins_[k]; + } + } else { + nintPt = iConfig.getUntrackedParameter("nintPt"); + maxPt = iConfig.getUntrackedParameter("maxPt"); + minPt = iConfig.getUntrackedParameter("minPt"); + ptBins = new float[nintPt + 1]; + for (int k = 0; k < nintPt + 1; k++) { + ptBins[k] = minPt + k * ((maxPt - minPt) / nintPt); + } + } + //Eta + int EtaSize = uEtaBins_.size(); + if (EtaSize > 1) { + etaBins = new float[EtaSize]; + nintEta = EtaSize - 1; + for (int k = 0; k < EtaSize; k++) { + etaBins[k] = uEtaBins_[k]; + } + } else { + nintEta = iConfig.getUntrackedParameter("nintEta"); + maxEta = iConfig.getUntrackedParameter("maxEta"); + minEta = iConfig.getUntrackedParameter("minEta"); + etaBins = new float[nintEta + 1]; + for (int k = 0; k < nintEta + 1; k++) { + etaBins[k] = minEta + k * ((maxEta - minEta) / nintEta); + } + } + + //Phi + int PhiSize = uPhiBins_.size(); + if (PhiSize > 1) { + phiBins = new float[PhiSize]; + nintPhi = PhiSize - 1; + for (int k = 0; k < PhiSize; k++) { + phiBins[k] = uPhiBins_[k]; + } + } else { + nintPhi = iConfig.getUntrackedParameter("nintPhi"); + maxPhi = iConfig.getUntrackedParameter("maxPhi"); + minPhi = iConfig.getUntrackedParameter("minPhi"); + phiBins = new float[nintPhi + 1]; + for (int k = 0; k < nintPhi + 1; k++) { + phiBins[k] = minPhi + k * ((maxPhi - minPhi) / nintPhi); + } + } + + //v2Pt + int v2PtSize = uV2ptBins_.size(); + if (v2PtSize > 1) { + v2ptBins = new float[v2PtSize]; + nintV2pt = v2PtSize - 1; + for (int k = 0; k < v2PtSize; k++) { + v2ptBins[k] = uV2ptBins_[k]; + } + } else { + nintV2pt = iConfig.getUntrackedParameter("nintV2pt"); + maxV2pt = iConfig.getUntrackedParameter("maxV2pt"); + minV2pt = iConfig.getUntrackedParameter("minV2pt"); + v2ptBins = new float[nintV2pt + 1]; + for (int k = 0; k < nintV2pt + 1; k++) { + v2ptBins[k] = minV2pt + k * ((maxV2pt - minV2pt) / nintV2pt); + } + } + + //v2Eta + int v2EtaSize = uV2etaBins_.size(); + if (v2EtaSize > 1) { + v2etaBins = new float[v2EtaSize]; + nintV2eta = v2EtaSize - 1; + for (int k = 0; k < v2EtaSize; k++) { + v2etaBins[k] = uV2etaBins_[k]; + } + } else { + nintV2eta = iConfig.getUntrackedParameter("nintV2eta"); + maxV2eta = iConfig.getUntrackedParameter("maxV2eta"); + minV2eta = iConfig.getUntrackedParameter("minV2eta"); + v2etaBins = new float[nintV2eta + 1]; + for (int k = 0; k < nintV2eta + 1; k++) { + v2etaBins[k] = minV2eta + k * ((maxV2eta - minV2eta) / nintV2eta); + } + } + + } //user histo + } //do histo +} HydjetAnalyzer::~HydjetAnalyzer() { // do anything here that needs to be done at desctruction time // (e.g. close files, deallocate resources etc.) @@ -176,50 +427,47 @@ HydjetAnalyzer::~HydjetAnalyzer() { // // member functions // - // ------------ method called to for each event ------------ -void HydjetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) { +void HydjetAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) { using namespace edm; using namespace HepMC; - hev_.event = iEvent.id().event(); for (int ieta = 0; ieta < ETABINS; ++ieta) { hev_.n[ieta] = 0; hev_.ptav[ieta] = 0; } hev_.mult = 0; - - double phi0 = 0; - double b = -1; - double scale = -1; + double phi0 = 0.; + double phi3 = 0.; + double b = -1.; + double v2, v3, v4, v5, v6; + double scale = -1.; int npart = -1; int ncoll = -1; int nhard = -1; - double vx = -99; - double vy = -99; - double vz = -99; - double vr = -99; - const GenEvent* evt; - + double vx = -99.; + double vy = -99.; + double vz = -99.; + double vr = -99.; + const GenEvent *evt; int nmix = -1; int np = 0; int sig = -1; int src = -1; - if (useHepMCProduct_) { + //_______________________________________________________________________________________ if (doCF_) { Handle > cf; - iEvent.getByLabel(InputTag("mix", "source"), cf); + iEvent.getByToken(srcTmix_, cf); MixCollection mix(cf.product()); nmix = mix.size(); - cout << "Mix Collection Size: " << mix << endl; - + //cout << "Mix Collection Size: " << mix.size() <<", pileup size: " <::iterator mbegin = mix.begin(); MixCollection::iterator mend = mix.end(); - for (MixCollection::iterator mixit = mbegin; mixit != mend; ++mixit) { - const GenEvent* subevt = (*mixit).GetEvent(); + const GenEvent *subevt = (*mixit).GetEvent(); int all = subevt->particles_size(); + //cout << "Subevent size: " << all << " Subevent type (1-signal): "<< (mixit).getTrigger()<<" Source type (pileup=0, cosmics=1, beam halo+ =2, beam halo- =3): "<< (mixit).getSourceType()<<" Bunchcrossing number: "<< mixit.bunch()<<" Impact: " <heavy_ion()->impact_parameter()<particles_begin(); HepMC::GenEvent::particle_const_iterator end = subevt->particles_end(); @@ -229,15 +477,15 @@ void HydjetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iS float eta = (*it)->momentum().eta(); float phi = (*it)->momentum().phi(); float pt = (*it)->momentum().perp(); - const ParticleData* part = pdt->particle(pdg_id); + const ParticleData *part = pdt->particle(pdg_id); int charge = static_cast(part->charge()); - hev_.pt[hev_.mult] = pt; hev_.eta[hev_.mult] = eta; hev_.phi[hev_.mult] = phi; hev_.pdg[hev_.mult] = pdg_id; hev_.chg[hev_.mult] = charge; + //cout << "Mix Particles: pt= " << pt<<" eta="<particle(pdg_id); + } +*/ + // int charge = static_cast(part->charge()); + hev_.px[hev_.mult] = px; + hev_.py[hev_.mult] = py; + hev_.pz[hev_.mult] = pz; + hev_.e[hev_.mult] = e; + hev_.pseudoRapidity[hev_.mult] = pseudoRapidity; hev_.pt[hev_.mult] = pt; hev_.eta[hev_.mult] = eta; hev_.phi[hev_.mult] = phi; hev_.pdg[hev_.mult] = pdg_id; hev_.chg[hev_.mult] = charge; + phi = phi - phi0; + /// + double phiTrue; + if (phi > pi) { + phiTrue = phi - (2 * pi); + } else if (phi < (-1 * pi)) { + phiTrue = phi + (2 * pi); + } else { + phiTrue = phi; + } + /// + v2 = TMath::Cos(2 * (phiTrue)); + v3 = TMath::Cos(3 * (phiTrue - phi3)); + v4 = TMath::Cos(4 * (phiTrue)); + v5 = TMath::Cos(5 * (phiTrue - phi3)); + v6 = TMath::Cos(6 * (phiTrue)); + + if (doHistos_) { + //common histos + if ((*it)->status() < 10) { //status 1 + dheta->Fill(eta); + + if (TMath::Abs(eta) < 0.8) { + dhpt->Fill(pt); + dhphi->Fill(phiTrue); + } + + if (TMath::Abs(pdg_id) == 211 || TMath::Abs(pdg_id) == 321 || TMath::Abs(pdg_id) == 2212) { //ch + + if (TMath::Abs(eta) < 0.8) { + dhv0pt->Fill(pt, 1.); + dhv2pt->Fill(pt, v2); + dhpt_ch->Fill(pt); + dhphi_ch->Fill(phiTrue); + } + + dhv0eta->Fill(eta, 1.); + dhv2eta->Fill(eta, v2); + dheta_ch->Fill(eta); + } //ch + } //status 1 + + //user histos + if (userHistos_ && ((uStatus_ == 3) || (((*it)->status() < 10) && (uStatus_ == 1)) || + (((*it)->status() > 10) && (uStatus_ == 2)))) { //user status + + //set1 + if (TMath::Abs(pdg_id) == uPDG_1 || TMath::Abs(pdg_id) == uPDG_2 || + TMath::Abs(pdg_id) == uPDG_3) { //uPDG + if ((uStatus_ == 3) && ((*it)->status() < 10)) + cout << "ustatus=3, but stab. part. found!!!" << endl; + + if (TMath::Abs(eta) > downTetaCut_ && TMath::Abs(eta) < upTetaCut_) { //eta cut + + uhv0pt->Fill(pt, 1.); + uhv2pt->Fill(pt, v2); + uhv3pt->Fill(pt, v3); + uhv4pt->Fill(pt, v4); + uhv5pt->Fill(pt, v5); + uhv6pt->Fill(pt, v6); + + uhv0pt_db->Fill(pt, 1.); + uhv2pt_db->Fill(pt, v2); + uhv3pt_db->Fill(pt, v3); + uhv4pt_db->Fill(pt, v4); + uhv5pt_db->Fill(pt, v5); + uhv6pt_db->Fill(pt, v6); + + if (pt >= 1.5 && pt < 10.) { + uhv2Npart->Fill(npart, v2); + uhv2Npart_db->Fill(npart, v2); + + uhPtNpart->Fill(npart, pt); + uhPtNpart_db->Fill(npart, pt); + + uhNpart->Fill(npart, 1.); + uhNpart_db->Fill(npart, 1.); + } + + uhpt->Fill(pt); + uhpt_db->Fill(pt); + uhphi->Fill(phiTrue); + + if (((*it)->status() == 16) || ((*it)->status() == 6)) { //hydro + uhv0pth->Fill(pt, 1.); + uhv2pth->Fill(pt, v2); + + uhv0pth_db->Fill(pt, 1.); + uhv2pth_db->Fill(pt, v2); + + if (pt >= 1.5 && pt < 10.) { + uhv2Nparth->Fill(npart, v2); + uhv2Nparth_db->Fill(npart, v2); + } + + uhPtNparth->Fill(npart, pt); + uhPtNparth_db->Fill(npart, pt); + + uhpth->Fill(pt); + uhpth_db->Fill(pt); + uhphih->Fill(phiTrue); + } + + if (((*it)->status() == 17) || ((*it)->status() == 7)) { //jet + uhv0ptj->Fill(pt, 1.); + uhv2ptj->Fill(pt, v2); + + uhv0ptj_db->Fill(pt, 1.); + uhv2ptj_db->Fill(pt, v2); + + if (pt >= 1.5 && pt < 10.) { + uhv2Npartj->Fill(npart, v2); + uhv2Npartj_db->Fill(npart, v2); + } + + uhPtNpartj->Fill(npart, pt); + uhPtNpartj_db->Fill(npart, pt); + + uhptj->Fill(pt); + uhptj_db->Fill(pt); + uhphij->Fill(phiTrue); + } + } //eta cut + + uheta->Fill(eta); + + uhv0eta->Fill(eta, 1.); + uhv2eta->Fill(eta, v2); + uhv3eta->Fill(eta, v3); + uhv4eta->Fill(eta, v4); + uhv5eta->Fill(eta, v5); + uhv6eta->Fill(eta, v6); + + uhv0eta_db->Fill(eta, 1.); + uhv2eta_db->Fill(eta, v2); + uhv3eta_db->Fill(eta, v3); + uhv4eta_db->Fill(eta, v4); + uhv5eta_db->Fill(eta, v5); + uhv6eta_db->Fill(eta, v6); + + if (((*it)->status() == 16) || ((*it)->status() == 6)) { //hydro + uhv2etah->Fill(eta, v2); + uhv0etah->Fill(eta, 1.); + + uhv2etah_db->Fill(eta, v2); + uhv0etah_db->Fill(eta, 1.); + + uhetah->Fill(eta); + } + if (((*it)->status() == 17) || ((*it)->status() == 7)) { //jet + uhv2etaj->Fill(eta, v2); + uhv0etaj->Fill(eta, 1.); + + uhv2etaj_db->Fill(eta, v2); + uhv0etaj_db->Fill(eta, 1.); + + uhetaj->Fill(eta); + } + } //uPDG + + } //user status + + } //doHistos_ + eta = fabs(eta); int etabin = 0; if (eta > 0.5) @@ -305,11 +744,11 @@ void HydjetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iS } } } - } else { + } else { // not HepMC edm::Handle parts; iEvent.getByLabel(genParticleSrc_, parts); for (unsigned int i = 0; i < parts->size(); ++i) { - const reco::GenParticle& p = (*parts)[i]; + const reco::GenParticle &p = (*parts)[i]; hev_.pt[hev_.mult] = p.pt(); hev_.eta[hev_.mult] = p.eta(); hev_.phi[hev_.mult] = p.phi(); @@ -366,14 +805,18 @@ void HydjetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iS hev_.vz = vz; hev_.vr = vr; - nt->Fill(nmix, np, src, sig); + if (doAnalysis_) { + nt->Fill(nmix, np, src, sig); + hydjetTree_->Fill(); + } - hydjetTree_->Fill(); + //event counter + if (doHistos_) { + hNev->Fill(1., 1); + } } - // ------------ method called once each job just before starting event loop ------------ -void HydjetAnalyzer::beginRun(const edm::Run&, const edm::EventSetup& iSetup) { iSetup.getData(pdt); } - +void HydjetAnalyzer::beginRun(const edm::Run &, const edm::EventSetup &iSetup) { iSetup.getData(pdt); } void HydjetAnalyzer::beginJob() { if (printLists_) { out_b.open(fBFileName.c_str()); @@ -387,9 +830,147 @@ void HydjetAnalyzer::beginJob() { throw cms::Exception("BadFile") << "Can\'t open file " << fMFileName; } + if (doHistos_) { + if (userHistos_) { + cout << "---------------------------------------------------------------INPUT------------------------------------" + "-------------------------------" + << endl; + cout << "etaCut for pT = " << downTetaCut_ << " - " << upTetaCut_ << endl; + + //pt + uhpt = new TH1D("uhpt", "uhpt", nintPt, ptBins); + uhptj = new TH1D("uhptj", "uhptj", nintPt, ptBins); + uhpth = new TH1D("uhpth", "uhpth", nintPt, ptBins); + + //pt_db + uhpt_db = new TH1D("uhpt_db", "uhpt_db", 1000, 0.0000000000001, 100.); + uhptj_db = new TH1D("uhptj_db", "uhptj_db", 1000, 0.0000000000001, 100.); + uhpth_db = new TH1D("uhpth_db", "uhpth_db", 1000, 0.0000000000001, 100.); + + //eta + uheta = new TH1D("uheta", "uheta", nintEta, etaBins); + uhetaj = new TH1D("uhetaj", "uhetaj", nintEta, etaBins); + uhetah = new TH1D("uhetah", "uhetah", nintEta, etaBins); + + //phi + uhphi = new TH1D("uhphi", "uhphi", nintPhi, phiBins); + uhphij = new TH1D("uhphij", "uhphij", nintPhi, phiBins); + uhphih = new TH1D("uhphih", "uhphih", nintPhi, phiBins); + + const int NbinNpar = 5; + const double BinsNpart[NbinNpar + 1] = {0., 29., 90., 202., 346., 400.}; + + //ptNpart + uhNpart = new TH1D("uhNpart", "uhNpart", NbinNpar, BinsNpart); + uhNpartj = new TH1D("uhNpartj", "uhNpartj", NbinNpar, BinsNpart); + uhNparth = new TH1D("uhNparth", "uhNparth", NbinNpar, BinsNpart); + + //ptNpart_db + uhNpart_db = new TH1D("uhNpart_db", "uhNpart_db", 400, 0., 400.); + uhNpartj_db = new TH1D("uhNpartj_db", "uhNpartj_db", 400, 0., 400.); + uhNparth_db = new TH1D("uhNparth_db", "uhNparth_db", 400, 0., 400.); + + //ptNpart + uhPtNpart = new TH1D("uhptNpart", "uhptNpart", NbinNpar, BinsNpart); + uhPtNpartj = new TH1D("uhptNpartj", "uhptNpartj", NbinNpar, BinsNpart); + uhPtNparth = new TH1D("uhptNparth", "uhptNparth", NbinNpar, BinsNpart); + + //ptNpart_db + uhPtNpart_db = new TH1D("uhptNpart_db", "uhptNpart_db", 400, 0., 400.); + uhPtNpartj_db = new TH1D("uhptNpartj_db", "uhptNpartj_db", 400, 0., 400.); + uhPtNparth_db = new TH1D("uhptNparth_db", "uhptNparth_db", 400, 0., 400.); + + //v2Npart + uhv2Npart = new TH1D("uhv2Npart", "uhv2Npart", NbinNpar, BinsNpart); + uhv2Npartj = new TH1D("uhv2Npartj", "uhv2Npartj", NbinNpar, BinsNpart); + uhv2Nparth = new TH1D("uhv2Nparth", "uhv2Nparth", NbinNpar, BinsNpart); + + //v2Npart_db + uhv2Npart_db = new TH1D("uhv2Npart_db", "uhv2Npart_db", 400, 0., 400.); + uhv2Npartj_db = new TH1D("uhv2Npartj_db", "uhv2Npartj_db", 400, 0., 400.); + uhv2Nparth_db = new TH1D("uhv2Nparth_db", "uhv2Nparth_db", 400, 0., 400.); + + //v0pt + uhv0pt = new TH1D("uhv0pt", "uhv0pt", nintV2pt, v2ptBins); + uhv0ptj = new TH1D("uhv0ptj", "uhv0ptj", nintV2pt, v2ptBins); + uhv0pth = new TH1D("uhv0pth", "uhv0pth", nintV2pt, v2ptBins); + + //v2pt + uhv2pt = new TH1D("uhv2pt", "uhv2pt", nintV2pt, v2ptBins); + uhv2ptj = new TH1D("uhv2ptj", "uhv2ptj", nintV2pt, v2ptBins); + uhv2pth = new TH1D("uhv2pth", "uhv2pth", nintV2pt, v2ptBins); + + uhv3pt = new TH1D("uhv3pt", "uhv3pt", nintV2pt, v2ptBins); + uhv4pt = new TH1D("uhv4pt", "uhv4pt", nintV2pt, v2ptBins); + uhv5pt = new TH1D("uhv5pt", "uhv5pt", nintV2pt, v2ptBins); + uhv6pt = new TH1D("uhv6pt", "uhv6pt", nintV2pt, v2ptBins); + + //v0pt + uhv0pt_db = new TH1D("uhv0pt_db", "uhv0pt_db", 200, 0.0, 10.); + uhv0ptj_db = new TH1D("uhv0ptj_db", "uhv0ptj_db", 200, 0.0, 10.); + uhv0pth_db = new TH1D("uhv0pth_db", "uhv0pth_db", 200, 0.0, 10.); + + //v2pt_db + uhv2pt_db = new TH1D("uhv2pt_db", "uhv2pt_db", 200, 0.0, 10.); + uhv2ptj_db = new TH1D("uhv2ptj_db", "uhv2ptj_db", 200, 0.0, 10.); + uhv2pth_db = new TH1D("uhv2pth_db", "uhv2pth_db", 200, 0.0, 10.); + + uhv3pt_db = new TH1D("uhv3pt_db", "uhv3pt_db", 200, 0.0, 10.); + uhv4pt_db = new TH1D("uhv4pt_db", "uhv4pt_db", 200, 0.0, 10.); + uhv5pt_db = new TH1D("uhv5pt_db", "uhv5pt_db", 200, 0.0, 10.); + uhv6pt_db = new TH1D("uhv6pt_db", "uhv6pt_db", 200, 0.0, 10.); + + //v0eta + uhv0eta = new TH1D("uhv0eta", "uhv0eta", nintV2eta, v2etaBins); + uhv0etaj = new TH1D("uhv0etaj", "uhv0etaj", nintV2eta, v2etaBins); + uhv0etah = new TH1D("uhv0etah", "uhv0etah", nintV2eta, v2etaBins); + + //v2eta + uhv2eta = new TH1D("uhv2eta", "uhv2eta", nintV2eta, v2etaBins); + uhv2etaj = new TH1D("uhv2etaj", "uhv2etaj", nintV2eta, v2etaBins); + uhv2etah = new TH1D("uhv2etah", "uhv2etah", nintV2eta, v2etaBins); + + uhv3eta = new TH1D("uhv3eta", "uhv3eta", nintV2eta, v2etaBins); + uhv4eta = new TH1D("uhv4eta", "uhv4eta", nintV2eta, v2etaBins); + uhv5eta = new TH1D("uhv5eta", "uhv5eta", nintV2eta, v2etaBins); + uhv6eta = new TH1D("uhv6eta", "uhv6eta", nintV2eta, v2etaBins); + + //v0eta_db + uhv0eta_db = new TH1D("uhv0eta_db", "uhv0eta_db", 200, -5, 5.); + uhv0etaj_db = new TH1D("uhv0etaj_db", "uhv0etaj_db", 200, -5, 5.); + uhv0etah_db = new TH1D("uhv0etah_db", "uhv0etah_db", 200, -5, 5.); + + //v2eta_db + uhv2eta_db = new TH1D("uhv2eta_db", "uhv2eta_db", 200, -5, 5.); + uhv2etaj_db = new TH1D("uhv2etaj_db", "uhv2etaj_db", 200, -5, 5.); + uhv2etah_db = new TH1D("uhv2etah_db", "uhv2etah_db", 200, -5, 5.); + + uhv3eta_db = new TH1D("uhv3eta_db", "uhv3eta_db", 200, -5, 5.); + uhv4eta_db = new TH1D("uhv4eta_db", "uhv4eta_db", 200, -5, 5.); + uhv5eta_db = new TH1D("uhv5eta_db", "uhv5eta_db", 200, -5, 5.); + uhv6eta_db = new TH1D("uhv6eta_db", "uhv6eta_db", 200, -5, 5.); + } + + dhpt = new TH1D("dhpt", "dhpt", 1000, 0.0000000000001, 100.); + dhpt_ch = new TH1D("dhpt_ch", "dhpt_ch", 1000, 0.0000000000001, 100.); + + dhphi = new TH1D("dhphi", "dhphi", 1000, -3.14159265358979, 3.14159265358979); + dhphi_ch = new TH1D("dhphi_ch", "dhphi_ch", 1000, -3.14159265358979, 3.14159265358979); + + dhv2pt = new TH1D("dhv2pt", "dhv2pt", 200, 0.0, 10.); + dhv0pt = new TH1D("dhv0pt", "dhv0pt", 200, 0.0, 10.); + + dhv2eta = new TH1D("dhv2eta", "dhv2eta", 200, -5, 5.); + dhv0eta = new TH1D("dhv0eta", "dhv0eta", 200, -5, 5.); + + dheta = new TH1D("dheta", "dheta", 1000, -10., 10.); + dheta_ch = new TH1D("dheta_ch", "dheta_ch", 1000, -10., 10.); + + hNev = new TH1D("hNev", "hNev", 1, 0., 2.); + } + if (doAnalysis_) { nt = f->make("nt", "Mixing Analysis", "mix:np:src:sig"); - hydjetTree_ = f->make("hi", "Tree of Hydjet Events"); hydjetTree_->Branch("event", &hev_.event, "event/I"); hydjetTree_->Branch("b", &hev_.b, "b/F"); @@ -398,12 +979,15 @@ void HydjetAnalyzer::beginJob() { hydjetTree_->Branch("nhard", &hev_.nhard, "nhard/F"); hydjetTree_->Branch("phi0", &hev_.phi0, "phi0/F"); hydjetTree_->Branch("scale", &hev_.scale, "scale/F"); - hydjetTree_->Branch("n", hev_.n, "n[3]/I"); hydjetTree_->Branch("ptav", hev_.ptav, "ptav[3]/F"); - if (doParticles_) { hydjetTree_->Branch("mult", &hev_.mult, "mult/I"); + hydjetTree_->Branch("px", hev_.px, "px[mult]/F"); + hydjetTree_->Branch("py", hev_.py, "py[mult]/F"); + hydjetTree_->Branch("pz", hev_.pz, "pz[mult]/F"); + hydjetTree_->Branch("e", hev_.e, "e[mult]/F"); + hydjetTree_->Branch("pseudoRapidity", hev_.pseudoRapidity, "pseudoRapidity[mult]/F"); hydjetTree_->Branch("pt", hev_.pt, "pt[mult]/F"); hydjetTree_->Branch("eta", hev_.eta, "eta[mult]/F"); hydjetTree_->Branch("phi", hev_.phi, "phi[mult]/F"); @@ -417,9 +1001,113 @@ void HydjetAnalyzer::beginJob() { } } } - // ------------ method called once each job just after ending the event loop ------------ -void HydjetAnalyzer::endJob() {} - +void HydjetAnalyzer::endJob() { + if (doHistos_) { + dhpt->Write(); + dheta->Write(); + dhphi->Write(); + dhv0pt->Write(); + dhv2pt->Write(); + dhv0eta->Write(); + dhv2eta->Write(); + + dhpt_ch->Write(); + dheta_ch->Write(); + hNev->Write(); + if (userHistos_) { + uhpt->Write(); + uhpth->Write(); + uhptj->Write(); + + uhpt_db->Write(); + uhpth_db->Write(); + uhptj_db->Write(); + + uhNpart->Write(); + uhNparth->Write(); + uhNpartj->Write(); + + uhNpart_db->Write(); + uhNparth_db->Write(); + uhNpartj_db->Write(); + + uhPtNpart->Write(); + uhPtNparth->Write(); + uhPtNpartj->Write(); + + uhPtNpart_db->Write(); + uhPtNparth_db->Write(); + uhPtNpartj_db->Write(); + + uhv2Npart->Write(); + uhv2Nparth->Write(); + uhv2Npartj->Write(); + + uhv2Npart_db->Write(); + uhv2Nparth_db->Write(); + uhv2Npartj_db->Write(); + + uheta->Write(); + uhetah->Write(); + uhetaj->Write(); + uhphi->Write(); + uhphih->Write(); + uhphij->Write(); + + uhv0eta->Write(); + uhv0etah->Write(); + uhv0etaj->Write(); + + uhv0eta_db->Write(); + uhv0etah_db->Write(); + uhv0etaj_db->Write(); + + uhv0pt->Write(); + uhv0pth->Write(); + uhv0ptj->Write(); + + uhv0pt_db->Write(); + uhv0pth_db->Write(); + uhv0ptj_db->Write(); + + uhv2eta->Write(); + uhv2etah->Write(); + uhv2etaj->Write(); + + uhv2eta_db->Write(); + uhv2etah_db->Write(); + uhv2etaj_db->Write(); + + uhv2pt->Write(); + uhv2pth->Write(); + uhv2ptj->Write(); + + uhv2pt_db->Write(); + uhv2pth_db->Write(); + uhv2ptj_db->Write(); + + uhv3eta->Write(); + uhv4eta->Write(); + uhv5eta->Write(); + uhv6eta->Write(); + + uhv3eta_db->Write(); + uhv4eta_db->Write(); + uhv5eta_db->Write(); + uhv6eta_db->Write(); + + uhv3pt->Write(); + uhv4pt->Write(); + uhv5pt->Write(); + uhv6pt->Write(); + + uhv3pt_db->Write(); + uhv4pt_db->Write(); + uhv5pt_db->Write(); + uhv6pt_db->Write(); + } + } +} //define this as a plug-in DEFINE_FWK_MODULE(HydjetAnalyzer); diff --git a/GeneratorInterface/HydjetInterface/test/testHydjet.py b/GeneratorInterface/HydjetInterface/test/testHydjet.py index 0dc8d7e638c8c..3150e80f35f51 100644 --- a/GeneratorInterface/HydjetInterface/test/testHydjet.py +++ b/GeneratorInterface/HydjetInterface/test/testHydjet.py @@ -6,7 +6,6 @@ process.load("Configuration.StandardSequences.Services_cff") process.load("GeneratorInterface.HydjetInterface.hydjetDefault_cfi") - process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService", generator = cms.PSet( initialSeed = cms.untracked.uint32(123456789), @@ -15,14 +14,14 @@ ) process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(100) -) + input = cms.untracked.int32(5) + ) -process.ana = cms.EDAnalyzer('Hydjet2Analyzer', +process.ana = cms.EDAnalyzer('HydjetAnalyzer', doHistos = cms.untracked.bool(True), userHistos = cms.untracked.bool(False), - + # Settings for USER histos uStatus = cms.untracked.int32(2), # 1 - it's 1,2,3,4,5 of Pythia status; 2 - 11,12,13,14,15; 3 - All @@ -62,6 +61,7 @@ ) +#process.generator.signalVtx = cms.untracked.vdouble(0.,0.,0.,0.) # Signal event vertex option, to set it by hand (instead of smearing) process.TFileService = cms.Service('TFileService', fileName = cms.string('Hydjet1_MB_5020GeV.root') diff --git a/GeneratorInterface/PyquenInterface/interface/PyquenGeneratorFilter.h b/GeneratorInterface/PyquenInterface/interface/PyquenGeneratorFilter.h index 752a84f85cfb8..45ecb1b859ea5 100644 --- a/GeneratorInterface/PyquenInterface/interface/PyquenGeneratorFilter.h +++ b/GeneratorInterface/PyquenInterface/interface/PyquenGeneratorFilter.h @@ -5,6 +5,14 @@ #include "GeneratorInterface/Core/interface/GeneratorFilter.h" #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h" +namespace edm { + template <> + GeneratorFilter::GeneratorFilter(ParameterSet const& ps) + : hadronizer_(ps, consumesCollector()) { + init(ps); + } +} // namespace edm + namespace gen { typedef edm::GeneratorFilter PyquenGeneratorFilter; } diff --git a/GeneratorInterface/PyquenInterface/interface/PyquenHadronizer.h b/GeneratorInterface/PyquenInterface/interface/PyquenHadronizer.h index 3e3e2e6e48a96..2f242bb62e908 100644 --- a/GeneratorInterface/PyquenInterface/interface/PyquenHadronizer.h +++ b/GeneratorInterface/PyquenInterface/interface/PyquenHadronizer.h @@ -1,19 +1,24 @@ #ifndef Pyquen_Hadronizer_h #define Pyquen_Hadronizer_h -/** \class PyquenHadronizer - * - * Generates PYTHIA+PYQUEN ==> HepMC events - * - * Camelia Mironov - * for the Generator Interface. March 2007 - ***************************************/ +/** + \class PyquenHadronizer + \brief Interface to the PYQUEN generator (since core v. 1.5.4), produces HepMC events + \version 2.0 + \authors Camelia Mironov, Andrey Belyaev +*/ #include "GeneratorInterface/Core/interface/BaseHadronizer.h" #include "GeneratorInterface/HiGenCommon/interface/BaseHiGenEvtSelector.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "SimDataFormats/CrossingFrame/interface/MixCollection.h" +#include "SimDataFormats/PileupSummaryInfo/interface/PileupMixingContent.h" + #include #include #include @@ -29,7 +34,7 @@ namespace gen { class PyquenHadronizer : public BaseHadronizer { public: - PyquenHadronizer(const edm::ParameterSet&); + PyquenHadronizer(const edm::ParameterSet&, edm::ConsumesCollector&&); ~PyquenHadronizer() override; bool generatePartonsAndHadronize(); @@ -60,37 +65,41 @@ namespace gen { void rotateEvtPlane(HepMC::GenEvent* evt, double angle); edm::ParameterSet pset_; - double abeamtarget_; //! beam/target atomic mass number - unsigned int angularspecselector_; //! angular emitted gluon spectrum selection - //! DEFAULT= 0 -- small angular emitted gluon spectrum - //! = 1 -- broad angular emitted gluon spectrum - //! = 2 -- collinear angular emitted gluon spectrum - double bmin_; //! min impact param (fm); valid only if cflag_!=0 - double bmax_; //! max impact param (fm); valid only if cflag_!=0 - double bfixed_; //! fixed impact param (fm); valid only if cflag_=0 - int cflag_; //! centrality flag =0 fixed impact param, <>0 minbias - double comenergy; //! collision energy - bool doquench_; //! if true perform quenching (default = true) - bool doradiativeenloss_; //! DEFAULT = true - bool docollisionalenloss_; //! DEFAULT = true - bool doIsospin_; //! Run n&p with proper ratios; if false, only p+p collisions + double abeamtarget_; ///< beam/target atomic mass number + unsigned int angularspecselector_; ///< angular emitted gluon spectrum selection + ///< DEFAULT= 0 -- small angular emitted gluon spectrum + ///< = 1 -- broad angular emitted gluon spectrum + ///< = 2 -- collinear angular emitted gluon spectrum + double bmin_; ///< min impact param (fm); valid only if cflag_!=0 + double bmax_; ///< max impact param (fm); valid only if cflag_!=0 + double bfixed_; ///< fixed impact param (fm); valid only if cflag_=0 + int cflag_; ///< centrality flag =0 fixed impact param, <>0 minbias + double comenergy; ///< collision energy + bool doquench_; ///< if true perform quenching (default = true) + bool doradiativeenloss_; ///< DEFAULT = true + bool docollisionalenloss_; ///< DEFAULT = true + bool doIsospin_; ///< Run n&p with proper ratios; if false, only p+p collisions int protonSide_; bool embedding_; double evtPlane_; - double pfrac_; //! Proton fraction in the nucleus - - unsigned int nquarkflavor_; //! number of active quark flavors in qgp - //! DEFAULT=0; allowed values: 0,1,2,3. - double qgpt0_; //! initial temperature of QGP - //! DEFAULT = 1GeV; allowed range [0.2,2.0]GeV; - double qgptau0_; //! proper time of QGP formation - //! DEFAULT = 0.1 fm/c; allowed range [0.01,10.0]fm/c; - unsigned int maxEventsToPrint_; //! Events to print if verbosity - bool pythiaHepMCVerbosity_; //! HepMC verbosity flag - unsigned int pythiaPylistVerbosity_; //! Pythia PYLIST Verbosity flag + double pfrac_; ///< Proton fraction in the nucleus + + unsigned int nquarkflavor_; ///< number of active quark flavors in qgp + ///< DEFAULT=0; allowed values: 0,1,2,3. + double qgpt0_; ///< initial temperature of QGP + ///< DEFAULT = 1GeV; allowed range [0.2,2.0]GeV; + double qgptau0_; ///< proper time of QGP formation + ///< DEFAULT = 0.1 fm/c; allowed range [0.01,10.0]fm/c; + unsigned int maxEventsToPrint_; ///< Events to print if verbosity + + HepMC::FourVector* fVertex_; ///< Event signal vertex + std::vector signalVtx_; ///< Pset double vector to set event signal vertex + + bool pythiaHepMCVerbosity_; ///< HepMC verbosity flag + unsigned int pythiaPylistVerbosity_; ///< Pythia PYLIST Verbosity flag // CLHEP::HepRandomEngine* fRandomEngine; - edm::InputTag src_; + edm::EDGetTokenT > src_; Pythia6Service* pythia6Service_; std::string filterType_; BaseHiGenEvtSelector* selector_; diff --git a/GeneratorInterface/PyquenInterface/python/pyquenDefault_cfi.py b/GeneratorInterface/PyquenInterface/python/pyquenDefault_cfi.py index 86887f94564b3..422c0931410b6 100644 --- a/GeneratorInterface/PyquenInterface/python/pyquenDefault_cfi.py +++ b/GeneratorInterface/PyquenInterface/python/pyquenDefault_cfi.py @@ -14,7 +14,7 @@ angularSpectrumSelector = cms.int32(0), ## angular emitted gluon spectrum : pythiaHepMCVerbosity = cms.untracked.bool(False), PythiaParameters = cms.PSet(pyquenPythiaDefaultBlock, - parameterSets = cms.vstring('pythiaUESettings','ppJets','pythiaPromptPhotons','kinematics'), + parameterSets = cms.vstring('pyquenMain','pythiaUESettings','ppJets','pythiaPromptPhotons','kinematics'), kinematics = cms.vstring('CKIN(3) = 50','CKIN(4) = 80') ), qgpProperTimeFormation = cms.double(0.1), ## proper time of QGP formation; allowed range [0.01,10.0]fm/c; diff --git a/GeneratorInterface/PyquenInterface/python/pyquenPythiaDefault_cff.py b/GeneratorInterface/PyquenInterface/python/pyquenPythiaDefault_cff.py index 78de12a2be837..d428fab88af4f 100644 --- a/GeneratorInterface/PyquenInterface/python/pyquenPythiaDefault_cff.py +++ b/GeneratorInterface/PyquenInterface/python/pyquenPythiaDefault_cff.py @@ -5,6 +5,11 @@ from Configuration.Generator.PythiaUESettings_cfi import * pyquenPythiaDefaultBlock = cms.PSet( pythiaUESettingsBlock, + pyquenMain = cms.vstring( + 'MSTP(122)=0', # ! no printout of Pythia initialization information hereinafter + 'MSTP(128)=2' # ! instructs pythia not to put multiple copies of resonances in the event record + ), + ppDefault = cms.vstring('MSEL=1 ! QCD hight pT processes', 'CKIN(3)=7.',# ! ptMin 'MSTP(81)=0' diff --git a/GeneratorInterface/PyquenInterface/src/PyquenHadronizer.cc b/GeneratorInterface/PyquenInterface/src/PyquenHadronizer.cc index 17ac575981f20..1a9c0745b7b98 100644 --- a/GeneratorInterface/PyquenInterface/src/PyquenHadronizer.cc +++ b/GeneratorInterface/PyquenInterface/src/PyquenHadronizer.cc @@ -1,9 +1,3 @@ -/* - * - * Generates PYQUEN HepMC events - * -*/ - #include #include @@ -35,7 +29,7 @@ HepMC::IO_HEPEVT hepevtio; const std::vector PyquenHadronizer::theSharedResources = {edm::SharedResourceNames::kPythia6, gen::FortranInstance::kFortranInstance}; -PyquenHadronizer ::PyquenHadronizer(const ParameterSet& pset) +PyquenHadronizer ::PyquenHadronizer(const ParameterSet& pset, edm::ConsumesCollector&& iC) : BaseHadronizer(pset), pset_(pset), abeamtarget_(pset.getParameter("aBeamTarget")), @@ -56,10 +50,23 @@ PyquenHadronizer ::PyquenHadronizer(const ParameterSet& pset) qgpt0_(pset.getParameter("qgpInitialTemperature")), qgptau0_(pset.getParameter("qgpProperTimeFormation")), maxEventsToPrint_(pset.getUntrackedParameter("maxEventsToPrint", 1)), + fVertex_(nullptr), pythiaHepMCVerbosity_(pset.getUntrackedParameter("pythiaHepMCVerbosity", false)), pythiaPylistVerbosity_(pset.getUntrackedParameter("pythiaPylistVerbosity", 0)), pythia6Service_(new Pythia6Service(pset)), filterType_(pset.getUntrackedParameter("filterType", "None")) { + if (pset.exists("signalVtx")) + signalVtx_ = pset.getUntrackedParameter >("signalVtx"); + + if (signalVtx_.size() == 4) { + if (!fVertex_) + fVertex_ = new HepMC::FourVector(); + LogDebug("EventSignalVertex") << "Setting event signal vertex " + << " x = " << signalVtx_.at(0) << " y = " << signalVtx_.at(1) + << " z= " << signalVtx_.at(2) << " t = " << signalVtx_.at(3) << endl; + fVertex_->set(signalVtx_.at(0), signalVtx_.at(1), signalVtx_.at(2), signalVtx_.at(3)); + } + // Verbosity Level // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13 LogDebug("PYLISTverbosity") << "Pythia PYLIST verbosity level = " << pythiaPylistVerbosity_ << endl; @@ -73,12 +80,14 @@ PyquenHadronizer ::PyquenHadronizer(const ParameterSet& pset) if (embedding_) { cflag_ = 0; - src_ = pset.getParameter("backgroundLabel"); + src_ = iC.consumes >( + pset.getUntrackedParameter("backgroundLabel", edm::InputTag("mix", "generatorSmeared"))); } selector_ = HiGenEvtSelectorFactory::get(filterType_, pset); int cm = 1, va, vb, vc; PYQVER(cm, va, vb, vc); + //HepMC::HEPEVT_Wrapper::set_max_number_entries(4000); } //_____________________________________________________________________ @@ -122,13 +131,47 @@ bool PyquenHadronizer::generatePartonsAndHadronize() { // at this part, need to overwrite filter() in // PyquenGeneratorFilter - const edm::Event& e = getEDMEvent(); - if (embedding_) { - Handle input; - e.getByLabel(src_, input); - const HepMC::GenEvent* inev = input->GetEvent(); + const edm::Event& e = getEDMEvent(); + HepMC::GenVertex* genvtx = nullptr; + const HepMC::GenEvent* inev = nullptr; + Handle > cf; + e.getByToken(src_, cf); + MixCollection mix(cf.product()); + if (mix.size() < 1) { + throw cms::Exception("MatchVtx") << "Mixing has " << mix.size() << " sub-events, should have been at least 1" + << endl; + } + const HepMCProduct& bkg = mix.getObject(0); + if (!(bkg.isVtxGenApplied())) { + throw cms::Exception("MatchVtx") << "Input background does not have smeared vertex!" << endl; + } else { + inev = bkg.GetEvent(); + } + + genvtx = inev->signal_process_vertex(); + + if (!genvtx) + throw cms::Exception("MatchVtx") << "Input background does not have signal process vertex!" << endl; + + double aX, aY, aZ, aT; + + aX = genvtx->position().x(); + aY = genvtx->position().y(); + aZ = genvtx->position().z(); + aT = genvtx->position().t(); + + if (!fVertex_) { + fVertex_ = new HepMC::FourVector(); + } + + //LogInfo("MatchVtx") + std::cout << " setting vertex " + << " aX " << aX << " aY " << aY << " aZ " << aZ << " aT " << aT << endl; + fVertex_->set(aX, aY, aZ, aT); + const HepMC::HeavyIon* hi = inev->heavy_ion(); + if (hi) { bfixed_ = hi->impact_parameter(); evtPlane_ = hi->event_plane_angle(); @@ -169,6 +212,7 @@ bool PyquenHadronizer::generatePartonsAndHadronize() { call_pyhepc(1); // event information + // hepevtio.set_trust_mothers_before_daughters(true); HepMC::GenEvent* evt = hepevtio.read_next_event(); evt->set_signal_process_id(pypars.msti[0]); // type of the process @@ -178,8 +222,15 @@ bool PyquenHadronizer::generatePartonsAndHadronize() { rotateEvtPlane(evt, evtPlane_); add_heavy_ion_rec(evt); - HepMC::HEPEVT_Wrapper::check_hepevt_consistency(); - //std::cout<<"Entries number: "< HepMCEvt(new edm::HepMCProduct(evt)); + + HepMCEvt->applyVtxGen(fVertex_); + evt = new HepMC::GenEvent((*HepMCEvt->GetEvent())); + } + + // HepMC::HEPEVT_Wrapper::check_hepevt_consistency(); event().reset(evt); diff --git a/GeneratorInterface/PyquenInterface/test/testPyquen.py b/GeneratorInterface/PyquenInterface/test/testPyquen.py index 5ed87135c34f9..41d46566f8542 100644 --- a/GeneratorInterface/PyquenInterface/test/testPyquen.py +++ b/GeneratorInterface/PyquenInterface/test/testPyquen.py @@ -8,7 +8,7 @@ process.source = cms.Source("EmptySource") -process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(100) +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(1000) ) process.SimpleMemoryCheck = cms.Service('SimpleMemoryCheck', @@ -24,7 +24,3 @@ ) process.p = cms.Path(process.generator*process.ana) - - - - diff --git a/SimGeneral/MixingModule/plugins/MixingModule.cc b/SimGeneral/MixingModule/plugins/MixingModule.cc index d655eaaa1bc00..c100a042a99fc 100644 --- a/SimGeneral/MixingModule/plugins/MixingModule.cc +++ b/SimGeneral/MixingModule/plugins/MixingModule.cc @@ -72,6 +72,11 @@ namespace edm { wrapLongTimes_ = ps_mix.getParameter("WrapLongTimes"); } + skipSignal_ = false; + if (ps_mix.exists("skipSignal")) { + skipSignal_ = ps_mix.getParameter("skipSignal"); + } + ParameterSet ps = ps_mix.getParameter("mixObjects"); std::vector names = ps.getParameterNames(); for (std::vector::iterator it = names.begin(); it != names.end(); ++it) { @@ -312,14 +317,14 @@ namespace edm { void MixingModule::checkSignal(const edm::Event& e) { if (adjusters_.empty()) { for (auto const& adjuster : adjustersObjects_) { - if (adjuster->checkSignal(e)) { + if (skipSignal_ or adjuster->checkSignal(e)) { adjusters_.push_back(adjuster); } } } if (workers_.empty()) { for (auto const& worker : workersObjects_) { - if (worker->checkSignal(e)) { + if (skipSignal_ or worker->checkSignal(e)) { workers_.push_back(worker); } } @@ -351,6 +356,10 @@ namespace edm { } void MixingModule::addSignals(const edm::Event& e, const edm::EventSetup& setup) { + if (skipSignal_) { + return; + } + LogDebug("MixingModule") << "===============> adding signals for " << e.id(); accumulateEvent(e, setup); diff --git a/SimGeneral/MixingModule/plugins/MixingModule.h b/SimGeneral/MixingModule/plugins/MixingModule.h index 3825202fc75b8..fe0121b173385 100644 --- a/SimGeneral/MixingModule/plugins/MixingModule.h +++ b/SimGeneral/MixingModule/plugins/MixingModule.h @@ -107,6 +107,7 @@ namespace edm { std::vector wantedBranches_; bool useCurrentProcessOnly_; bool wrapLongTimes_; + bool skipSignal_; // Digi-producing algorithms Accumulators digiAccumulators_; diff --git a/SimGeneral/MixingModule/python/HiMixEmbGEN_cff.py b/SimGeneral/MixingModule/python/HiMixEmbGEN_cff.py new file mode 100644 index 0000000000000..2f37361906e49 --- /dev/null +++ b/SimGeneral/MixingModule/python/HiMixEmbGEN_cff.py @@ -0,0 +1,104 @@ +import FWCore.ParameterSet.Config as cms + +# configuration to model pileup for initial physics phase +#from SimGeneral.MixingModule.mixObjects_cfi import theMixObjects#, run2_GEM_2017, premix_stage1 +from SimGeneral.MixingModule.mixPoolSource_cfi import * +#from SimGeneral.MixingModule.digitizers_cfi import theDigitizers + +FileNames = cms.untracked.vstring(['/store/relval/CMSSW_7_2_0_pre7/RelValQCD_Pt_80_120_13/GEN-SIM/PRE_LS172_V11-v1/00000/16547ECB-9C4B-E411-A815-0025905964BC.root', '/store/relval/CMSSW_7_2_0_pre7/RelValQCD_Pt_80_120_13/GEN-SIM/PRE_LS172_V11-v1/00000/86C3C326-9F4B-E411-903D-0025905A48EC.root', '/store/relval/CMSSW_7_2_0_pre7/RelValQCD_Pt_80_120_13/GEN-SIM/PRE_LS172_V11-v1/00000/C48D8223-9F4B-E411-BC37-0026189438DC.root', '/store/relval/CMSSW_7_2_0_pre7/RelValQCD_Pt_80_120_13/GEN-SIM/PRE_LS172_V11-v1/00000/D070AB62-9D4B-E411-9766-002618FDA207.root']) + +mix = cms.EDProducer("MixingModule", + skipSignal = cms.bool(True), + + digitizers = cms.PSet(),#theDigitizers), + LabelPlayback = cms.string(''), + maxBunch = cms.int32(0), + minBunch = cms.int32(0), ## in terms of 25 nsec + bunchspace = cms.int32(1), ##ns + mixProdStep1 = cms.bool(False), + mixProdStep2 = cms.bool(False), + + playback = cms.untracked.bool(False), + useCurrentProcessOnly = cms.bool(False), + + input = cms.SecSource("EmbeddedRootSource", + nbPileupEvents = cms.PSet( + averageNumber = cms.double(1.0) + ), + type = cms.string('fixed'), + sequential = cms.untracked.bool(False), + fileNames = FileNames + ), + + mixObjects = cms.PSet( + +# theMixObjects + + mixHepMC = cms.PSet( + input = cms.VInputTag( + cms.InputTag("generatorSmeared","",cms.InputTag.skipCurrentProcess()), + cms.InputTag("generator","unsmeared",cms.InputTag.skipCurrentProcess()), + cms.InputTag("generator","",cms.InputTag.skipCurrentProcess()) + ), + + makeCrossingFrame = cms.untracked.bool(True), + type = cms.string('HepMCProduct') + ) + ) +) + +''' +#mix.digitizers.castor.hitsProducer = cms.InputTag("g4SimHits","CastorFI",cms.InputTag.skipCurrentProcess()) +#mix.digitizers.puVtx.vtxTag = cms.InputTag("generatorSmeared","",cms.InputTag.skipCurrentProcess()) +#mix.digitizers.puVtx.vtxFallbackTag = cms.InputTag("generator","",cms.InputTag.skipCurrentProcess()) + +mix.mixObjects.mixCH.input = cms.VInputTag( + #cms.InputTag("g4SimHits","CaloHitsTk"), cms.InputTag("g4SimHits","CastorBU"), cms.InputTag("g4SimHits","CastorPL"), cms.InputTag("g4SimHits","CastorTU"), + cms.InputTag("g4SimHits","CastorFI",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","EcalHitsEB",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","EcalHitsEE",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","EcalHitsES",cms.InputTag.skipCurrentProcess()), + #cms.InputTag("g4SimHits","EcalTBH4BeamHits"), cms.InputTag("g4SimHits","HcalTB06BeamHits"), + cms.InputTag("g4SimHits","HcalHits",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","ZDCHITS",cms.InputTag.skipCurrentProcess()) + ) + +mix.mixObjects.mixTracks.input = cms.VInputTag( + cms.InputTag("g4SimHits","",cms.InputTag.skipCurrentProcess()) + ) + +mix.mixObjects.mixVertices.input = cms.VInputTag( + cms.InputTag("g4SimHits","",cms.InputTag.skipCurrentProcess()) + ) + +mix.mixObjects.mixSH.input = cms.VInputTag( + #cms.InputTag("g4SimHits","BSCHits"), cms.InputTag("g4SimHits","BCM1FHits"), cms.InputTag("g4SimHits","PLTHits"), cms.InputTag("g4SimHits","FP420SI"), + cms.InputTag("g4SimHits","MuonCSCHits",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","MuonDTHits",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","MuonRPCHits",cms.InputTag.skipCurrentProcess()), + #cms.InputTag("g4SimHits","TotemHitsRP"), cms.InputTag("g4SimHits","TotemHitsT1"), cms.InputTag("g4SimHits","TotemHitsT2Gem"), + cms.InputTag("g4SimHits","TrackerHitsPixelBarrelHighTof",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","TrackerHitsPixelBarrelLowTof",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","TrackerHitsPixelEndcapHighTof",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","TrackerHitsPixelEndcapLowTof",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","TrackerHitsTECHighTof",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","TrackerHitsTECLowTof",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","TrackerHitsTIBHighTof",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","TrackerHitsTIBLowTof",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","TrackerHitsTIDHighTof",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","TrackerHitsTIDLowTof",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","TrackerHitsTOBHighTof",cms.InputTag.skipCurrentProcess()), + cms.InputTag("g4SimHits","TrackerHitsTOBLowTof",cms.InputTag.skipCurrentProcess()) + ) + +mix.mixObjects.mixHepMC.input = cms.VInputTag( + cms.InputTag("generatorSmeared","",cms.InputTag.skipCurrentProcess()), + cms.InputTag("generator","unsmeared",cms.InputTag.skipCurrentProcess()), + cms.InputTag("generator","",cms.InputTag.skipCurrentProcess()) + ) + +mix.mixObjects.mixHepMC.makeCrossingFrame = True + +''' + +mixGen = cms.Sequence(mix)