ATLAS Offline Software
Loading...
Searching...
No Matches
Pythia8B_i Class Reference

Authors: James Catmore and Maria Smizanska James.nosp@m..Cat.nosp@m.more@.nosp@m.cern.nosp@m..ch / Maria.nosp@m..Smi.nosp@m.zansk.nosp@m.a@ce.nosp@m.rn.ch Inherits from Pythia8_i by James Monk. More...

#include <Pythia8B_i.h>

Inheritance diagram for Pythia8B_i:
Collaboration diagram for Pythia8B_i:

Public Member Functions

 Pythia8B_i (const std::string &name, ISvcLocator *pSvcLocator)
 ~Pythia8B_i ()
virtual StatusCode genInitialize ()
 For initializing the generator, if required.
virtual StatusCode genuserInitialize ()
 For initialization of user code, if required. Called after genInitialize.
virtual StatusCode callGenerator ()
 For calling the generator on each iteration of the event loop.
virtual StatusCode genFinalize ()
 For finalising the generator, if required.
virtual StatusCode fillEvt (HepMC::GenEvent *)
 For filling the HepMC event object.
bool leptonSelect (Pythia8::Event &, const std::vector< double > &, double, const std::vector< int > &, int, double, bool)
bool cleanUndecayed (Pythia8::Event &, const std::vector< int > &)
bool pairProperties (Pythia8::Event &, const std::vector< int > &, double, bool)
void descendThroughDecay (Pythia8::Event &, std::vector< Pythia8::Particle > &, int) const
std::vector< int > getCodes (const std::vector< Pythia8::Particle > &) const
bool compare (std::vector< int >, std::vector< int >) const
bool passesPTCuts (const std::vector< Pythia8::Particle > &) const
bool passesEtaCuts (const std::vector< Pythia8::Particle > &) const
bool signalAccept (Pythia8::Event &, const std::vector< int > &, unsigned int) const
bool userSelection (Pythia8::Event &, const std::string &, const std::vector< double > &)
void printSignalSelections (const std::vector< int > &, const std::vector< double > &, const std::vector< double > &, unsigned int) const
virtual StatusCode fillWeights (HepMC::GenEvent *evt)
double pythiaVersion () const
virtual StatusCode sysInitialize () override
 Override sysInitialize.
virtual bool isClonable () const override
 Specify if the algorithm is clonable.
virtual StatusCode sysExecute (const EventContext &ctx) override
 Execute an algorithm.
virtual const DataObjIDColl & extraOutputDeps () const override
 Return the list of extra output dependencies.
virtual bool filterPassed (const EventContext &ctx) const
 Get filter decision:
virtual void setFilterPassed (bool state, const EventContext &ctx) const
 Set filter decision:
ServiceHandle< StoreGateSvc > & evtStore ()
 The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
const ServiceHandle< StoreGateSvc > & detStore () const
 The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
virtual StatusCode sysStart () override
 Handle START transition.
virtual std::vector< Gaudi::DataHandle * > inputHandles () const override
 Return this algorithm's input handles.
virtual std::vector< Gaudi::DataHandle * > outputHandles () const override
 Return this algorithm's output handles.
Gaudi::Details::PropertyBase & declareProperty (Gaudi::Property< T, V, H > &t)
void updateVHKA (Gaudi::Details::PropertyBase &)
MsgStream & msg () const
bool msgLvl (const MSG::Level lvl) const
Event loop algorithm methods: not to be overloaded
StatusCode initialize ()
StatusCode execute (const EventContext &ctx)
 Execute method.
StatusCode finalize ()
Event collection accessors (const and non-const)
HepMC::GenEvent *event ATLAS_NOT_CONST_THREAD_SAFE ()
 Access the current signal event (first in the McEventCollection).
McEventCollection *events ATLAS_NOT_CONST_THREAD_SAFE ()
 Access the current event's McEventCollection.
const HepMC::GenEventevent_const (const EventContext &ctx) const
 Access the current signal event (const).
const McEventCollectionevents_const (const EventContext &ctx) const
 Access the current event's McEventCollection (const).

Static Public Member Functions

static const std::string & pythia_stream ()
static std::string xmlpath ()

Static Public Attributes

static CLHEP::HepRandomEngine * p_rndmEngine = nullptr

Protected Member Functions

bool useRndmGenSvc () const
bool useReseed () const
virtual bool isReEntrant () const override final
 Legacy algorithms are not thread-safe.
void renounceArray (SG::VarHandleKeyArray &handlesArray)
 remove all handles from I/O resolution
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce (T &h)
void extraDeps_update_handler (Gaudi::Details::PropertyBase &ExtraDeps)
 Add StoreName to extra input/output deps as needed.

Protected Attributes

std::unique_ptr< Pythia8::Pythia > m_pythia {}
HepMC::Pythia8ToHepMC m_pythiaToHepMC
UnsignedIntegerProperty m_maxFailures {this, "MaxFailures", 10}
BooleanProperty m_useRndmGenSvc {this, "useRndmGenSvc", true, "the max number of consecutive failures"}
std::shared_ptr< customRndmm_atlasRndmEngine {}
BooleanProperty m_useReseed {this,"useReseed", false}
IntegerProperty m_dsid {this, "Dsid", 999999, "Dataset ID number"}
StringArrayProperty m_userHooks {this, "UserHooks", {} }
DoubleProperty m_pt0timesMPI {this,"pT0timesMPI", 1.0}
DoubleProperty m_numberAlphaS {this,"numberAlphaS", 3.0}
BooleanProperty m_sameAlphaSAsMPI {this,"useSameAlphaSasMPI", false}

Private Types

enum  PDGID {
  PROTON =2212 , ANTIPROTON =-2212 , LEAD =1000822080 , OXYGEN =1000080160 ,
  HELIUM = 1000020040 , NEUTRON =2112 , ANTINEUTRON =-2112 , MUON =13 ,
  ANTIMUON =-13 , ELECTRON =11 , POSITRON =-11 , INVALID =0
}
typedef ServiceHandle< StoreGateSvcStoreGateSvc_t

Private Member Functions

Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
 specialization for handling Gaudi::Property<SG::VarHandleKey>

Static Private Member Functions

static std::string findValue (const std::string &command, const std::string &key)
static int s_allowedTunes (double version)

Private Attributes

unsigned int m_had
unsigned int m_dec
int m_trigCode
unsigned int m_nSignalRequired
std::vector< int > m_bcodes
std::vector< int > m_sigCodes
std::vector< int > m_cutCount
std::vector< double > m_trigPtCut
std::vector< double > m_sigPtCuts
std::vector< double > m_sigEtaCuts
int m_totalPythiaCalls
int m_totalBQuark
int m_totalCQuark
int m_totalBBarQuark
int m_totalCBarQuark
int m_totalClone
int m_passingTriggerCuts
int m_internal_event_number
int m_totalHard
int m_atLeastOneAcc
std::map< int, int > m_speciesCount
std::vector< long int > m_seeds
bool m_and
bool m_oppCharges
bool m_vetoDoubleB
bool m_vetoDoubleC
bool m_selectBQuarks
bool m_selectCQuarks
double m_qPtCut
double m_aqPtCut
double m_qEtaCut
double m_aqEtaCut
double m_invMass
double m_trigEtaCut
std::string m_userString
std::vector< double > m_userVar
std::vector< Pythia8::Event > m_BEventBuffer
std::vector< int > m_internalEventNumbers
bool m_doSuppressSmallPT
Pythia8::SuppressSmallPT * m_SuppressSmallPT
unsigned int m_failureCount
double m_version {-1.}
StringArrayProperty m_commands
std::vector< std::string > m_userParams
std::vector< std::string > m_userModes
DoubleProperty m_collisionEnergy {this, "CollisionEnergy", 14000.0}
StringProperty m_beam1 {this, "Beam1", "PROTON"}
StringProperty m_beam2 {this, "Beam2", "PROTON"}
bool m_override_transform_beamenergy {false}
StringProperty m_lheFile {this, "LHEFile", ""}
BooleanProperty m_doCKKWLAcceptance {this, "CKKWLAcceptance", false}
BooleanProperty m_doFxFxXS {this, "FxFxXS", false}
BooleanProperty m_computeEfficiency {this, "computeEfficiency", false}
double m_nAccepted {0.}
double m_nMerged {0.}
double m_sigmaTotal {0.}
double m_conversion {1.}
std::map< std::string, PDGIDm_particleIDs
StringProperty m_userProcess {this, "UserProcess", ""}
std::shared_ptr< Pythia8::Sigma2Process > m_procPtr {}
std::vector< UserHooksPtrTypem_userHooksPtrs {}
StringProperty m_userResonances {this, "UserResonances", ""}
std::vector< std::shared_ptr< Pythia8::ResonanceWidths > > m_userResonancePtrs
BooleanProperty m_useLHAPDF {this, "UseLHAPDF", true}
StringProperty m_particleDataFile {this, "ParticleData", ""}
StringProperty m_outputParticleDataFile {this, "OutputParticleData", "ParticleData.local.xml"}
double m_mergingWeight {1.0}
double m_enhanceWeight {1.0}
std::vector< std::string > m_weightIDs {}
std::vector< std::string > m_weightNames {}
bool m_doLHE3Weights {false}
std::vector< std::string > m_weightCommands {}
std::vector< std::string > m_showerWeightNames {"nominal"}
StringArrayProperty m_showerWeightNamesProp {this, "ShowerWeightNames", {} }
PublicToolHandle< IPythia8Customm_athenaTool {this, "CustomInterface", ""}
BooleanProperty m_saveLHE {this, "SaveLHERecord", false}
DataObjIDColl m_extendedExtraObjects
 Extra output dependency collection, extended by AthAlgorithmDHUpdate to add symlinks.
StoreGateSvc_t m_evtStore
 Pointer to StoreGate (event store by default).
StoreGateSvc_t m_detStore
 Pointer to StoreGate (detector store by default).
std::vector< SG::VarHandleKeyArray * > m_vhka
bool m_varHandleArraysDeclared

Features for derived classes to use internally

ServiceHandle< IAthRNGSvcm_rndmSvc {this, "RndmSvc", "AthRNGSvc"}
 Data members.
ServiceHandle< IIncidentSvc > m_incidentSvc {this, "IncidentSvc", "IncidentSvc"}
 Handle on the incident service.
IntegerProperty m_randomSeed {this, "RandomSeed", 1234567, "Random seed for the built-in random engine"}
 Seed for random number engine.
BooleanProperty m_isAfterburner {this, "IsAfterburner", false, "Set true if generator modifies existing events rather than creating new ones"}
 Flag for normal vs. afterburner generators.
std::shared_ptr< HepMC3::GenRunInfo > m_runinfo {}
 The run info for HepMC3.
CLHEP::HepRandomEngine * getRandomEngine (const std::string &streamName, const EventContext &ctx) const
CLHEP::HepRandomEngine * getRandomEngine (const std::string &streamName, unsigned long int randomSeedOffset, const EventContext &ctx) const
CLHEP::HepRandomEngine * getRandomEngineDuringInitialize (const std::string &streamName, unsigned long int randomSeedOffset, unsigned int conditionsRun=1, unsigned int lbn=1) const

Properties

SG::ReadHandleKey< McEventCollectionm_mcevents_const { this, "McEventKey", "GEN_EVENT", "StoreGate key of the MC event collection" }
 Const handle to the MC event collection.
std::string m_mcEventKey {}
 StoreGate key for the MC event collection (defaults to GEN_EVENT).
BooleanProperty m_mkMcEvent {this, "MakeMcEvent", false, "Create a new MC event collection if it doesn't exist"}
 Flag to determine if a new MC event collection should be made if it doesn't exist.

Detailed Description

Authors: James Catmore and Maria Smizanska James.nosp@m..Cat.nosp@m.more@.nosp@m.cern.nosp@m..ch / Maria.nosp@m..Smi.nosp@m.zansk.nosp@m.a@ce.nosp@m.rn.ch Inherits from Pythia8_i by James Monk.

Definition at line 16 of file Pythia8B_i.h.

Member Typedef Documentation

◆ StoreGateSvc_t

typedef ServiceHandle<StoreGateSvc> AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::StoreGateSvc_t
privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Member Enumeration Documentation

◆ PDGID

enum Pythia8_i::PDGID
privateinherited
Enumerator
PROTON 
ANTIPROTON 
LEAD 
OXYGEN 
HELIUM 
NEUTRON 
ANTINEUTRON 
MUON 
ANTIMUON 
ELECTRON 
POSITRON 
INVALID 

Definition at line 124 of file Pythia8_i.h.

124{PROTON=2212, ANTIPROTON=-2212, LEAD=1000822080, OXYGEN=1000080160, HELIUM= 1000020040, NEUTRON=2112, ANTINEUTRON=-2112, MUON=13, ANTIMUON=-13, ELECTRON=11, POSITRON=-11, INVALID=0};

Constructor & Destructor Documentation

◆ Pythia8B_i()

Pythia8B_i::Pythia8B_i ( const std::string & name,
ISvcLocator * pSvcLocator )

Definition at line 27 of file Pythia8B_i.cxx.

28 : Pythia8_i(name,pSvcLocator) {
29 declareProperty("NHadronizationLoops", m_had=1);
30 declareProperty("NDecayLoops", m_dec=1);
31 declareProperty("SelectBQuarks",m_selectBQuarks=true);
32 declareProperty("SelectCQuarks",m_selectCQuarks=false);
33 declareProperty("QuarkPtCut", m_qPtCut= 0.0);
34 declareProperty("AntiQuarkPtCut", m_aqPtCut= 0.0);
35 declareProperty("QuarkEtaCut", m_qEtaCut=102.5);
36 declareProperty("AntiQuarkEtaCut", m_aqEtaCut=102.5);
37 declareProperty("RequireBothQuarksPassCuts", m_and=false);
38 declareProperty("VetoDoubleBEvents", m_vetoDoubleB=true);
39 declareProperty("VetoDoubleCEvents", m_vetoDoubleC=true);
40 declareProperty("BPDGCodes", m_bcodes);
41 declareProperty("TriggerPDGCode", m_trigCode=0); // 0 = no selection applied
42 declareProperty("TriggerStatePtCut", m_trigPtCut);
43 declareProperty("TriggerStateEtaCut", m_trigEtaCut=2.5);
44 declareProperty("MinimumCountPerCut", m_cutCount);
45 declareProperty("TriggerInvariantMassCut",m_invMass=0.0);
46 declareProperty("TriggerOppositeCharges",m_oppCharges=false);
47 declareProperty("SignalPDGCodes", m_sigCodes);
48 declareProperty("SignalPtCuts", m_sigPtCuts);
49 declareProperty("SignalEtaCuts", m_sigEtaCuts);
50 declareProperty("NumberOfSignalsRequiredPerEvent", m_nSignalRequired=1);
51 declareProperty("UserSelection", m_userString="NONE");
52 declareProperty("UserSelectionVariables", m_userVar);
53 declareProperty("SuppressSmallPT", m_doSuppressSmallPT=false);
54
55 m_totalBQuark = 0;
57 m_totalCQuark = 0;
60 m_totalHard = 0;
61 m_totalClone = 0;
65 m_speciesCount.clear();
68
69 for (std::vector<int>::iterator iit=m_bcodes.begin(); iit!=m_bcodes.end(); ++iit) {
70 m_speciesCount[*iit] = 0;
71 }
72 m_runinfo = std::make_shared<HepMC3::GenRunInfo>();
73 std::vector<std::string> names = {"Default"};
74 m_runinfo->set_weight_names(names);
75}
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
std::shared_ptr< HepMC3::GenRunInfo > m_runinfo
The run info for HepMC3.
Definition GenModule.h:90
int m_totalClone
Definition Pythia8B_i.h:51
std::vector< int > m_sigCodes
Definition Pythia8B_i.h:48
bool m_selectCQuarks
Definition Pythia8B_i.h:54
bool m_doSuppressSmallPT
Definition Pythia8B_i.h:60
int m_totalBQuark
Definition Pythia8B_i.h:51
int m_totalHard
Definition Pythia8B_i.h:51
double m_aqPtCut
Definition Pythia8B_i.h:55
std::vector< double > m_sigPtCuts
Definition Pythia8B_i.h:50
double m_aqEtaCut
Definition Pythia8B_i.h:55
unsigned int m_nSignalRequired
Definition Pythia8B_i.h:47
int m_totalCQuark
Definition Pythia8B_i.h:51
Pythia8::SuppressSmallPT * m_SuppressSmallPT
Definition Pythia8B_i.h:61
unsigned int m_had
Definition Pythia8B_i.h:45
int m_atLeastOneAcc
Definition Pythia8B_i.h:51
int m_trigCode
Definition Pythia8B_i.h:46
int m_totalBBarQuark
Definition Pythia8B_i.h:51
bool m_vetoDoubleC
Definition Pythia8B_i.h:54
bool m_oppCharges
Definition Pythia8B_i.h:54
unsigned int m_dec
Definition Pythia8B_i.h:45
std::vector< int > m_bcodes
Definition Pythia8B_i.h:48
double m_qEtaCut
Definition Pythia8B_i.h:55
std::map< int, int > m_speciesCount
Definition Pythia8B_i.h:52
double m_qPtCut
Definition Pythia8B_i.h:55
std::string m_userString
Definition Pythia8B_i.h:56
bool m_vetoDoubleB
Definition Pythia8B_i.h:54
std::vector< int > m_cutCount
Definition Pythia8B_i.h:48
bool m_selectBQuarks
Definition Pythia8B_i.h:54
double m_invMass
Definition Pythia8B_i.h:55
unsigned int m_failureCount
Definition Pythia8B_i.h:62
std::vector< double > m_trigPtCut
Definition Pythia8B_i.h:49
std::vector< double > m_userVar
Definition Pythia8B_i.h:57
int m_totalPythiaCalls
Definition Pythia8B_i.h:51
std::vector< double > m_sigEtaCuts
Definition Pythia8B_i.h:50
int m_passingTriggerCuts
Definition Pythia8B_i.h:51
int m_totalCBarQuark
Definition Pythia8B_i.h:51
int m_internal_event_number
Definition Pythia8B_i.h:51
double m_trigEtaCut
Definition Pythia8B_i.h:55
Pythia8_i(const std::string &name, ISvcLocator *pSvcLocator)
Definition Pythia8_i.cxx:68

◆ ~Pythia8B_i()

Pythia8B_i::~Pythia8B_i ( )

Definition at line 78 of file Pythia8B_i.cxx.

78 {
79}

Member Function Documentation

◆ ATLAS_NOT_CONST_THREAD_SAFE() [1/2]

McEventCollection *events GenBase::ATLAS_NOT_CONST_THREAD_SAFE ( )
inherited

Access the current event's McEventCollection.

Note
This function will make a new McEventCollection if there is not already a valid one and MakeMcEvent=True.

◆ ATLAS_NOT_CONST_THREAD_SAFE() [2/2]

HepMC::GenEvent *event GenBase::ATLAS_NOT_CONST_THREAD_SAFE ( )
inlineinherited

Access the current signal event (first in the McEventCollection).

Note
This function will make a new McEventCollection if there is not already a valid one and MakeMcEvent=True.

Definition at line 74 of file GenBase.h.

74 {
75 if (events()->empty())
76 ATH_MSG_ERROR("McEventCollection is empty during first event access");
77 return *(events()->begin());
78 }
#define ATH_MSG_ERROR(x)
static const Attributes_t empty

◆ callGenerator()

StatusCode Pythia8B_i::callGenerator ( )
virtual

For calling the generator on each iteration of the event loop.

Reimplemented from Pythia8_i.

Definition at line 148 of file Pythia8B_i.cxx.

148 {
149
150 ATH_MSG_DEBUG(">>> Pythia8B_i from callGenerator");
151
152 // Check the buffers are empty - if not return so that the GenModule
153 // calls fillEvt again and takes another event from the buffer
154 if (m_BEventBuffer.size()!=0) {
155 ATH_MSG_DEBUG("BEventBuffer not empty; skipping callGenerator");
156 return StatusCode::SUCCESS;
157 }
158
160 //Re-seed the random number stream
161 long seeds[7];
162 const EventContext& ctx = Gaudi::Hive::currentContext();
163 ATHRNG::calculateSeedsMC21(seeds, "PYTHIA8B", ctx.eventID().event_number(), m_dsid, m_randomSeed);
164 m_atlasRndmEngine->getEngine()->setSeeds(seeds, 0); // NOT THREAD-SAFE
165 m_seeds.clear();
166 m_seeds.push_back(seeds[0]);
167 m_seeds.push_back(seeds[1]);
168 }
169
170 // Generator. Shorthand for event.
171 Pythia8::Event& event = Pythia8_i::m_pythia->event;
172
173 // Begin event loop.
174 int iEvent(0);
175 while (iEvent < 1) { // keep going until an event is accepted
176
177 ATH_MSG_DEBUG("Throwing the dice....");
178// if (!Pythia8_i::m_pythia->next()) continue;
179 if ( !Pythia8_i::m_pythia->next() ) {
180 // First check if it failed because it ran out of events.
181 if ( Pythia8_i::m_pythia->info.atEndOfFile() ) {
182 return StatusCode::FAILURE;
183 }
184
185 // Otherwise, just make sure that it's not failing too many times in a row.
187 if ( m_failureCount >= m_maxFailures ) {
188 ATH_MSG_ERROR("Exceeded the max number of consecutive event failures.");
189 return StatusCode::FAILURE;
190 } else {
191 ATH_MSG_INFO("Event generation failed - re-trying.");
192 continue;
193 }
194 }
195 // Reset failure counter and increment event counter after a successful event
196 m_failureCount = 0;
198
199 // Find b(c)/antib(c) quarks and enforce cuts as required
200 int nbBeforeSelection(0);
201 int nbbarBeforeSelection(0);
202 int ncBeforeSelection(0);
203 int ncbarBeforeSelection(0);
204 int nbQuark(0);
205 int nbbarQuark(0);
206 int ncQuark(0);
207 int ncbarQuark(0);
208 int stat;
209 for (int i = 0; i < event.size(); ++i) {
210 stat = event[i].statusAbs();
211 bool isBQuark(false);
212 bool isCQuark(false);
213 bool isAntiBQuark(false);
214 bool isAntiCQuark(false);
215 // b-quarks (=anti-B hadrons)
216 if (event[i].id() == 5 && (stat == 62 || stat == 63)) {isBQuark=true; ++nbBeforeSelection;}
217 if (event[i].id() == 4 && (stat == 62 || stat == 63)) {isCQuark=true; ++ncBeforeSelection;}
218 if ( (isBQuark && m_selectBQuarks) || ((isCQuark && m_selectCQuarks)) ) {
219 bool passesPtCut(false); bool passesEtaCut(false);
220 std::string accString = " : REJECTED";
221 double qpt = event[i].pT(); double qeta = std::abs(event[i].eta());
222 if (qpt>m_qPtCut) passesPtCut=true;
223 if (qeta<m_qEtaCut) passesEtaCut=true;
224 if (passesPtCut && passesEtaCut) {
225 if (isBQuark) ++nbQuark;
226 if (isCQuark) ++ncQuark;
227 accString = " : ACCEPTED";
228 }
229 if (isBQuark) ATH_MSG_DEBUG("bQuark pt/eta = " << qpt << "/" << qeta << accString);
230 if (isCQuark) ATH_MSG_DEBUG("cQuark pt/eta = " << qpt << "/" << qeta << accString);
231 }
232 // anti-b quarks (=B-hadrons)
233 if (event[i].id() == -5 && (stat == 62 || stat == 63)) {isAntiBQuark=true; ++nbbarBeforeSelection;}
234 if (event[i].id() == -4 && (stat == 62 || stat == 63)) {isAntiCQuark=true; ++ncbarBeforeSelection;}
235 if ( (isAntiBQuark && m_selectBQuarks) || ((isAntiCQuark && m_selectCQuarks)) ) {
236 bool passesPtCut(false); bool passesEtaCut(false);
237 std::string accString = " : REJECTED";
238 double aqpt = event[i].pT(); double aqeta = std::abs(event[i].eta());
239 if (aqpt>m_aqPtCut) passesPtCut=true;
240 if (aqeta<m_aqEtaCut) passesEtaCut=true;
241 if (passesPtCut && passesEtaCut) {
242 if (isAntiBQuark) ++nbbarQuark;
243 if (isAntiCQuark) ++ncbarQuark;
244 accString = " : ACCEPTED";
245 }
246 if (isAntiBQuark) ATH_MSG_DEBUG("bbarQuark pt/eta = " << aqpt << "/" << aqeta << accString);
247 if (isAntiCQuark) ATH_MSG_DEBUG("ccbarQuark pt/eta = " << aqpt << "/" << aqeta << accString);
248 }
249 }
250 if (nbBeforeSelection+nbbarBeforeSelection>=4 && m_vetoDoubleB) {
251 ATH_MSG_DEBUG("At user request, rejecting double b-bbar event; throwing dice again");
252 continue;
253 }
254 if (ncBeforeSelection+ncbarBeforeSelection>=4 && m_vetoDoubleC) {
255 ATH_MSG_DEBUG("At user request, rejecting double c-cbar event; throwing dice again");
256 continue;
257 }
258 bool rejectBBbar(false);
259 bool rejectCCbar(false);
260 // Enforce user selections
261 if (!m_and && (nbbarQuark==0 && nbQuark==0) && m_selectBQuarks) rejectBBbar=true;
262 if (m_and && (nbbarQuark==0 || nbQuark==0) && m_selectBQuarks) rejectBBbar=true;
263 if (!m_and && (ncbarQuark==0 && ncQuark==0) && m_selectCQuarks) rejectCCbar=true;
264 if (m_and && (ncbarQuark==0 || ncQuark==0) && m_selectCQuarks) rejectCCbar=true;
265 if (m_selectBQuarks && m_selectCQuarks && rejectBBbar && rejectCCbar) {
266 ATH_MSG_DEBUG("No b- or c- quarks accepted; throwing the dice again");
267 continue;
268 }
269 if (m_selectBQuarks && !m_selectCQuarks && rejectBBbar) {
270 ATH_MSG_DEBUG("No b-quarks accepted; throwing the dice again");
271 continue;
272 }
273 if (!m_selectBQuarks && m_selectCQuarks && rejectCCbar) {
274 ATH_MSG_DEBUG("No c-quarks accepted; throwing the dice again");
275 continue;
276 }
277 m_totalBQuark += nbQuark;
278 m_totalBBarQuark += nbbarQuark;
279 m_totalCQuark += ncQuark;
280 m_totalCBarQuark += ncbarQuark;
281 ++m_totalHard;
282
283 // Switch off decays of species of interest if requested (only for repeated decays)
284 bool doRepeatedDecays(false); if (m_dec>1) doRepeatedDecays=true;
285 if (doRepeatedDecays) {
286 for (unsigned int iC = 0; iC < m_bcodes.size(); ++iC) Pythia8_i::m_pythia->particleData.mayDecay( m_bcodes[iC], false);
287 }
288
289 // REPEATED HADRONIZATION LOOP
290 // Make a spare copy of event
291 Pythia8::Event eventCopy;
292 std::vector<Pythia8::Event> repeatHadronizedEvents;
293 std::vector<Pythia8::Event>::iterator eventItr;
294 eventCopy = event;
295 // Begin loop over rounds of hadronization for same event.
296 for (unsigned int iRepeat = 0; iRepeat < m_had; ++iRepeat) {
297 ++m_totalClone;
298 // Repeated hadronization: restore saved event record.
299 if (iRepeat > 0) event = eventCopy;
300 // Repeated hadronization: do HadronLevel (repeatedly).
301 if (!Pythia8_i::m_pythia->forceHadronLevel()) continue;
302 repeatHadronizedEvents.push_back(event); // save the events for selections
303 }
304
305 // REPEATED DECAY LOOP
306 std::vector<Pythia8::Event> savedEvents;
307 if (doRepeatedDecays) {
308 // switch decays back on
309 for (unsigned int iC = 0; iC < m_bcodes.size(); ++iC) Pythia8_i::m_pythia->particleData.mayDecay( m_bcodes[iC], true);
310 // Loop over repeatedly hadronized events
311 for (eventItr=repeatHadronizedEvents.begin(); eventItr!=repeatHadronizedEvents.end(); ++eventItr) {
312 for (unsigned int iRepeat = 0; iRepeat < m_dec; ++iRepeat) {
313 event = (*eventItr);
314 if (!Pythia8_i::m_pythia->moreDecays()) break;
315 savedEvents.push_back(event);
316 }
317 }
318 }
319 if (!doRepeatedDecays) savedEvents = std::move(repeatHadronizedEvents);
320
321 //event.list();
322
323 // Make final selections on savedEvents
324 std::vector<Pythia8::Event> finalEvents;
325 bool signalSelect(false);
326 if (m_sigCodes.size()>0) signalSelect=true;
327 for (eventItr=savedEvents.begin(); eventItr!=savedEvents.end(); ++eventItr) {
328 // Does event contain basic leptons needed to fire triggers?
331 // Does the events contain a signal decay as required by the user
332 if (signalSelect) {
333 if (!cleanUndecayed(*eventItr,m_bcodes)) continue;
334 if (!signalAccept(*eventItr,m_sigCodes,m_nSignalRequired)) continue;
335 }
336 // User-defined selections if any
337 if (!m_userString.empty()) {
338 if (!userSelection(*eventItr,m_userString,m_userVar)) continue;
339 }
340 finalEvents.push_back(*eventItr);
341 ++iEvent;
342 }
343 savedEvents.clear();
344 if (finalEvents.size()>0) ++m_atLeastOneAcc;
345
346 // Species counting and preparing vector for filling
347 for (eventItr=finalEvents.begin(); eventItr!=finalEvents.end(); ++eventItr) {
348 for (int i = 0; i < (*eventItr).size(); ++i) {
349 int id = (*eventItr)[i].id();
350 std::map<int,int>::iterator it;
351 it = m_speciesCount.find(id);
352 if (it != m_speciesCount.end()) {
353 it->second++;
354 }
355 }
356 m_BEventBuffer.push_back(*eventItr);
359 //(*eventItr).list();
360 }
361
362 // End of event loop.
363 }
364
365 StatusCode returnCode = StatusCode::SUCCESS;
366
367 return returnCode;
368}
Scalar eta() const
pseudorapidity method
#define ATH_MSG_INFO(x)
#define ATH_MSG_DEBUG(x)
IntegerProperty m_randomSeed
Seed for random number engine.
Definition GenModule.h:84
std::vector< long int > m_seeds
Definition Pythia8B_i.h:53
static CLHEP::HepRandomEngine * p_rndmEngine
Definition Pythia8B_i.h:41
std::vector< int > m_internalEventNumbers
Definition Pythia8B_i.h:59
bool leptonSelect(Pythia8::Event &, const std::vector< double > &, double, const std::vector< int > &, int, double, bool)
std::vector< Pythia8::Event > m_BEventBuffer
Definition Pythia8B_i.h:58
bool signalAccept(Pythia8::Event &, const std::vector< int > &, unsigned int) const
bool cleanUndecayed(Pythia8::Event &, const std::vector< int > &)
bool userSelection(Pythia8::Event &, const std::string &, const std::vector< double > &)
IntegerProperty m_dsid
Definition Pythia8_i.h:97
std::shared_ptr< customRndm > m_atlasRndmEngine
Definition Pythia8_i.h:93
std::unique_ptr< Pythia8::Pythia > m_pythia
Definition Pythia8_i.h:88
UnsignedIntegerProperty m_maxFailures
Definition Pythia8_i.h:90
BooleanProperty m_useRndmGenSvc
Definition Pythia8_i.h:92
void calculateSeedsMC21(long *seeds, const std::string &algName, uint64_t ev, uint64_t run, uint64_t offset=0)
Set the random seed using a string (e.g.
::StatusCode StatusCode
StatusCode definition for legacy code.

◆ cleanUndecayed()

bool Pythia8B_i::cleanUndecayed ( Pythia8::Event & theEvent,
const std::vector< int > & bCodes )

Definition at line 539 of file Pythia8B_i.cxx.

539 {
540
541 bool cleanEvent(true);
542 std::string accString(" : ACCEPTED");
543 for (int i = 0; i<theEvent.size(); ++i) {
544 const Pythia8::Particle &theParticle = theEvent[i];
545 int id = theParticle.id();
546 int status = theParticle.status();
547 for (auto iItr = bCodes.begin(); iItr!=bCodes.end(); ++iItr) {
548 if ( (id == *iItr) && (status>0) ) {accString=" : REJECTED"; cleanEvent = false;}
549 }
550 }
551 ATH_MSG_DEBUG("Check event for undecayed signal particles" << accString);
552
553 return cleanEvent;
554
555}
status
Definition merge.py:16

◆ compare()

bool Pythia8B_i::compare ( std::vector< int > vect1,
std::vector< int > vect2 ) const

Definition at line 631 of file Pythia8B_i.cxx.

631 {
632
633 if (vect1.size()!=vect2.size()) return false;
634
635 bool isTheSame(true);
636 int size = vect1.size();
637 std::sort(vect1.rbegin(), vect1.rend());
638 std::sort(vect2.rbegin(), vect2.rend());
639 for (int i=0; i<size; ++i) {
640 if (vect1[i] != vect2[i]) {isTheSame = false; break;}
641 }
642 return isTheSame;
643
644}
size_t size() const
Number of registered mappings.
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.

◆ declareGaudiProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::declareGaudiProperty ( Gaudi::Property< T, V, H > & hndl,
const SG::VarHandleKeyType &  )
inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKey>

Definition at line 156 of file AthCommonDataStore.h.

158 {
160 hndl.value(),
161 hndl.documentation());
162
163 }

◆ declareProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::declareProperty ( Gaudi::Property< T, V, H > & t)
inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

145 {
146 typedef typename SG::HandleClassifier<T>::type htype;
148 }
Gaudi::Details::PropertyBase & declareGaudiProperty(Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
specialization for handling Gaudi::Property<SG::VarHandleKey>

◆ descendThroughDecay()

void Pythia8B_i::descendThroughDecay ( Pythia8::Event & theEvent,
std::vector< Pythia8::Particle > & list,
int i ) const

Definition at line 602 of file Pythia8B_i.cxx.

602 {
603
604 list.push_back(theEvent[i]);
605 std::vector<int> childrenIndices = theEvent.daughterList(i);
606 std::vector<int>::iterator iItr;
607 for (iItr = childrenIndices.begin(); iItr != childrenIndices.end(); ++iItr) {
608 descendThroughDecay(theEvent,list,*iItr);
609 }
610
611 return;
612
613}
void descendThroughDecay(Pythia8::Event &, std::vector< Pythia8::Particle > &, int) const
list(name, path='/')
Definition histSizes.py:38

◆ detStore()

const ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::detStore ( ) const
inlineinherited

The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 95 of file AthCommonDataStore.h.

◆ event_const()

const HepMC::GenEvent * GenBase::event_const ( const EventContext & ctx) const
inlineinherited

Access the current signal event (const).

Definition at line 81 of file GenBase.h.

81 {
82 const McEventCollection* coll = events_const(ctx);
83 if (coll->empty())
84 ATH_MSG_ERROR("Const McEventCollection is empty during first event access");
85 return *(coll->begin());
86 }
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
bool empty() const noexcept
Returns true if the collection is empty.
const McEventCollection * events_const(const EventContext &ctx) const
Access the current event's McEventCollection (const).
Definition GenBase.h:95

◆ events_const()

const McEventCollection * GenBase::events_const ( const EventContext & ctx) const
inlineinherited

Access the current event's McEventCollection (const).

Definition at line 95 of file GenBase.h.

95 {
96 SG::ReadHandle<McEventCollection> ret = SG::makeHandle(m_mcevents_const, ctx);
97 if (!ret.isValid())
98 ATH_MSG_ERROR("No McEventCollection found in StoreGate with key " << m_mcevents_const.key());
99 return ret.cptr();
100 }
SG::ReadHandleKey< McEventCollection > m_mcevents_const
Const handle to the MC event collection.
Definition GenBase.h:119
virtual bool isValid() override final
Can the handle be successfully dereferenced?
const_pointer_type cptr()
Dereference the pointer.
SG::ReadCondHandle< T > makeHandle(const SG::ReadCondHandleKey< T > &key, const EventContext &ctx=Gaudi::Hive::currentContext())

◆ evtStore()

ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::evtStore ( )
inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 85 of file AthCommonDataStore.h.

◆ execute()

StatusCode GenModule::execute ( const EventContext & ctx)
virtualinherited

Execute method.

Todo
Remove hard-coded alg name checking (already incomplete)

Reimplemented from GenBase.

Definition at line 70 of file GenModule.cxx.

70 {
71 // Examples of how to retrieve the random number engine for a given
72 // stream.
73 // NB getRandomEngine should only be called once per event for a
74 // given stream, as it causes the stream to be re-seeded each time
75 // it is called.
76
77 // Example 1 - seeded based on the current event number (+ slot, run, streamName)
78 //CLHEP::HepRandomEngine* rndmEngine = this->getRandomEngine("MyStream", ctx);
79
80 // Example 2 - seeded based on the m_randomSeed property (+ slot, run, streamName)
81 //CLHEP::HepRandomEngine* rndmEngine = this->getRandomEngine("MyStream", m_randomSeed.value(), ctx);
82
83 // Call the code that generates an event
84 CHECK(this->callGenerator());
85
86 // Create the MC event and send the GeneratorEvent stored in it to fillEvt
88 CHECK(this->fillEvt(evt));
90
91 // Add the event to the MC event collection
92 if (events()) {
93 // If this is an "afterburner" generator, replace the last event rather than add a new one
95 if (m_isAfterburner.value() || name() == "Tauola" || name() == "Photos") {
96 events()->pop_back();
97 }
98 // Add the event to the end of the collection
99 events()->push_back(evt);
100 ATH_MSG_DEBUG("MC event added to McEventCollection");
101
102 // remove the empty event in case of ParticleDecayer
103 if (name() == "ParticleDecayer") {
104 events()->pop_back();
105 }
106 }
107
108 // Call the incident service to notify that an event has been made
109 m_incidentSvc->fireIncident( Incident(name(), "McEventGenerated") );
110 return StatusCode::SUCCESS;
111}
#define CHECK(...)
Evaluate an expression and check for errors.
ServiceHandle< IIncidentSvc > m_incidentSvc
Handle on the incident service.
Definition GenModule.h:99
virtual StatusCode callGenerator()
For calling the generator on each iteration of the event loop.
Definition GenModule.h:66
virtual StatusCode fillEvt(HepMC::GenEvent *evt)=0
For filling the HepMC event object.
BooleanProperty m_isAfterburner
Flag for normal vs. afterburner generators.
Definition GenModule.h:87
void fillBarcodesAttribute(GenEvent *e)
Definition GenEvent.h:393
GenEvent * newGenEvent(const int signal_process_id, const int event_number)
Definition GenEvent.h:360
HepMC3::GenEvent GenEvent
Definition GenEvent.h:39

◆ extraDeps_update_handler()

void AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::extraDeps_update_handler ( Gaudi::Details::PropertyBase & ExtraDeps)
protectedinherited

Add StoreName to extra input/output deps as needed.

use the logic of the VarHandleKey to parse the DataObjID keys supplied via the ExtraInputs and ExtraOuputs Properties to add the StoreName if it's not explicitly given

◆ extraOutputDeps()

const DataObjIDColl & AthCommonAlgorithm< Gaudi::Algorithm >::extraOutputDeps ( ) const
overridevirtualinherited

Return the list of extra output dependencies.

This list is extended to include symlinks implied by inheritance relations.

Definition at line 89 of file AthCommonAlgorithm.cxx.

54{
55 // If we didn't find any symlinks to add, just return the collection
56 // from the base class. Otherwise, return the extended collection.
57 if (!m_extendedExtraObjects.empty()) {
59 }
61}
Common base class for algorithms.

◆ fillEvt()

StatusCode Pythia8B_i::fillEvt ( HepMC::GenEvent * evt)
virtual

For filling the HepMC event object.

Reimplemented from Pythia8_i.

Definition at line 371 of file Pythia8B_i.cxx.

371 {
372
373 ATH_MSG_DEBUG(">>> Pythia8B_i from fillEvt");
374 ATH_MSG_DEBUG("BEventBuffer contains " << m_BEventBuffer.size() << " events ");
375 if (m_BEventBuffer.size()==0 ) {
376 ATH_MSG_DEBUG("BEventBuffer now empty - going to next Pythia event");
377 return StatusCode::SUCCESS;
378 }
379
380 Pythia8::Event &pyev = *(m_BEventBuffer.begin());
381 evt->set_event_number(*(m_internalEventNumbers.begin()));
382 m_pythiaToHepMC.fill_next_event(pyev, evt, 1, &Pythia8_i::m_pythia->info);
383
384
385 // set the randomseeds
388
389 // set the event weight
390 if (!evt->run_info()) evt->set_run_info(m_runinfo);
391 evt->set_units(HepMC3::Units::MEV, HepMC3::Units::MM);
392
393 // Remove event/number from buffer
394 m_BEventBuffer.erase(m_BEventBuffer.begin());
396
397
398 return StatusCode::SUCCESS;
399}
bool useRndmGenSvc() const
Definition Pythia8_i.h:85
HepMC::Pythia8ToHepMC m_pythiaToHepMC
Definition Pythia8_i.h:89
void set_random_states(GenEvent *e, std::vector< long int > &a)
Definition GenEvent.h:588

◆ fillWeights()

StatusCode Pythia8_i::fillWeights ( HepMC::GenEvent * evt)
virtualinherited

mismatch in weight name!

Definition at line 482 of file Pythia8_i.cxx.

482 {
483
484 // reset
485 m_mergingWeight = m_pythia->info.mergingWeightNLO(); // first initialisation, good for the nominal weight in any case (see Merging:includeWeightInXsection).
486 m_enhanceWeight = 1.0; // better to keep the enhancement from UserHooks in a separate variable, to keep the code clear
487
488#ifndef PYTHIA8_304SERIES
489 // include Enhance userhook weight
490 for(const auto &hook: m_userHooksPtrs) {
491 if (hook->canEnhanceEmission()) {
492 m_enhanceWeight *= hook->getEnhancedEventWeight();
493 }
494 }
495#endif // not PYTHIA8_304SERIES
496
497 // DO NOT try to distinguish between phase space and merging contributions: only their product contains both, so always multiply them in order to be robust
498 double eventWeight = m_pythia->info.weight() * m_mergingWeight * m_enhanceWeight;
499
500 std::map<std::string,double> fWeights; // temporary weight vector, needed while managing hepmc2 and 3 simultaneously
501
502 // save as Default the usual pythia8 weight, including enhancements
503 // NB: info.weight() is built upon info.eventWeightLHEF (but includes shower "decorations"), not from the rwgt LHE vector.
504 size_t atlas_specific_weights = 1;
505 fWeights["Default"]=eventWeight;
506 if(m_internal_event_number == 1) {
507 ATH_MSG_INFO("found LHEF version "<< m_pythia->info.LHEFversion() );
508 if (m_pythia->info.LHEFversion() > 1) m_doLHE3Weights = true;
509 m_weightIDs.clear();
510 m_weightIDs.push_back("Default"); // notice that we are assigning a weightname = weightID
511 }
512
513 // Now systematic weights: a) generator and b) shower variations
514 // merging implies reading from lhe files (while the opposite is not true, ok, but still works)
515 if (m_lheFile!="") {
516 // remember that the true merging component can be disentangled only when there are LHE input events
517
518 // make sure to get the merging weight correctly: this is needed -only- for systematic weights (generator and parton shower variations)
519 // NB: most robust way, regardless of the specific value of Merging:includeWeightInXsection (if on, info.mergingWeight() is always 1.0 ! )
520 m_mergingWeight *= ( m_pythia->info.weight() / m_pythia->info.eventWeightLHEF ); // this is the actual merging component, regardless of conventions,
521 // needed by systematic variation weights
522 // this includes rare compensation or selection-biased weights
523 }
524
525
526 ATH_MSG_DEBUG("Event weights: enhancing weight, merging weight, total weight = "<<m_enhanceWeight<<", "<<m_mergingWeight<<", "<<eventWeight);
527
528 // we already filled the "Default" weight
529 std::vector<std::string>::const_iterator id = m_weightIDs.begin()+atlas_specific_weights;
530
531 // a) fill generator weights
532 if(m_pythia->info.getWeightsDetailedSize() != 0){
533 for(std::map<std::string, Pythia8::LHAwgt>::const_iterator wgt = m_pythia->info.rwgt->wgts.begin();
534 wgt != m_pythia->info.rwgt->wgts.end(); ++wgt){
535
537 m_weightIDs.push_back(wgt->first);
538 }else{
539 if(*id != wgt->first){
540 ATH_MSG_ERROR("Mismatch in LHE3 weight id. Found "<<wgt->first<<", expected "<<*id);
541 return StatusCode::FAILURE;
542 }
543 ++id;
544 }
545
546 std::map<std::string, Pythia8::LHAweight>::const_iterator weightName = m_pythia->info.init_weights->find(wgt->first);
547 if(weightName != m_pythia->info.init_weights->end()){
548 fWeights[weightName->second.contents] = m_mergingWeight * m_enhanceWeight * wgt->second.contents;
549 }else{
550 // if no name was assigned/found, assign a weightname = weightid
551 fWeights[wgt->first] = m_mergingWeight * m_enhanceWeight * wgt->second.contents;
552 }
553
554 }
555 }
556
557 // b) shower weight variations
558 // always start from second shower weight: the first was saved already as "Default"
559
560 for(int iw = 1; iw < m_pythia->info.PYTHIA8_NWEIGHTS(); ++iw){
561
562 std::string wtName = ((int)m_showerWeightNames.size() == m_pythia->info.PYTHIA8_NWEIGHTS())? m_showerWeightNames[iw]: "ShowerWt_" +std::to_string(iw);
563
564 fWeights[wtName] = m_mergingWeight * m_enhanceWeight * m_pythia->info.PYTHIA8_WEIGHT(iw)*m_conversion;
565
566 if(m_internal_event_number == 1) {
567 m_weightIDs.push_back(wtName);
568 ATH_MSG_DEBUG("Shower weight name, value, conversion: "<<wtName<<", "<< fWeights[wtName] <<","<<m_conversion );
569 }
570 }
571
572 // Now save, as a backup, also the LHEF weight without enhancements (useful for possible, rare, non-default cases)
573 // to make clear that it should not be used in analyses, save it always with a -10 factor (so that its goal is clear even if not checking its name)
574 if (m_lheFile!="")
575 {
576 fWeights["EXTRA_bare_LHE_weight"]=(-10.0)*m_pythia->info.eventWeightLHEF;
577 if(m_internal_event_number == 1) m_weightIDs.push_back("EXTRA_bare_LHE_weight");
578 }
579
580 // Sad, but needed: create a string vector with acceptable order of weight names
582 m_weightNames.clear();
583 for(const std::string &id : m_weightIDs){
584 if(m_doLHE3Weights) {
585 std::map<std::string, Pythia8::LHAweight>::const_iterator weight = m_pythia->info.init_weights->find(id);
586 if(weight != m_pythia->info.init_weights->end()) m_weightNames.push_back(weight->second.contents);
587 else m_weightNames.push_back(id);
588 }
589 else {
590 m_weightNames.push_back(id);
591 }
592 }
593 if (m_weightNames.size() != fWeights.size() ) {
594 ATH_MSG_ERROR("Something wrong when building list of weight names: " << m_weightNames.size() << " vs "<< fWeights.size() << ", exiting ...");
595 return StatusCode::FAILURE;
596 }
597 }
598
599 // The following depends on the specific hepmc2/3 implementation
600 if (!evt->run_info()) {
601 evt->set_run_info(m_runinfo);
602 }
603 evt->run_info()->set_weight_names(m_weightNames);
604
605 // for the first event, weight EXTRA_bare_LHE_weight is not present in evt->weights(), so we need to book a place for it by hand
606 if (m_internal_event_number == 1 && evt->run_info()->weight_names().size() == evt->weights().size()+1 ) {
607 evt->weights().push_back(1.0);
608 }
609
610 // added conversion GeV -> MeV to ensure correct units
611 evt->set_units(HepMC3::Units::MEV, HepMC3::Units::MM);
612
613 evt->weights().resize(fWeights.size(), 1.0);
614 for (const auto & w: fWeights) {
615 evt->weight(w.first)=w.second;
616 }
617
618 auto beams=evt->beams();
619 ATH_MSG_DEBUG( " Energy of the beams " << beams[0]->momentum().e() );
620
621//uncomment to list HepMC events
622// std::cout << " print::listing Pythia8 " << std::endl;
623// HepMC3::Print::listing(std::cout, *evt);
624
625
626 return StatusCode::SUCCESS;
627}
bool m_doLHE3Weights
Definition Pythia8_i.h:169
StringProperty m_lheFile
Definition Pythia8_i.h:133
std::vector< std::string > m_showerWeightNames
Definition Pythia8_i.h:171
std::vector< std::string > m_weightIDs
Definition Pythia8_i.h:167
std::vector< std::string > m_weightNames
Definition Pythia8_i.h:168
int m_internal_event_number
Definition Pythia8_i.h:108
double m_enhanceWeight
Definition Pythia8_i.h:166
double m_conversion
Definition Pythia8_i.h:141
double m_mergingWeight
Definition Pythia8_i.h:165
std::vector< UserHooksPtrType > m_userHooksPtrs
Definition Pythia8_i.h:154

◆ filterPassed()

virtual bool AthCommonAlgorithm< Gaudi::Algorithm >::filterPassed ( const EventContext & ctx) const
inlinevirtualinherited

Get filter decision:

Definition at line 93 of file AthCommonAlgorithm.h.

93 {
94 return execState( ctx ).filterPassed();
95 }
virtual bool filterPassed(const EventContext &ctx) const
Get filter decision:

◆ finalize()

StatusCode GenModule::finalize ( )
inlineinherited

Definition at line 55 of file GenModule.h.

55{ return genFinalize(); }
virtual StatusCode genFinalize()
For finalising the generator, if required.
Definition GenModule.h:70

◆ findValue()

std::string Pythia8_i::findValue ( const std::string & command,
const std::string & key )
staticprivateinherited

◆ genFinalize()

StatusCode Pythia8B_i::genFinalize ( )
virtual

For finalising the generator, if required.

Reimplemented from Pythia8_i.

Definition at line 403 of file Pythia8B_i.cxx.

403 {
404
405 ATH_MSG_INFO(">>> Pythia8B_i from genFinalize");
406
407 Pythia8_i::m_pythia->stat();
408
409 Pythia8::Info info = Pythia8_i::m_pythia->info;
410 double xs = info.sigmaGen();// in mb
411 xs *= 1000. * 1000.;//convert to nb
412
413 std::cout << "MetaData: cross-section (nb)= " << xs << std::endl;
414 std::cout << "Total events from Pythia = " << m_totalPythiaCalls<< std::endl;
415 std::cout << "Number of accepted b quarks = " << m_totalBQuark<< std::endl;
416 std::cout << "Number of accepted bbar quarks = " << m_totalBBarQuark<< std::endl;
417 std::cout << "Number of accepted c quarks = " << m_totalCQuark<< std::endl;
418 std::cout << "Number of accepted cbar quarks = " << m_totalCBarQuark<< std::endl;
419 std::cout << "Number of accepted b/c events before hadronization = " << m_totalHard<< std::endl;
420 std::cout << "Number of hadronization loops per b/c-event = " << m_had<< std::endl;
421 std::cout << "Total number of hadronization loops = " << m_totalClone<< std::endl;
422 std::cout << "Number of accepted b/c events yielding at least one finally accepted event = " << m_atLeastOneAcc<< std::endl;
423 std::cout << "Average number of accepted events per hard event = " << static_cast<float>(m_internal_event_number)/static_cast<float>(m_atLeastOneAcc)<< std::endl;
424 std::cout << "Number of events passing trigger cuts = " << m_passingTriggerCuts<< std::endl;
425 std::cout << "Number of accepted events = " << m_internal_event_number<< std::endl;
426 std::cout << "Summary of cuts: " << std::endl;
428 if (m_selectBQuarks) std::cout << "Quarks cuts apply to b-quarks" << std::endl;
429 if (m_selectCQuarks) std::cout << "Quarks cuts apply to c-quarks" << std::endl;
430 std::cout << "Quark pt > " << m_qPtCut << std::endl;
431 std::cout << "Antiquark pt > " << m_aqPtCut << std::endl;
432 std::cout << "Quark eta < " << m_qEtaCut << std::endl;
433 std::cout << "Antiquark eta < " << m_aqEtaCut << std::endl;
434 if (m_and) std::cout << "Require quark and anti-quark pass cuts" << std::endl;
435 }
436 std::cout << "Trigger lepton type = " << m_trigCode << std::endl;
437 std::cout << "Trigger lepton pt cuts: ";
438 for (unsigned int prCntr=0; prCntr<m_trigPtCut.size(); ++prCntr) {std::cout << m_trigPtCut[prCntr] << " ";}
439 std::cout << std::endl;
440 std::cout << "Trigger lepton eta cut: " << m_trigEtaCut << std::endl;
441 std::cout << "Required number of leptons passing each trigger cut = ";
442 for (unsigned int prCntr=0; prCntr<m_cutCount.size(); ++prCntr) {std::cout << m_cutCount[prCntr] << " ";}
443 std::cout << std::endl;
444 std::cout << "Invariant mass of trigger leptons > " << m_invMass << std::endl;
445 if (m_oppCharges) std::cout << "Trigger leptons required to have opposite charges" << std::endl;
447 if (!m_userString.empty()) std::cout << "User selection: " << m_userString << std::endl;
448
449 if (m_speciesCount.size()>0) {
450 std::cout <<"\nSpecies\tCount\n"<< std::endl;
451 for ( std::map<int,int>::iterator iter = m_speciesCount.begin();
452 iter != m_speciesCount.end(); ++iter )
453 std::cout << iter->first << '\t' << iter->second << '\n'<< std::endl;
454 }
455 if (m_dec>1) std::cout << "Number of times each of these states were copied and decayed: " << m_dec << std::endl;
456
457 float cloningFactor = (float)(m_internal_event_number)/(float)(m_atLeastOneAcc);
458 float finalXS = ((float)xs*((float)m_internal_event_number/(float)m_totalPythiaCalls))/(float)m_had;
459
460 std::cout << " I===================================================================================== " << std::endl;
461 std::cout << " I CROSSSECTION OF YOUR B-CHANNEL IS I " << std::endl;
462 std::cout << " I BX= PX*NB/AC/MHAD= I " << finalXS << " nb" << std::endl;
463 std::cout << " I I " << std::endl;
464 std::cout << " I IN CASE YOU FORCED ANY DECAY YOU SHOULD I " << std::endl;
465 std::cout << " I CORRECT CROSS SECTION BX FURTHER, MULTIPLYING I " << std::endl;
466 std::cout << " I BX BY BRANCHING RATIO(S) OF YOUR FORCED I " << std::endl;
467 std::cout << " I DECAY(S) AND BY A FACTOR OF 2 FOR SYMMETRY I " << std::endl;
468 std::cout << " I I " << std::endl;
469 std::cout << " I MORE DETAILS ON CROSS SECTION I " << std::endl;
470 std::cout << " I PYTHIA CROSS SECTION IS PX= I " << xs <<" nb" << std::endl;
471 std::cout << " I NUMBER OF ACCEPTED HARD EVENTS AC= I " << m_totalPythiaCalls << std::endl;
472 std::cout << " I NUMBER OF ACCEPTED B-EVENTS IS NB= I " << m_internal_event_number<< std::endl;
473 std::cout << " I REPEATED HADRONIZATIONS IN EACH EVENT MHAD= I " << m_had << std::endl;
474 std::cout << " I AVERAGE NUM OF ACCEPTED EVTS IN HADRONIZATION LOOP I " << cloningFactor << std::endl;
475 std::cout << " I IN CASE YOU FORCED ANY DECAY YOU SHOULD I " << std::endl;
476 std::cout << " I CORRECT CROSS SECTION BX FURTHER, MULTIPLYING I " << std::endl;
477 std::cout << " I BX BY BRANCHING RATIO(S) OF YOUR FORCED I " << std::endl;
478 std::cout << " I DECAY(S) AND BY A FACTOR OF 2 FOR SYMMETRY I " << std::endl;
479 std::cout << " I I " << std::endl;
480 std::cout << " I===================================================================================== " << std::endl;
481 std::cout << "" << std::endl;
482 std::cout << "MetaData: cross-section (nb)= " << finalXS << std::endl;
483
484 return StatusCode::SUCCESS;
485}
void printSignalSelections(const std::vector< int > &, const std::vector< double > &, const std::vector< double > &, unsigned int) const

◆ genInitialize()

StatusCode Pythia8B_i::genInitialize ( )
virtual

For initializing the generator, if required.

Here one can fill extra information, e.g. the used tools in a format generator name, version string, comment.

Reimplemented from Pythia8_i.

Definition at line 82 of file Pythia8B_i.cxx.

82 {
83
84 // Logic checks
85 unsigned int trigPtCutsSize = m_trigPtCut.size();
86 if (trigPtCutsSize!=m_cutCount.size()) {
87 ATH_MSG_ERROR("You are requesting " << trigPtCutsSize << " trigger-like pt cuts but are providing required counts for " << m_cutCount.size() << " cuts. This doesn't make sense.");
88 return StatusCode::FAILURE;
89 }
90
91 // This over-rides the genInitialize in the base class Pythia8_i, but then calls it
92 // Sets the built-in UserHook called SuppressLowPT
93 // FOR ONIA USE ONLY
94 ATH_MSG_INFO("genInitialize() from Pythia8B_i");
95
96 bool canSetHook=true;
98 Pythia8_i::m_userHooks=std::vector<std::string>(1, "SuppressSmallPT");
99 }
100
101 if (!canSetHook) {
102 ATH_MSG_ERROR(" *** Unable to initialise PythiaB !! ***");
103 return StatusCode::FAILURE;
104 }
105
106 if(m_userString == "NONE") m_userString.clear();
107 // Call the base class genInitialize()
109
110 ATH_MSG_DEBUG("... seeding Athena random number generator");
111 if (m_useRndmGenSvc){
112 p_rndmEngine = m_atlasRndmEngine->getEngine(); // NOT THREAD-SAFE
113 if (!p_rndmEngine) {
114 ATH_MSG_FATAL("Unable to retrieve HepRandomEngine. Bailing out.");
115 return StatusCode::FAILURE;
116 }
117 }
118 return StatusCode::SUCCESS;
119}
#define ATH_MSG_FATAL(x)
virtual StatusCode genInitialize()
For initializing the generator, if required.
Definition Pythia8_i.cxx:98
StringArrayProperty m_userHooks
Definition Pythia8_i.h:98

◆ genuserInitialize()

StatusCode Pythia8B_i::genuserInitialize ( )
virtual

For initialization of user code, if required. Called after genInitialize.

Reimplemented from GenModule.

Definition at line 123 of file Pythia8B_i.cxx.

123 {
124
125 // Just counter setting
126 ATH_MSG_INFO("Pythia8B_i from genuserInitialize()");
127
128 // Initialize global counters
129 m_totalBQuark = 0;
131 m_totalCQuark = 0;
134 m_totalHard = 0;
135 m_totalClone = 0;
136 m_atLeastOneAcc = 0;
139 m_speciesCount.clear();
140 for (std::vector<int>::iterator iit=m_bcodes.begin(); iit!=m_bcodes.end(); ++iit) {
141 m_speciesCount[*iit] = 0;
142 }
143
144 return StatusCode::SUCCESS;
145}

◆ getCodes()

std::vector< int > Pythia8B_i::getCodes ( const std::vector< Pythia8::Particle > & theParticles) const

Definition at line 617 of file Pythia8B_i.cxx.

617 {
618
619 std::vector<int> codes;
620 codes.reserve(theParticles.size());
621 for (auto pItr = theParticles.begin(); pItr!=theParticles.end(); ++pItr ) {
622 codes.push_back( (*pItr).id() );
623 }
624 return codes;
625
626}
dict codes
helper to get a human-readable string
Definition ExitCodes.py:49

◆ getRandomEngine() [1/2]

CLHEP::HepRandomEngine * GenModule::getRandomEngine ( const std::string & streamName,
const EventContext & ctx ) const
protectedinherited

Definition at line 34 of file GenModule.cxx.

36{
37 ATHRNG::RNGWrapper* rngWrapper = m_rndmSvc->getEngine(this, streamName);
38 std::string rngName = name()+streamName;
39 rngWrapper->setSeed( rngName, ctx );
40 return rngWrapper->getEngine(ctx);
41}
void setSeed(const std::string &algName, const EventContext &ctx)
Set the random seed using a string (e.g.
Definition RNGWrapper.h:154
CLHEP::HepRandomEngine * getEngine(const EventContext &ctx) const
Retrieve the random engine corresponding to the provided EventContext.
Definition RNGWrapper.h:108
ServiceHandle< IAthRNGSvc > m_rndmSvc
Data members.
Definition GenModule.h:97

◆ getRandomEngine() [2/2]

CLHEP::HepRandomEngine * GenModule::getRandomEngine ( const std::string & streamName,
unsigned long int randomSeedOffset,
const EventContext & ctx ) const
protectedinherited

Definition at line 44 of file GenModule.cxx.

46{
47 ATHRNG::RNGWrapper* rngWrapper = m_rndmSvc->getEngine(this, streamName);
48 rngWrapper->setSeed( streamName, ctx.slot(), randomSeedOffset, ctx.eventID().run_number() );
49 return rngWrapper->getEngine(ctx);
50}

◆ getRandomEngineDuringInitialize()

CLHEP::HepRandomEngine * GenModule::getRandomEngineDuringInitialize ( const std::string & streamName,
unsigned long int randomSeedOffset,
unsigned int conditionsRun = 1,
unsigned int lbn = 1 ) const
protectedinherited

Definition at line 53 of file GenModule.cxx.

54{
55 const size_t slot=0;
56 EventContext ctx;
57 ctx.setSlot( slot );
58 ctx.setEventID (EventIDBase (conditionsRun,
59 EventIDBase::UNDEFEVT, // event
60 EventIDBase::UNDEFNUM, // timestamp
61 EventIDBase::UNDEFNUM, // timestamp ns
62 lbn));
64 Atlas::ExtendedEventContext( evtStore()->hiveProxyDict(),
65 conditionsRun) );
66 return getRandomEngine(streamName, randomSeedOffset, ctx);
67}
ServiceHandle< StoreGateSvc > & evtStore()
CLHEP::HepRandomEngine * getRandomEngine(const std::string &streamName, const EventContext &ctx) const
Definition GenModule.cxx:34
void setExtendedEventContext(EventContext &ctx, ExtendedEventContext &&ectx)
Move an extended context into a context object.

◆ initialize()

StatusCode GenModule::initialize ( )
virtualinherited

Reimplemented from GenBase.

Definition at line 21 of file GenModule.cxx.

21 {
22 // Base class initializations
24 // Get the random number service
25 CHECK(m_rndmSvc.retrieve());
26 // Get the incident service
27 CHECK(m_incidentSvc.retrieve());
30 return StatusCode::SUCCESS;
31}
virtual StatusCode initialize() override
Definition GenBase.cxx:17
virtual StatusCode genuserInitialize()
For initialization of user code, if required. Called after genInitialize.
Definition GenModule.h:64
virtual StatusCode genInitialize()
For initializing the generator, if required.
Definition GenModule.h:62

◆ inputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::inputHandles ( ) const
overridevirtualinherited

Return this algorithm's input handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ isClonable()

virtual bool AthCommonAlgorithm< Gaudi::Algorithm >::isClonable ( ) const
inlineoverridevirtualinherited

Specify if the algorithm is clonable.

Only relevant for non-reentrant algorithms. Actual number of clones needs to be set via the "Cardinality" property.

Reimplemented in AFP_DigiTop, AlgB, AlgT, BCM_Digitization, CscDigitBuilder, CscDigitToCscRDO, G4AtlasAlg, G4RunAlg, HGTD_Digitization, HiveAlgBase, InDet::GNNSeedingTrackMaker, InDet::SCT_Clusterization, InDet::SiSPGNNTrackMaker, InDet::SiSPSeededTrackFinder, InDet::SiTrackerSpacePointFinder, ISF::SimKernelMT, ITk::StripDigitization, ITkPixelCablingAlg, ITkStripCablingAlg, LArHitEMapMaker, LArTTL1Maker, LUCID_DigiTop, LVL1::L1TopoSimulation, MergeCalibHits, MergeGenericMuonSimHitColl, MergeHijingPars, MergeMcEventCollection, MergeTrackRecordCollection, MergeTruthJets, MergeTruthParticles, MuonDigitizer, PileUpMTAlg, PixelDigitization, RoIBResultToxAOD, SCT_ByteStreamErrorsTestAlg, SCT_CablingCondAlgFromCoraCool, SCT_CablingCondAlgFromText, SCT_ConditionsParameterTestAlg, SCT_ConditionsSummaryTestAlg, SCT_ConfigurationConditionsTestAlg, SCT_Digitization, SCT_FlaggedConditionTestAlg, SCT_LinkMaskingTestAlg, SCT_MajorityConditionsTestAlg, SCT_ModuleVetoTestAlg, SCT_MonitorConditionsTestAlg, SCT_PrepDataToxAOD, SCT_RawDataToxAOD, SCT_ReadCalibChipDataTestAlg, SCT_ReadCalibDataTestAlg, SCT_RODVetoTestAlg, SCT_SensorsTestAlg, SCT_SiliconConditionsTestAlg, SCT_StripVetoTestAlg, SCT_TdaqEnabledTestAlg, SCT_TestCablingAlg, SCTEventFlagWriter, SCTRawDataProvider, SCTSiLorentzAngleTestAlg, SCTSiPropertiesTestAlg, SGInputLoader, Simulation::BeamEffectsAlg, TileHitVecToCnt, TileMuonFitter, TilePulseForTileMuonReceiver, TileRawChannelMaker, TRTDigitization, and ZDC_DigiTop.

Definition at line 68 of file AthCommonAlgorithm.h.

68 {
69 return true;
70 }

◆ isReEntrant()

virtual bool AthAlgorithm::isReEntrant ( ) const
inlinefinaloverrideprotectedvirtualinherited

Legacy algorithms are not thread-safe.

Definition at line 47 of file AthAlgorithm.h.

47{ return false; }

◆ leptonSelect()

bool Pythia8B_i::leptonSelect ( Pythia8::Event & theEvent,
const std::vector< double > & ptCut,
double etaCut,
const std::vector< int > & counts,
int type_id,
double massCut,
bool opposite )

Definition at line 493 of file Pythia8B_i.cxx.

493 {
494
495 if (type_id==0) return true;
496 bool passed(false);
497 std::string accString(" : REJECTED");
498
499 std::vector<int> leptonIDs;
500 int nCuts=ptCut.size();
501 std::vector<int> countGood(nCuts, 0);
502
503 for (int i = 0; i<theEvent.size(); ++i) {
504 const Pythia8::Particle &theParticle = theEvent[i];
505 int id = theParticle.idAbs();
506 if ( id == type_id ) { // find lepton flavour requested by user
507 double pt = theParticle.pT();
508 double eta = theParticle.eta();
509 ATH_MSG_DEBUG("Lepton of type " << id << " with pt/eta " << pt << "/" << eta);
510 for (int cutCntr=0; cutCntr<nCuts; ++cutCntr) {
511 if ( (pt>ptCut[cutCntr]) && (std::abs(eta)<etaCut) ) {
512 countGood[cutCntr] += 1;
513 leptonIDs.push_back(i); // remember leptons
514 }
515 }
516 }
517 // can the loop stop?
518 int nPassed(0);
519 for (int cutCntr=0; cutCntr<nCuts; ++cutCntr) {
520 if (countGood[cutCntr] >= counts[cutCntr]) ++nPassed;
521 }
522 if (nPassed==nCuts) break;
523 } // end of loop over particles
524
525 // Check the accumulated counts
526 int nPassed(0);
527 for (int cutCntr=0; cutCntr<nCuts; ++cutCntr) {
528 if (countGood[cutCntr] >= counts[cutCntr]) ++nPassed;
529 }
530 if (nPassed >= nCuts && nCuts==1) {accString=" : ACCEPTED"; passed = true;} // only 1 lepton required so no pair properties
531 if (nPassed >= nCuts && nCuts>1 && pairProperties(theEvent,leptonIDs,massCut,opposite)) {accString=" : ACCEPTED"; passed = true;}
532 ATH_MSG_DEBUG("Trigger-like selection" << accString);
533 return passed;
534}
bool passed(DecisionID id, const DecisionIDContainer &)
checks if required decision ID is in the set of IDs in the container
bool pairProperties(Pythia8::Event &, const std::vector< int > &, double, bool)

◆ msg()

MsgStream & AthCommonMsg< Gaudi::Algorithm >::msg ( ) const
inlineinherited

Definition at line 24 of file AthCommonMsg.h.

24 {
25 return this->msgStream();
26 }

◆ msgLvl()

bool AthCommonMsg< Gaudi::Algorithm >::msgLvl ( const MSG::Level lvl) const
inlineinherited

Definition at line 30 of file AthCommonMsg.h.

30 {
31 return this->msgLevel(lvl);
32 }

◆ outputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::outputHandles ( ) const
overridevirtualinherited

Return this algorithm's output handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ pairProperties()

bool Pythia8B_i::pairProperties ( Pythia8::Event & theEvent,
const std::vector< int > & leptonIDs,
double massCut,
bool opposite )

Definition at line 563 of file Pythia8B_i.cxx.

563 {
564
565 bool passesCuts(false);
566 std::string accString(" : REJECTED");
567 for (auto iit = leptonIDs.begin(); iit!=leptonIDs.end(); ++iit) {
568 for (auto iit2 = iit+1; iit2!=leptonIDs.end(); ++iit2) {
569 int q1=theEvent[*iit].charge();
570 int q2=theEvent[*iit2].charge();
571 if (opposite && (q1*q2>0)) continue;
572 double px1=theEvent[*iit].px();
573 double py1=theEvent[*iit].py();
574 double pz1=theEvent[*iit].pz();
575 double mass1=theEvent[*iit].mSel();
576 double e1=std::sqrt(px1*px1+py1*py1+pz1*pz1+mass1*mass1);
577 double px2=theEvent[*iit2].px();
578 double py2=theEvent[*iit2].py();
579 double pz2=theEvent[*iit2].pz();
580 double mass2=theEvent[*iit2].mSel();
581 double e2=std::sqrt(px2*px2+py2*py2+pz2*pz2+mass2*mass2);
582 double eSum=e1+e2;
583 double pxSum=px1+px2;
584 double pySum=py1+py2;
585 double pzSum=pz1+pz2;
586 double M=std::sqrt(eSum*eSum-pxSum*pxSum-pySum*pySum-pzSum*pzSum);
587 if (M>massCut) {
588 passesCuts=true;
589 ATH_MSG_DEBUG("Acceptable lepton pair with invariant mass : " << M);
590 break;
591 }
592 }
593 if (passesCuts) {accString=" : ACCEPTED"; break;}
594 }
595 ATH_MSG_DEBUG("Dilepton selection" << accString);
596 return passesCuts;
597}
double e2(const xAOD::CaloCluster &cluster)
return the uncorrected cluster energy in 2nd sampling
double e1(const xAOD::CaloCluster &cluster)
return the uncorrected cluster energy in 1st sampling

◆ passesEtaCuts()

bool Pythia8B_i::passesEtaCuts ( const std::vector< Pythia8::Particle > & theParticles) const

Definition at line 662 of file Pythia8B_i.cxx.

662 {
663
664 bool pass(true);
665 unsigned int i(0);
666 for (auto pItr=theParticles.cbegin(); pItr!=theParticles.cend(); ++pItr,++i) {
667 if (std::abs((*pItr).eta()) > m_sigEtaCuts[i]) pass = false;
668 if (!pass) break;
669 }
670
671 return pass;
672
673}

◆ passesPTCuts()

bool Pythia8B_i::passesPTCuts ( const std::vector< Pythia8::Particle > & theParticles) const

Definition at line 650 of file Pythia8B_i.cxx.

650 {
651
652 bool pass(true);
653 unsigned int i(0);
654 for (auto pItr=theParticles.cbegin(); pItr!=theParticles.cend(); ++pItr,++i) {
655 if ((*pItr).pT() < m_sigPtCuts[i]) pass = false;
656 if (!pass) break;
657 }
658 return pass;
659
660}

◆ printSignalSelections()

void Pythia8B_i::printSignalSelections ( const std::vector< int > & signalProcess,
const std::vector< double > & ptCuts,
const std::vector< double > & etaCuts,
unsigned int nRequired ) const

Definition at line 745 of file Pythia8B_i.cxx.

745 {
746 std::cout << "Signal PDG codes required: ";
747 for (unsigned int k=0; k<m_sigCodes.size(); ++k) std::cout << signalProcess[k] << " ";
748 if (signalProcess.size()==ptCuts.size()) {
749 std::cout << "Cuts on the pt of the signal states: " << std::endl;
750 for (unsigned int l=0; l<ptCuts.size(); ++l) std::cout << ptCuts[l] << " ";
751 }
752 if (signalProcess.size()==etaCuts.size()) {
753 std::cout << "Cuts on the eta of the signal states: " << std::endl;
754 for (unsigned int l=0; l<etaCuts.size(); ++l) std::cout << etaCuts[l] << " ";
755 }
756 std::cout << "Number of such decays required per event: " << nRequired << std::endl;
757 std::cout << std::endl;
758}
l
Printing final latex table to .tex output file.

◆ pythia_stream()

const std::string & Pythia8_i::pythia_stream ( )
staticinherited

Definition at line 688 of file Pythia8_i.cxx.

688 {
689 return s_pythia_stream;
690}

◆ pythiaVersion()

double Pythia8_i::pythiaVersion ( ) const
inherited

Definition at line 683 of file Pythia8_i.cxx.

683 {
684 return m_version;
685}
double m_version
Definition Pythia8_i.h:110

◆ renounce()

std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::renounce ( T & h)
inlineprotectedinherited

Definition at line 380 of file AthCommonDataStore.h.

381 {
382 h.renounce();
384 }
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce(T &h)

◆ renounceArray()

void AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::renounceArray ( SG::VarHandleKeyArray & handlesArray)
inlineprotectedinherited

remove all handles from I/O resolution

Definition at line 364 of file AthCommonDataStore.h.

364 {
366 }

◆ s_allowedTunes()

int Pythia8_i::s_allowedTunes ( double version)
staticprivateinherited

◆ setFilterPassed()

virtual void AthCommonAlgorithm< Gaudi::Algorithm >::setFilterPassed ( bool state,
const EventContext & ctx ) const
inlinevirtualinherited

Set filter decision:

Reimplemented in AthFilterAlgorithm.

Definition at line 99 of file AthCommonAlgorithm.h.

99 {
101 }
virtual void setFilterPassed(bool state, const EventContext &ctx) const
Set filter decision:

◆ signalAccept()

bool Pythia8B_i::signalAccept ( Pythia8::Event & theEvent,
const std::vector< int > & requiredDecay,
unsigned int nRequired ) const

Definition at line 678 of file Pythia8B_i.cxx.

679 {
680
681 bool acceptEvent(false);
682 std::vector<int> parentsIndices;
683 for (int i = 0; i<theEvent.size(); ++i) {
684 const Pythia8::Particle &theParticle = theEvent[i];
685 int id = theParticle.id();
686 if (id==requiredDecay[0]) parentsIndices.push_back(i);
687 }
688
689 // For the special case where the user merely requires a particular
690 // species to be in the event (mainly for EvtGen use)
691 if ( (requiredDecay.size()==1) && (parentsIndices.size()>0) ) {
692 acceptEvent = true;
693 return acceptEvent;
694 }
695
696 unsigned int goodDecayCounter(0);
697 for (std::vector<int>::iterator iItr = parentsIndices.begin(); iItr!=parentsIndices.end(); ++iItr) {
698 std::vector<Pythia8::Particle> decayMembers;
699 descendThroughDecay(theEvent,decayMembers,*iItr);
700 std::vector<int> pdgCodes = getCodes(decayMembers);
701 if (!compare(requiredDecay,std::move(pdgCodes))) {
702 ATH_MSG_DEBUG("Signal event REJECTED as does not contain required decay chain");
703 continue;
704 }
705 ATH_MSG_DEBUG("Event contains required signal decay chain");
706
707 if (decayMembers.size()==m_sigPtCuts.size()) {
708 if (!passesPTCuts(decayMembers)) {
709 ATH_MSG_DEBUG("Signal event REJECTED as signal chain does not pass pt cuts");
710 continue;
711 }
712 }
713
714 if (decayMembers.size()==m_sigEtaCuts.size()) {
715 if (!passesEtaCuts(decayMembers)) {
716 ATH_MSG_DEBUG("Signal event REJECTED as signal chain does not pass eta cuts");
717 continue;
718 }
719 }
720 if (nRequired==1) {
721 ATH_MSG_DEBUG("Signal decay good; event ACCEPTED");
722 acceptEvent = true;
723 break;
724 }
725 else {++goodDecayCounter;}
726 }
727 if (nRequired>1) {
728 if (goodDecayCounter>=nRequired) {
729 ATH_MSG_DEBUG(nRequired << "signal decays required; " << goodDecayCounter << " found; event ACCEPTED");
730 acceptEvent = true;
731 }
732 else if (goodDecayCounter<nRequired) {
733 ATH_MSG_DEBUG(nRequired << "signal decays required; " << goodDecayCounter << " found; event REJECTED");
734 acceptEvent = false;
735 }
736 }
737
738 return acceptEvent;
739
740}
bool compare(std::vector< int >, std::vector< int >) const
bool passesPTCuts(const std::vector< Pythia8::Particle > &) const
bool passesEtaCuts(const std::vector< Pythia8::Particle > &) const
std::vector< int > getCodes(const std::vector< Pythia8::Particle > &) const

◆ sysExecute()

StatusCode AthCommonAlgorithm< Gaudi::Algorithm >::sysExecute ( const EventContext & ctx)
overridevirtualinherited

Execute an algorithm.

We override this in order to work around an issue with the Algorithm base class storing the event context in a member variable that can cause crashes in MT jobs.

Reimplemented in AthAnalysisAlgorithm.

Definition at line 80 of file AthCommonAlgorithm.cxx.

41{
42 return BaseAlg::sysExecute (ctx);
43}

◆ sysInitialize()

StatusCode AthCommonAlgorithm< Gaudi::Algorithm >::sysInitialize ( )
overridevirtualinherited

Override sysInitialize.

Override sysInitialize from the base class.

Loop through all output handles, and if they're WriteCondHandles, automatically register them and this Algorithm with the CondSvc

Scan through all outputHandles, and if they're WriteCondHandles, register them with the CondSvc

Reimplemented from AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >.

Reimplemented in AthAnalysisAlgorithm, AthFilterAlgorithm, AthHistogramAlgorithm, HypoBase, InputMakerBase, and PyAthena::Alg.

Definition at line 60 of file AthCommonAlgorithm.cxx.

71 {
73
74 if (sc.isFailure()) {
75 return sc;
76 }
77
78 ServiceHandle<ICondSvc> cs("CondSvc",name());
79 for (auto h : outputHandles()) {
80 if (h->isCondition() && h->mode() == Gaudi::DataHandle::Writer) {
81 // do this inside the loop so we don't create the CondSvc until needed
82 if ( cs.retrieve().isFailure() ) {
83 ATH_MSG_WARNING("no CondSvc found: won't autoreg WriteCondHandles");
85 }
86 if (cs->regHandle(this,*h).isFailure()) {
88 ATH_MSG_ERROR("unable to register WriteCondHandle " << h->fullKey()
89 << " with CondSvc");
90 }
91 }
92 }
93 return sc;
94}
#define ATH_MSG_WARNING(x)
virtual StatusCode sysInitialize() override
virtual std::vector< Gaudi::DataHandle * > outputHandles() const override

◆ sysStart()

virtual StatusCode AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::sysStart ( )
overridevirtualinherited

Handle START transition.

We override this in order to make sure that conditions handle keys can cache a pointer to the conditions container.

◆ updateVHKA()

void AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::updateVHKA ( Gaudi::Details::PropertyBase & )
inlineinherited

Definition at line 308 of file AthCommonDataStore.h.

308 {
309 // debug() << "updateVHKA for property " << p.name() << " " << p.toString()
310 // << " size: " << m_vhka.size() << endmsg;
311 for (auto &a : m_vhka) {
313 for (auto k : keys) {
314 k->setOwner(this);
315 }
316 }
317 }

◆ useReseed()

bool Pythia8_i::useReseed ( ) const
inlineprotectedinherited

Definition at line 86 of file Pythia8_i.h.

86{return m_useReseed; }
BooleanProperty m_useReseed
Definition Pythia8_i.h:95

◆ useRndmGenSvc()

bool Pythia8_i::useRndmGenSvc ( ) const
inlineprotectedinherited

Definition at line 85 of file Pythia8_i.h.

85{ return m_useRndmGenSvc; }

◆ userSelection()

bool Pythia8B_i::userSelection ( Pythia8::Event & event,
const std::string & userString,
const std::vector< double > & userVars )

Definition at line 237 of file UserSelections.h.

238 {
239
240 using CLHEP::HepLorentzVector;
241 bool accept(false);
242
243 // ==================================================
244 // code up your selection here; as many blocks as
245 // needed can be inserted
246 if (userString == "EXAMPLE" && event.size() > 0 && userVars.size() > 0) {
247
248 accept = true;
249
250 }
251
252 // ==================================================
253 // selects decays where a B-hadron decays to a J/psi
254 // by any route
255 // Accepts event if event contains a B which (eventually)
256 // goes to a J/psi
257
258 else if (userString == "BJPSIINCLUSIVE" && event.size() > 0) {
259 int chargeConj = 1;
260 if (userVars.size()==0) {
261 ATH_MSG_INFO("User selection BJPSIINCLUSIVE with B-state");
262 }
263 if (userVars.size()>0) {
264 if (userVars.size()>1) ATH_MSG_WARNING("User selection BJPSIINCLUSIVE with more than one argument! Check job options");
265 if (userVars.at(0)>0) ATH_MSG_DEBUG("User selection BJPSIINCLUSIVE with B-state");
266 if (userVars.at(0)<0) {ATH_MSG_DEBUG("User selection BJPSIINCLUSIVE with anti-B-state"); chargeConj = -1;}
267 }
268
269 // B-decay codes which can go to charmonium
270 std::vector<int> bToCharmoniaCodes;
271 bToCharmoniaCodes.push_back(511*chargeConj);
272 bToCharmoniaCodes.push_back(521*chargeConj);
273 bToCharmoniaCodes.push_back(531*chargeConj);
274 bToCharmoniaCodes.push_back(541*chargeConj);
275 bToCharmoniaCodes.push_back(-5122*chargeConj);
276 bToCharmoniaCodes.push_back(-5132*chargeConj);
277 bToCharmoniaCodes.push_back(-5232*chargeConj);
278 bToCharmoniaCodes.push_back(-5332*chargeConj);
279
280 int eventSize = event.size();
281
282 bool isBtoJpsi(false);
283 bool containsB(false);
284 // loop over all particles in the event
285 for (int i = 0; i < eventSize; i++) {
286 int pdgID = event[i].id();
287 std::vector < Pythia8::Particle > decayMembers;
288 bool isB(false);
289 // Does event contain B which can go to J/psi?
290 for (unsigned int k = 0; k < bToCharmoniaCodes.size(); ++k) {
291 if (pdgID == bToCharmoniaCodes[k]) {
292 containsB = true;
293 isB = true;
294 break;
295 }
296 }
297 // Get decay chain of the B; see if there is a J/psi
298 if (isB) {
299 descendThroughDecay(event, decayMembers, i);
300 const std::vector<int> & pdgCodes = getCodes(decayMembers);
301 for (unsigned int k = 0; k < pdgCodes.size(); ++k) {
302 if (pdgCodes[k] == 443)
303 isBtoJpsi = true;
304 }
305 }
306 if (isBtoJpsi)
307 break;
308 }
309 if (containsB && isBtoJpsi)
310 accept = true;
311 if (containsB && !isBtoJpsi)
312 accept = false;
313 if (!containsB)
314 accept = false;
315
316 }
317
318 // ==================================================
319 // prototype for Bs->J/psi phi angular distribution
320 // (Adam Barton)
321
322 else if ((userString == "BJPSIPHI_TRANS" || userString == "BJPSIPHI_HEL")
323 && event.size() > 0) {
324
325 const bool debug = false;
326 bool flat;
327
328 //Read
329 //userVars[0] 0(flat)/1(angle)
330 //[1] prob_limit
331 //[2] tag mode
332 //[3] A0
333 //[4] Al
334 //[5] As
335 //[6] GammaS
336 //[7] DeltaGamma
337 //[8] DeltaM
338 //[9] phiS
339 //[10] delta_p
340 //[11] delta_l
341 //[12]delta_s
342
343 if (userVars.size() < 13) {
345 "REQUIRED userVARs not provided for BsJpsiPhi pdf defaulting to flat");
346 flat = true;
347 } else {
348 flat = userVars[0] == 0.0;
349
350 }
351
352 // for(int i=0; i<10;i++){
353 // ATH_MSG_INFO("BJPSIPHI_TRANS: " << i << " random number: " << Rdmengine->flat());
354 //
355 // }
356
357 // std::vector<Pythia8::Particle> decayMembers;
358
359 int i_Bs = 0;
360 int i_Muplus = 0;
361 int i_Muminus = 0;
362 int i_Kplus = 0;
363 int i_Kminus = 0;
364 int i_Phi = 0;
365 int i_Jpsi = 0;
366
367 bool isBsJpsiPhi = false;
368 int eventSize = event.size();
369 for (int i = 0; i < eventSize; i++) {
370
371 int pID = event[i].id();
372 if (std::abs(pID) == 531) { //NOTE THIS WILL FIND BS AND ANTIBS
373 i_Bs = i;
374 std::vector<int> daughterlist = event.daughterList(i);
375
376 if (daughterlist.size() != 2)
377 continue;
378 bool isjpsi = false;
379 bool isphi = false;
380 if (event[daughterlist[0]].id() == 443) {
381 isjpsi = true;
382 i_Jpsi = daughterlist[0];
383 }
384 if (event[daughterlist[1]].id() == 443) {
385 isjpsi = true;
386 i_Jpsi = daughterlist[1];
387 }
388 if (event[daughterlist[0]].id() == 333) {
389 isphi = true;
390 i_Phi = daughterlist[0];
391 }
392 if (event[daughterlist[1]].id() == 333) {
393 isphi = true;
394 i_Phi = daughterlist[1];
395 }
396 if (!isphi || !isjpsi)
397 continue;
398 std::vector<int> daughterlistJpsi = event.daughterList(i_Jpsi);
399 std::vector<int> daughterlistPhi = event.daughterList(i_Phi);
400 if (daughterlistJpsi.size() != 2 || daughterlistPhi.size() != 2)
401 continue;
402 //resets values to avoid possible bug
403
404 if (event[daughterlistJpsi[0]].id() == 13)
405 i_Muminus = daughterlistJpsi[0];
406 else if (event[daughterlistJpsi[1]].id() == 13)
407 i_Muminus = daughterlistJpsi[1];
408 else
409 i_Muminus = 0;
410
411 if (event[daughterlistJpsi[0]].id() == -13)
412 i_Muplus = daughterlistJpsi[0];
413 else if (event[daughterlistJpsi[1]].id() == -13)
414 i_Muplus = daughterlistJpsi[1];
415 else
416 i_Muplus = 0;
417
418 if (event[daughterlistPhi[0]].id() == 321)
419 i_Kplus = daughterlistPhi[0];
420 else if (event[daughterlistPhi[1]].id() == 321)
421 i_Kplus = daughterlistPhi[1];
422 else
423 i_Kplus = 0;
424
425 if (event[daughterlistPhi[0]].id() == -321)
426 i_Kminus = daughterlistPhi[0];
427 else if (event[daughterlistPhi[1]].id() == -321)
428 i_Kminus = daughterlistPhi[1];
429 else
430 i_Kminus = 0;
431 if (i_Muminus == 0 || i_Muplus == 0 || i_Kminus == 0 || i_Kplus
432 == 0)
433 continue;
434 else {
435 if (debug)
437 "found entire Bs->Jpsi(mumu)phi(KK) decay ");
438 isBsJpsiPhi = true;
439 break;
440 }
441 }
442 }
443
444 if (!isBsJpsiPhi)
445 return false;
446 if (flat)
447 return true;
448
449 Pythia8::Particle &p_Bs = event[i_Bs];
450 Pythia8::Particle &p_Muplus = event[i_Muplus];
451 Pythia8::Particle &p_Muminus = event[i_Muminus];
452 Pythia8::Particle &p_Kplus = event[i_Kplus];
453 Pythia8::Particle &p_Kminus = event[i_Kminus];
454 Pythia8::Particle &p_Phi = event[i_Phi];
455 Pythia8::Particle &p_Jpsi = event[i_Jpsi];
456
457 // cout << "Bs " << p_Bs.id() << endl;
458 // cout << "|> " << p_Jpsi.id() << " + " << p_Phi.id() << endl;
459 // cout << "|> " << p_Muplus.id() << " + " << p_Muminus.id() << " + " << p_Kplus.id() << " + " << p_Kminus.id() << endl;
460
461
462 HepLorentzVector v_Bs = convertVector(p_Bs.p());
463 HepLorentzVector v_Muplus = convertVector(p_Muplus.p());
464 HepLorentzVector v_Muminus = convertVector(p_Muminus.p());
465 HepLorentzVector v_Kplus = convertVector(p_Kplus.p());
466 HepLorentzVector v_Kminus = convertVector(p_Kminus.p());
467 HepLorentzVector v_Phi = convertVector(p_Phi.p());
468 HepLorentzVector v_Jpsi = convertVector(p_Jpsi.p());
469
470 BsJpsiPhiAngles angles(v_Kplus, v_Muplus, v_Phi, v_Jpsi, v_Bs);
471
472 CLHEP::HepRandomEngine* Rdmengine = Pythia8B_i::p_rndmEngine;
473 const double gentau = Pythia8_i::m_pythia->particleData.tau0(531);
474 const double correctionfactor = 0.299792458;
475 const double gentauCorrect = gentau / correctionfactor;
476 if (debug) {
477 ATH_MSG_INFO("Lifetime for 531: " << gentau);
478 ATH_MSG_INFO("correct lifetime " << gentauCorrect);
479 }
480 const double Bstau = p_Bs.tau() / correctionfactor;
481 double prob1 = 1000;
482
483 //PUT PDFS HERE
484 if (userString == "BJPSIPHI_TRANS") {
485 double x[5];
486 x[0] = Bstau;
487 x[1] = angles.thetatrfix();
488 x[2] = angles.thetaKfix();
489 x[3] = angles.phitrfix();
490 x[4] = userVars[2];
491
492 prob1 = BsJpsiPhi_PDF(&userVars[3], x, false);
493 } else if (userString == "BJPSIPHI_HEL") {
494 double x[5];
495 x[0] = Bstau;
496 x[1] = angles.thetaLfix();
497 x[2] = angles.thetaKfix();
498 x[3] = angles.chifix();
499 x[4] = userVars[2];
500
501 prob1 = BsJpsiPhi_PDF(&userVars[3], x, true);
502 }
503
504 const double prob2 = exp(Bstau / gentauCorrect) * gentauCorrect;
505 if (Bstau > 20)
506 ATH_MSG_WARNING("Warning Bstau > 20 ps, this could indicate a bug");
507 const double prob_norm = userVars[1];
508 if (debug) {
509 ATH_MSG_INFO("prob limit set to " << prob_norm);
510 ATH_MSG_INFO("This Bs has lifetime " << Bstau);
511 }
512 double rand = Rdmengine->flat() * prob_norm;
513 if (prob1 * prob2 > prob_norm) {
515 "Max prob exceeded, too many of these indicates a problem, a few is fine");
516 }
517 if (debug)
518 std::cout << "totalprob " << prob1 * prob2 << std::endl;
519 if (rand < (prob1 * prob2)) {
520
521 if (debug)
522 ATH_MSG_INFO("Passed PDF filter");
523 accept = true;
524 } else {
525 if (debug)
526 ATH_MSG_INFO("Rejected PDF filter");
527 accept = false;
528 }
529
530 } // end of Bs->J/psiphi
531 else if(userString == "BD_BPLUS_TAUCUT"){
532 const double correctionfactor = 0.299792458;
533 bool pass = true;
534 int eventSize = event.size();
535 int foundcount=0;
536 for (int i = 0; i < eventSize; i++) {
537 int pID = event[i].id();
538 if (pID == 511) {
539 const double Bdtau = event[i].tau() / correctionfactor;
540 pass &= Bdtau > userVars[0] && Bdtau < userVars[1];
541 foundcount++;
542 }
543 if (pID == 521) {
544 const double Bptau = event[i].tau() / correctionfactor;
545 pass &= Bptau > userVars[0] && Bptau < userVars[1];
546 foundcount++;
547 }
548 }
549 accept = (foundcount > 0) & pass;
550 }
551
552 else if ((userString == "BDJPSIKSTAR_TRANS") && event.size() > 0) {
553 const bool debug = false;
554 bool flat;
555 if (userVars.size() < 13) {
557 "REQUIRED userVARs not provided for BdJpsiKstar pdf defaulting to flat");
558 flat = true;
559 } else {
560 flat = userVars[0] == 0.0;
561 }
562
563 int i_Bd = 0;
564 int i_Muplus = 0;
565 int i_Muminus = 0;
566 int i_Kplus = 0;
567 int i_piminus = 0;
568 int i_Kstar = 0;
569 int i_Jpsi = 0;
570
571 bool isBsJpsiKstar = false;
572 const int eventSize = event.size();
573 for (int i = 0; i < eventSize; i++) {
574
575 const int pID = event[i].id();
576 if (std::abs(pID) == 511) { //NOTE THIS FIND BD and Anti-Bd
577 i_Bd = i;
578 std::vector<int> daughterlist = event.daughterList(i);
579
580 if (daughterlist.size() != 2)
581 continue;
582 bool isjpsi = false;
583 bool iskstar = false;
584 if (event[daughterlist[0]].id() == 443) {
585 isjpsi = true;
586 i_Jpsi = daughterlist[0];
587 }
588 if (event[daughterlist[1]].id() == 443) {
589 isjpsi = true;
590 i_Jpsi = daughterlist[1];
591 }
592 if (std::abs(event[daughterlist[0]].id()) == 313) { //This will find kstar or KstarBar
593 iskstar = true;
594 i_Kstar = daughterlist[0];
595 }
596 if (std::abs(event[daughterlist[1]].id()) == 313) { //This will find kstar or KstarBar
597 iskstar = true;
598 i_Kstar = daughterlist[1];
599 }
600 if (!iskstar || !isjpsi)
601 continue;
602 std::vector<int> daughterlistJpsi = event.daughterList(i_Jpsi);
603 std::vector<int> daughterlistKstar = event.daughterList(i_Kstar);
604 if (daughterlistJpsi.size() != 2 || daughterlistKstar.size() != 2)
605 continue;
606 //resets values to avoid possible bug
607
608 if (event[daughterlistJpsi[0]].id() == 13)
609 i_Muminus = daughterlistJpsi[0];
610 else if (event[daughterlistJpsi[1]].id() == 13)
611 i_Muminus = daughterlistJpsi[1];
612 else
613 i_Muminus = 0;
614
615 if (event[daughterlistJpsi[0]].id() == -13)
616 i_Muplus = daughterlistJpsi[0];
617 else if (event[daughterlistJpsi[1]].id() == -13)
618 i_Muplus = daughterlistJpsi[1];
619 else
620 i_Muplus = 0;
621
622 if (std::abs(event[daughterlistKstar[0]].id()) == 321)
623 i_Kplus = daughterlistKstar[0];
624 else if (std::abs(event[daughterlistKstar[1]].id()) == 321)
625 i_Kplus = daughterlistKstar[1];
626 else
627 i_Kplus = 0;
628
629 if (std::abs(event[daughterlistKstar[0]].id()) == 211)
630 i_piminus = daughterlistKstar[0];
631 else if (std::abs(event[daughterlistKstar[1]].id()) == 211)
632 i_piminus = daughterlistKstar[1];
633 else
634 i_piminus = 0;
635
636 if (i_Muminus == 0 || i_Muplus == 0 || i_piminus == 0 || i_Kplus
637 == 0)
638 continue;
639
640 if (debug)
642 "found entire Bd->Jpsi(mumu)Kstar(KPi) decay ");
643 isBsJpsiKstar = true;
644 break;
645
646 }
647 }
648
649 if (!isBsJpsiKstar)
650 return false;
651 if (flat)
652 return true;
653
654 Pythia8::Particle &p_Bd = event[i_Bd];
655 Pythia8::Particle &p_Muplus = event[i_Muplus];
656 Pythia8::Particle &p_Muminus = event[i_Muminus];
657 Pythia8::Particle &p_Kplus = event[i_Kplus];
658 //Pythia8::Particle &p_piminus = event[i_piminus];
659 Pythia8::Particle &p_Kstar = event[i_Kstar];
660 Pythia8::Particle &p_Jpsi = event[i_Jpsi];
661
662 // cout << "Bs " << p_Bs.id() << endl;
663 // cout << "|> " << p_Jpsi.id() << " + " << p_Phi.id() << endl;
664 // cout << "|> " << p_Muplus.id() << " + " << p_Muminus.id() << " + " << p_Kplus.id() << " + " << p_Kminus.id() << endl;
665
666
667 HepLorentzVector v_Bd = convertVector(p_Bd.p());
668 HepLorentzVector v_Muplus = convertVector(p_Muplus.p());
669 HepLorentzVector v_Muminus = convertVector(p_Muminus.p());
670 HepLorentzVector v_Kplus = convertVector(p_Kplus.p());
671 //HepLorentzVector v_piminus = convertVector(p_piminus.p());
672 HepLorentzVector v_Kstar = convertVector(p_Kstar.p());
673 HepLorentzVector v_Jpsi = convertVector(p_Jpsi.p());
674
675 BsJpsiPhiAngles angles(v_Kplus, v_Muplus, v_Kstar, v_Jpsi, v_Bd);
676
677 CLHEP::HepRandomEngine* Rdmengine = Pythia8B_i::p_rndmEngine;
678 const double gentau = Pythia8_i::m_pythia->particleData.tau0(511);
679 const double correctionfactor = 0.299792458;
680 const double gentauCorrect = gentau / correctionfactor;
681 if (debug) {
682 ATH_MSG_INFO("Lifetime for 511: " << gentau);
683 ATH_MSG_INFO("correct lifetime " << gentauCorrect);
684 }
685 const double Bdtau = p_Bd.tau() / correctionfactor;
686 double prob1;
687
688 //PUT PDFS HERE
689
690 double x[5];
691 x[0] = Bdtau;
692 x[1] = angles.thetatrfix();
693 x[2] = angles.thetaKfix();
694 x[3] = angles.phitrfix();
695 x[4] = userVars[2];
696 prob1 = BsJpsiPhi_PDF(&userVars[3], x, false);
697
698
699 const double prob2 = exp(Bdtau / gentauCorrect) * gentauCorrect;
700 if (Bdtau > 20)
701 ATH_MSG_WARNING("Warning Bdtau > 20 ps, this could indicate a bug");
702 const double prob_norm = userVars[1];
703 if (debug) {
704 ATH_MSG_INFO("prob limit set to " << prob_norm);
705 ATH_MSG_INFO("This Bd has lifetime " << Bdtau);
706 }
707 const double rand = Rdmengine->flat() * prob_norm;
708 if (prob1 * prob2 > prob_norm) {
710 "Max prob exceeded, too many of these indicates a problem, a few is fine");
711 }
712 if (debug)
713 std::cout << "totalprob " << prob1 * prob2 << std::endl;
714 if (rand < (prob1 * prob2)) {
715 if (debug)
716 ATH_MSG_INFO("Passed PDF filter");
717 accept = true;
718 } else {
719 if (debug)
720 ATH_MSG_INFO("Rejected PDF filter");
721 accept = false;
722 }
723 }
724
725 return (accept);
726
727}
const bool debug
double BsJpsiPhi_PDF(const double *params, double *x, bool useHelicity)
CLHEP::HepLorentzVector convertVector(const Pythia8::Vec4 v)
#define x
StatusCode accept(const xAOD::Muon *mu)

◆ xmlpath()

std::string Pythia8_i::xmlpath ( )
staticinherited

Definition at line 693 of file Pythia8_i.cxx.

693 {
694
695
696 std::string foundpath = "";
697
698// Try to find the xmldoc directory using PathResolver:
699 foundpath = PathResolverFindCalibDirectory( "Pythia8/xmldoc" );
700
701 return foundpath;
702}
std::string PathResolverFindCalibDirectory(const std::string &logical_file_name)

Member Data Documentation

◆ m_and

bool Pythia8B_i::m_and
private

Definition at line 54 of file Pythia8B_i.h.

◆ m_aqEtaCut

double Pythia8B_i::m_aqEtaCut
private

Definition at line 55 of file Pythia8B_i.h.

◆ m_aqPtCut

double Pythia8B_i::m_aqPtCut
private

Definition at line 55 of file Pythia8B_i.h.

◆ m_athenaTool

PublicToolHandle<IPythia8Custom> Pythia8_i::m_athenaTool {this, "CustomInterface", ""}
privateinherited

Definition at line 174 of file Pythia8_i.h.

174{this, "CustomInterface", ""};

◆ m_atlasRndmEngine

std::shared_ptr<customRndm> Pythia8_i::m_atlasRndmEngine {}
protectedinherited

Definition at line 93 of file Pythia8_i.h.

93{};

◆ m_atLeastOneAcc

int Pythia8B_i::m_atLeastOneAcc
private

Definition at line 51 of file Pythia8B_i.h.

◆ m_bcodes

std::vector<int> Pythia8B_i::m_bcodes
private

Definition at line 48 of file Pythia8B_i.h.

◆ m_beam1

StringProperty Pythia8_i::m_beam1 {this, "Beam1", "PROTON"}
privateinherited

Definition at line 129 of file Pythia8_i.h.

129{this, "Beam1", "PROTON"};

◆ m_beam2

StringProperty Pythia8_i::m_beam2 {this, "Beam2", "PROTON"}
privateinherited

Definition at line 130 of file Pythia8_i.h.

130{this, "Beam2", "PROTON"};

◆ m_BEventBuffer

std::vector<Pythia8::Event> Pythia8B_i::m_BEventBuffer
private

Definition at line 58 of file Pythia8B_i.h.

◆ m_collisionEnergy

DoubleProperty Pythia8_i::m_collisionEnergy {this, "CollisionEnergy", 14000.0}
privateinherited

Definition at line 126 of file Pythia8_i.h.

126{this, "CollisionEnergy", 14000.0};

◆ m_commands

StringArrayProperty Pythia8_i::m_commands
privateinherited
Initial value:
{this,
"Commands",
{},
"List of Pythia8 commands",
"GeneratorSettings<std::string>"}

Definition at line 116 of file Pythia8_i.h.

116 {this,
117 "Commands",
118 {},
119 "List of Pythia8 commands",
120 "GeneratorSettings<std::string>"};

◆ m_computeEfficiency

BooleanProperty Pythia8_i::m_computeEfficiency {this, "computeEfficiency", false}
privateinherited

Definition at line 137 of file Pythia8_i.h.

137{this, "computeEfficiency", false};

◆ m_conversion

double Pythia8_i::m_conversion {1.}
privateinherited

Definition at line 141 of file Pythia8_i.h.

141{1.};

◆ m_cutCount

std::vector<int> Pythia8B_i::m_cutCount
private

Definition at line 48 of file Pythia8B_i.h.

◆ m_dec

unsigned int Pythia8B_i::m_dec
private

Definition at line 45 of file Pythia8B_i.h.

◆ m_detStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::m_detStore
privateinherited

Pointer to StoreGate (detector store by default).

Definition at line 393 of file AthCommonDataStore.h.

◆ m_doCKKWLAcceptance

BooleanProperty Pythia8_i::m_doCKKWLAcceptance {this, "CKKWLAcceptance", false}
privateinherited

Definition at line 135 of file Pythia8_i.h.

135{this, "CKKWLAcceptance", false};

◆ m_doFxFxXS

BooleanProperty Pythia8_i::m_doFxFxXS {this, "FxFxXS", false}
privateinherited

Definition at line 136 of file Pythia8_i.h.

136{this, "FxFxXS", false};

◆ m_doLHE3Weights

bool Pythia8_i::m_doLHE3Weights {false}
privateinherited

Definition at line 169 of file Pythia8_i.h.

169{false};

◆ m_doSuppressSmallPT

bool Pythia8B_i::m_doSuppressSmallPT
private

Definition at line 60 of file Pythia8B_i.h.

◆ m_dsid

IntegerProperty Pythia8_i::m_dsid {this, "Dsid", 999999, "Dataset ID number"}
protectedinherited

Definition at line 97 of file Pythia8_i.h.

97{this, "Dsid", 999999, "Dataset ID number"};

◆ m_enhanceWeight

double Pythia8_i::m_enhanceWeight {1.0}
privateinherited

Definition at line 166 of file Pythia8_i.h.

166{1.0};

◆ m_evtStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::m_evtStore
privateinherited

Pointer to StoreGate (event store by default).

Definition at line 390 of file AthCommonDataStore.h.

◆ m_extendedExtraObjects

DataObjIDColl AthCommonAlgorithm< Gaudi::Algorithm >::m_extendedExtraObjects
privateinherited

Extra output dependency collection, extended by AthAlgorithmDHUpdate to add symlinks.

Empty if no symlinks were found.

Definition at line 108 of file AthCommonAlgorithm.h.

◆ m_failureCount

unsigned int Pythia8B_i::m_failureCount
private

Definition at line 62 of file Pythia8B_i.h.

◆ m_had

unsigned int Pythia8B_i::m_had
private

Definition at line 45 of file Pythia8B_i.h.

◆ m_incidentSvc

ServiceHandle<IIncidentSvc> GenModule::m_incidentSvc {this, "IncidentSvc", "IncidentSvc"}
privateinherited

Handle on the incident service.

Definition at line 99 of file GenModule.h.

99{this, "IncidentSvc", "IncidentSvc"};

◆ m_internal_event_number

int Pythia8B_i::m_internal_event_number
private

Definition at line 51 of file Pythia8B_i.h.

◆ m_internalEventNumbers

std::vector<int> Pythia8B_i::m_internalEventNumbers
private

Definition at line 59 of file Pythia8B_i.h.

◆ m_invMass

double Pythia8B_i::m_invMass
private

Definition at line 55 of file Pythia8B_i.h.

◆ m_isAfterburner

BooleanProperty GenModule::m_isAfterburner {this, "IsAfterburner", false, "Set true if generator modifies existing events rather than creating new ones"}
protectedinherited

Flag for normal vs. afterburner generators.

Definition at line 87 of file GenModule.h.

87{this, "IsAfterburner", false, "Set true if generator modifies existing events rather than creating new ones"};

◆ m_lheFile

StringProperty Pythia8_i::m_lheFile {this, "LHEFile", ""}
privateinherited

Definition at line 133 of file Pythia8_i.h.

133{this, "LHEFile", ""};

◆ m_maxFailures

UnsignedIntegerProperty Pythia8_i::m_maxFailures {this, "MaxFailures", 10}
protectedinherited

Definition at line 90 of file Pythia8_i.h.

90{this, "MaxFailures", 10};

◆ m_mcEventKey

std::string GenBase::m_mcEventKey {}
protectedinherited

StoreGate key for the MC event collection (defaults to GEN_EVENT).

Definition at line 110 of file GenBase.h.

110{};

◆ m_mcevents_const

SG::ReadHandleKey<McEventCollection> GenBase::m_mcevents_const { this, "McEventKey", "GEN_EVENT", "StoreGate key of the MC event collection" }
privateinherited

Const handle to the MC event collection.

Definition at line 119 of file GenBase.h.

119{ this, "McEventKey", "GEN_EVENT", "StoreGate key of the MC event collection" };

◆ m_mergingWeight

double Pythia8_i::m_mergingWeight {1.0}
privateinherited

Definition at line 165 of file Pythia8_i.h.

165{1.0};

◆ m_mkMcEvent

BooleanProperty GenBase::m_mkMcEvent {this, "MakeMcEvent", false, "Create a new MC event collection if it doesn't exist"}
protectedinherited

Flag to determine if a new MC event collection should be made if it doesn't exist.

Definition at line 112 of file GenBase.h.

112{this, "MakeMcEvent", false, "Create a new MC event collection if it doesn't exist"};

◆ m_nAccepted

double Pythia8_i::m_nAccepted {0.}
privateinherited

Definition at line 138 of file Pythia8_i.h.

138{0.};

◆ m_nMerged

double Pythia8_i::m_nMerged {0.}
privateinherited

Definition at line 139 of file Pythia8_i.h.

139{0.};

◆ m_nSignalRequired

unsigned int Pythia8B_i::m_nSignalRequired
private

Definition at line 47 of file Pythia8B_i.h.

◆ m_numberAlphaS

DoubleProperty Pythia8_i::m_numberAlphaS {this,"numberAlphaS", 3.0}
protectedinherited

Definition at line 101 of file Pythia8_i.h.

101{this,"numberAlphaS", 3.0};

◆ m_oppCharges

bool Pythia8B_i::m_oppCharges
private

Definition at line 54 of file Pythia8B_i.h.

◆ m_outputParticleDataFile

StringProperty Pythia8_i::m_outputParticleDataFile {this, "OutputParticleData", "ParticleData.local.xml"}
privateinherited

Definition at line 163 of file Pythia8_i.h.

163{this, "OutputParticleData", "ParticleData.local.xml"};

◆ m_override_transform_beamenergy

bool Pythia8_i::m_override_transform_beamenergy {false}
privateinherited

Definition at line 131 of file Pythia8_i.h.

131{false};

◆ m_particleDataFile

StringProperty Pythia8_i::m_particleDataFile {this, "ParticleData", ""}
privateinherited

Definition at line 162 of file Pythia8_i.h.

162{this, "ParticleData", ""};

◆ m_particleIDs

std::map<std::string, PDGID> Pythia8_i::m_particleIDs
privateinherited

Definition at line 145 of file Pythia8_i.h.

◆ m_passingTriggerCuts

int Pythia8B_i::m_passingTriggerCuts
private

Definition at line 51 of file Pythia8B_i.h.

◆ m_procPtr

std::shared_ptr<Pythia8::Sigma2Process> Pythia8_i::m_procPtr {}
privateinherited

Definition at line 152 of file Pythia8_i.h.

152{};

◆ m_pt0timesMPI

DoubleProperty Pythia8_i::m_pt0timesMPI {this,"pT0timesMPI", 1.0}
protectedinherited

Definition at line 100 of file Pythia8_i.h.

100{this,"pT0timesMPI", 1.0};

◆ m_pythia

std::unique_ptr<Pythia8::Pythia> Pythia8_i::m_pythia {}
protectedinherited

Definition at line 88 of file Pythia8_i.h.

88{};

◆ m_pythiaToHepMC

HepMC::Pythia8ToHepMC Pythia8_i::m_pythiaToHepMC
protectedinherited

Definition at line 89 of file Pythia8_i.h.

◆ m_qEtaCut

double Pythia8B_i::m_qEtaCut
private

Definition at line 55 of file Pythia8B_i.h.

◆ m_qPtCut

double Pythia8B_i::m_qPtCut
private

Definition at line 55 of file Pythia8B_i.h.

◆ m_randomSeed

IntegerProperty GenModule::m_randomSeed {this, "RandomSeed", 1234567, "Random seed for the built-in random engine"}
protectedinherited

Seed for random number engine.

Definition at line 84 of file GenModule.h.

84{this, "RandomSeed", 1234567, "Random seed for the built-in random engine"}; // FIXME make this into an unsigned long int?

◆ m_rndmSvc

ServiceHandle<IAthRNGSvc> GenModule::m_rndmSvc {this, "RndmSvc", "AthRNGSvc"}
privateinherited

Data members.

Definition at line 97 of file GenModule.h.

97{this, "RndmSvc", "AthRNGSvc"};

◆ m_runinfo

std::shared_ptr<HepMC3::GenRunInfo> GenModule::m_runinfo {}
protectedinherited

The run info for HepMC3.

Definition at line 90 of file GenModule.h.

90{};

◆ m_sameAlphaSAsMPI

BooleanProperty Pythia8_i::m_sameAlphaSAsMPI {this,"useSameAlphaSasMPI", false}
protectedinherited

Definition at line 102 of file Pythia8_i.h.

102{this,"useSameAlphaSasMPI", false};

◆ m_saveLHE

BooleanProperty Pythia8_i::m_saveLHE {this, "SaveLHERecord", false}
privateinherited

Definition at line 176 of file Pythia8_i.h.

176{this, "SaveLHERecord", false};

◆ m_seeds

std::vector<long int> Pythia8B_i::m_seeds
private

Definition at line 53 of file Pythia8B_i.h.

◆ m_selectBQuarks

bool Pythia8B_i::m_selectBQuarks
private

Definition at line 54 of file Pythia8B_i.h.

◆ m_selectCQuarks

bool Pythia8B_i::m_selectCQuarks
private

Definition at line 54 of file Pythia8B_i.h.

◆ m_showerWeightNames

std::vector<std::string> Pythia8_i::m_showerWeightNames {"nominal"}
privateinherited

Definition at line 171 of file Pythia8_i.h.

171{"nominal"};

◆ m_showerWeightNamesProp

StringArrayProperty Pythia8_i::m_showerWeightNamesProp {this, "ShowerWeightNames", {} }
privateinherited

Definition at line 172 of file Pythia8_i.h.

172{this, "ShowerWeightNames", {} };

◆ m_sigCodes

std::vector<int> Pythia8B_i::m_sigCodes
private

Definition at line 48 of file Pythia8B_i.h.

◆ m_sigEtaCuts

std::vector<double> Pythia8B_i::m_sigEtaCuts
private

Definition at line 50 of file Pythia8B_i.h.

◆ m_sigmaTotal

double Pythia8_i::m_sigmaTotal {0.}
privateinherited

Definition at line 140 of file Pythia8_i.h.

140{0.};

◆ m_sigPtCuts

std::vector<double> Pythia8B_i::m_sigPtCuts
private

Definition at line 50 of file Pythia8B_i.h.

◆ m_speciesCount

std::map<int,int> Pythia8B_i::m_speciesCount
private

Definition at line 52 of file Pythia8B_i.h.

◆ m_SuppressSmallPT

Pythia8::SuppressSmallPT* Pythia8B_i::m_SuppressSmallPT
private

Definition at line 61 of file Pythia8B_i.h.

◆ m_totalBBarQuark

int Pythia8B_i::m_totalBBarQuark
private

Definition at line 51 of file Pythia8B_i.h.

◆ m_totalBQuark

int Pythia8B_i::m_totalBQuark
private

Definition at line 51 of file Pythia8B_i.h.

◆ m_totalCBarQuark

int Pythia8B_i::m_totalCBarQuark
private

Definition at line 51 of file Pythia8B_i.h.

◆ m_totalClone

int Pythia8B_i::m_totalClone
private

Definition at line 51 of file Pythia8B_i.h.

◆ m_totalCQuark

int Pythia8B_i::m_totalCQuark
private

Definition at line 51 of file Pythia8B_i.h.

◆ m_totalHard

int Pythia8B_i::m_totalHard
private

Definition at line 51 of file Pythia8B_i.h.

◆ m_totalPythiaCalls

int Pythia8B_i::m_totalPythiaCalls
private

Definition at line 51 of file Pythia8B_i.h.

◆ m_trigCode

int Pythia8B_i::m_trigCode
private

Definition at line 46 of file Pythia8B_i.h.

◆ m_trigEtaCut

double Pythia8B_i::m_trigEtaCut
private

Definition at line 55 of file Pythia8B_i.h.

◆ m_trigPtCut

std::vector<double> Pythia8B_i::m_trigPtCut
private

Definition at line 49 of file Pythia8B_i.h.

◆ m_useLHAPDF

BooleanProperty Pythia8_i::m_useLHAPDF {this, "UseLHAPDF", true}
privateinherited

Definition at line 160 of file Pythia8_i.h.

160{this, "UseLHAPDF", true};

◆ m_useReseed

BooleanProperty Pythia8_i::m_useReseed {this,"useReseed", false}
protectedinherited

Definition at line 95 of file Pythia8_i.h.

95{this,"useReseed", false};

◆ m_userHooks

StringArrayProperty Pythia8_i::m_userHooks {this, "UserHooks", {} }
protectedinherited

Definition at line 98 of file Pythia8_i.h.

98{this, "UserHooks", {} };

◆ m_userHooksPtrs

std::vector<UserHooksPtrType> Pythia8_i::m_userHooksPtrs {}
privateinherited

Definition at line 154 of file Pythia8_i.h.

154{};

◆ m_userModes

std::vector<std::string> Pythia8_i::m_userModes
privateinherited

Definition at line 122 of file Pythia8_i.h.

◆ m_useRndmGenSvc

BooleanProperty Pythia8_i::m_useRndmGenSvc {this, "useRndmGenSvc", true, "the max number of consecutive failures"}
protectedinherited

Definition at line 92 of file Pythia8_i.h.

92{this, "useRndmGenSvc", true, "the max number of consecutive failures"};

◆ m_userParams

std::vector<std::string> Pythia8_i::m_userParams
privateinherited

Definition at line 121 of file Pythia8_i.h.

◆ m_userProcess

StringProperty Pythia8_i::m_userProcess {this, "UserProcess", ""}
privateinherited

Definition at line 149 of file Pythia8_i.h.

149{this, "UserProcess", ""};

◆ m_userResonancePtrs

std::vector<std::shared_ptr<Pythia8::ResonanceWidths> > Pythia8_i::m_userResonancePtrs
privateinherited

Definition at line 158 of file Pythia8_i.h.

◆ m_userResonances

StringProperty Pythia8_i::m_userResonances {this, "UserResonances", ""}
privateinherited

Definition at line 156 of file Pythia8_i.h.

156{this, "UserResonances", ""};

◆ m_userString

std::string Pythia8B_i::m_userString
private

Definition at line 56 of file Pythia8B_i.h.

◆ m_userVar

std::vector<double> Pythia8B_i::m_userVar
private

Definition at line 57 of file Pythia8B_i.h.

◆ m_varHandleArraysDeclared

bool AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::m_varHandleArraysDeclared
privateinherited

Definition at line 399 of file AthCommonDataStore.h.

◆ m_version

double Pythia8_i::m_version {-1.}
privateinherited

Definition at line 110 of file Pythia8_i.h.

110{-1.};

◆ m_vetoDoubleB

bool Pythia8B_i::m_vetoDoubleB
private

Definition at line 54 of file Pythia8B_i.h.

◆ m_vetoDoubleC

bool Pythia8B_i::m_vetoDoubleC
private

Definition at line 54 of file Pythia8B_i.h.

◆ m_vhka

std::vector<SG::VarHandleKeyArray*> AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::m_vhka
privateinherited

Definition at line 398 of file AthCommonDataStore.h.

◆ m_weightCommands

std::vector<std::string> Pythia8_i::m_weightCommands {}
privateinherited

Definition at line 170 of file Pythia8_i.h.

170{};

◆ m_weightIDs

std::vector<std::string> Pythia8_i::m_weightIDs {}
privateinherited

Definition at line 167 of file Pythia8_i.h.

167{};

◆ m_weightNames

std::vector<std::string> Pythia8_i::m_weightNames {}
privateinherited

Definition at line 168 of file Pythia8_i.h.

168{};

◆ p_rndmEngine

CLHEP::HepRandomEngine * Pythia8B_i::p_rndmEngine = nullptr
static

Definition at line 41 of file Pythia8B_i.h.


The documentation for this class was generated from the following files: