diff --git a/src/EvtGenExternal/EvtHmePythia.cpp b/src/EvtGenExternal/EvtHmePythia.cpp index 55a9013..002d53c 100644 --- a/src/EvtGenExternal/EvtHmePythia.cpp +++ b/src/EvtGenExternal/EvtHmePythia.cpp @@ -1,90 +1,91 @@ /*********************************************************************** * Copyright 1998-2020 CERN for the benefit of the EvtGen authors * * * * This file is part of EvtGen. * * * * EvtGen is free software: you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation, either version 3 of the License, or * * (at your option) any later version. * * * * EvtGen is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with EvtGen. If not, see . * ***********************************************************************/ #include "EvtGenExternal/EvtHmePythia.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenExternal/EvtExternalGenFactory.hh" #include "EvtGenExternal/EvtPythiaEngine.hh" // Name std::string EvtHmePythia::getName() { return "HMEPYTHIA"; } // Clone the model EvtDecayBase* EvtHmePythia::clone() { return new EvtHmePythia(); } // Initialize the model void EvtHmePythia::init() { // Retrieve the Pythia engine if ( !_pythiaEngine ) { _pythiaEngine = dynamic_cast( EvtExternalGenFactory::getInstance()->getGenerator( EvtExternalGenFactory::PythiaGenId ) ); } // Check the matrix element is specified and is available. // Allowed values are (all implicitly include the tau neutrino): // 1521 = HMETau2Meson // 1531 = HMETau2TwoLeptons // 1532 = HMETau2TwoMesonsViaVector // 1533 = HMETau2TwoMesonsViaVectorScalar // 1541 = HMETau2ThreePions // 1542 = HMETau2ThreeMesonsWithKaons // 1543 = HMETau2ThreeMesonsGeneric // 1544 = HMETau2TwoPionsGamma // 1551 = HMETau2FourPions // 1561 = HMETau2FivePions checkNArg( 1 ); _hmeMode = getArg( 0 ); if ( _pythiaEngine->hmeAmplitudes( _hmeMode ) == 0 ) { EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Helicity matrix element mode " << _hmeMode << " is not available." << std::endl; + ::abort(); } // Check the number of daughters matches the matrix element // e.g. tau- -> e- anti-nu_e nu_tau, HMETau2TwoLeptons (1531), 3 daughters checkNDaug( ( _hmeMode % 100 ) / 10 ); } void EvtHmePythia::initProbMax() { setProbMax( _pythiaEngine->hmeAmplitudes( _hmeMode ) ); } void EvtHmePythia::decay( EvtParticle* p ) { // Initialize the phase-space p->initializePhaseSpace( getNDaug(), getDaugs() ); // Call the Pythia helicity matrix element if ( _pythiaEngine ) { _pythiaEngine->hmeAmplitudes( _hmeMode, p, this ); } } diff --git a/src/EvtGenExternal/EvtPythiaEngine.cpp b/src/EvtGenExternal/EvtPythiaEngine.cpp index ebe314d..60261b0 100644 --- a/src/EvtGenExternal/EvtPythiaEngine.cpp +++ b/src/EvtGenExternal/EvtPythiaEngine.cpp @@ -1,952 +1,972 @@ /*********************************************************************** * Copyright 1998-2020 CERN for the benefit of the EvtGen authors * * * * This file is part of EvtGen. * * * * EvtGen is free software: you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation, either version 3 of the License, or * * (at your option) any later version. * * * * EvtGen is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with EvtGen. If not, see . * ***********************************************************************/ #ifdef EVTGEN_PYTHIA #include "EvtGenExternal/EvtPythiaEngine.hh" #include "EvtGenBase/EvtDecayTable.hh" #include "EvtGenBase/EvtExtGeneratorCommandsTable.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtParticleFactory.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtSpinType.hh" #include "EvtGenExternal/EvtPythia6CommandConverter.hh" #include "Pythia8/Event.h" #include "Pythia8/ParticleData.h" #include #include #include using std::endl; EvtPythiaEngine::EvtPythiaEngine( std::string xmlDir, bool convertPhysCodes, bool useEvtGenRandom ) { // Create two Pythia generators. One will be for generic // Pythia decays in the decay.dec file. The other one will be to // only decay aliased particles, which are in general "signal" // decays different from those in the decay.dec file. // Even though it is not possible to have two different particle // versions in one Pythia generator, we can use two generators to // get the required behaviour. EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Creating generic Pythia generator" << endl; _genericPythiaGen = std::make_unique( xmlDir ); EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Creating alias Pythia generator" << endl; _aliasPythiaGen = std::make_unique( xmlDir, false ); _thePythiaGenerator = 0; _daugPDGVector.clear(); _daugP4Vector.clear(); _convertPhysCodes = convertPhysCodes; // Specify if we are going to use the random number generator (engine) // from EvtGen for Pythia 8. _useEvtGenRandom = useEvtGenRandom; _evtgenRandom = std::make_unique(); _initialised = false; } EvtPythiaEngine::~EvtPythiaEngine() { _thePythiaGenerator = nullptr; this->clearDaughterVectors(); this->clearPythiaModeMap(); } void EvtPythiaEngine::clearDaughterVectors() { _daugPDGVector.clear(); _daugP4Vector.clear(); } void EvtPythiaEngine::clearPythiaModeMap() { PythiaModeMap::iterator iter; for ( iter = _pythiaModeMap.begin(); iter != _pythiaModeMap.end(); ++iter ) { std::vector modeVector = iter->second; modeVector.clear(); } _pythiaModeMap.clear(); } void EvtPythiaEngine::initialise() { if ( _initialised ) { return; } this->clearPythiaModeMap(); this->updateParticleLists(); // Hadron-level processes only (hadronized, string fragmentation and secondary decays). // We do not want to generate the full pp or e+e- event structure etc.. _genericPythiaGen->readString( "ProcessLevel:all = off" ); _aliasPythiaGen->readString( "ProcessLevel:all = off" ); // Turn off Pythia warnings, e.g. changes to particle properties _genericPythiaGen->readString( "Print:quiet = on" ); _aliasPythiaGen->readString( "Print:quiet = on" ); // Apply any other physics (or special particle) requirements/cuts etc.. this->updatePhysicsParameters(); // Set the random number generator if ( _useEvtGenRandom == true ) { _genericPythiaGen->setRndmEnginePtr( _evtgenRandom.get() ); _aliasPythiaGen->setRndmEnginePtr( _evtgenRandom.get() ); } _genericPythiaGen->init(); _aliasPythiaGen->init(); // Create the tau decay model map _hmes[0] = std::make_unique(); _hmes[1521] = std::make_unique(); _hmes[1531] = std::make_unique(); _hmes[1532] = std::make_unique(); _hmes[1533] = std::make_unique(); _hmes[1541] = std::make_unique(); _hmes[1542] = std::make_unique(); _hmes[1543] = std::make_unique(); _hmes[1544] = std::make_unique(); _hmes[1551] = std::make_unique(); _hmes[1561] = std::make_unique(); // Initialize the tau decay model pointers for ( auto& model : _hmes ) { model.second->initPointers( &_genericPythiaGen->particleData, &_genericPythiaGen->couplings ); } _initialised = true; } bool EvtPythiaEngine::doDecay( EvtParticle* theParticle ) { // Store the mother particle within a Pythia8 Event object. // Then do the hadron level decays. // The EvtParticle must be a colour singlet (meson/baryon/lepton), i.e. not a gluon or quark // We delete any daughters the particle may have, since we are asking Pythia // to generate the decay anew. Also note that _any_ Pythia decay allowed for the particle // will be generated and not the specific Pythia decay mode that EvtGen has already // specified. This is necessary since we only want to initialise the Pythia decay table // once; all Pythia branching fractions for a given mother particle are renormalised to sum to 1.0. // In EvtGen decay.dec files, it may be the case that Pythia decays are only used // for some of the particle decays (i.e. Pythia BF sum < 1.0). As we loop over many events, // the total frequency for each Pythia decay mode will normalise correctly to what // we wanted via the specifications made to the decay.dec file, even though event-by-event // the EvtGen decay channel and the Pythia decay channel may be different. if ( _initialised == false ) { this->initialise(); } if ( theParticle == 0 ) { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Error in EvtPythiaEngine::doDecay. The mother particle is null. Not doing any Pythia decay." << endl; return false; } // Delete EvtParticle daughters if they already exist if ( theParticle->getNDaug() != 0 ) { bool keepChannel( false ); theParticle->deleteDaughters( keepChannel ); } EvtId particleId = theParticle->getId(); int isAlias = particleId.isAlias(); // Choose the generator depending if we have an aliased (parent) particle or not _thePythiaGenerator = ( isAlias == 1 ? _aliasPythiaGen.get() : _genericPythiaGen.get() ); // Need to use the reference to the Pythia8::Event object, // otherwise it will just return a new empty, default event object. Pythia8::Event& theEvent = _thePythiaGenerator->event; theEvent.reset(); // Initialise the event to be the particle rest frame int PDGCode = EvtPDL::getStdHep( particleId ); int status( 1 ); int colour( 0 ), anticolour( 0 ); double px( 0.0 ), py( 0.0 ), pz( 0.0 ); double m0 = theParticle->mass(); double E = m0; theEvent.append( PDGCode, status, colour, anticolour, px, py, pz, E, m0 ); // Generate the Pythia event int iTrial( 0 ); bool generatedEvent( false ); for ( iTrial = 0; iTrial < 10; iTrial++ ) { generatedEvent = _thePythiaGenerator->next(); if ( generatedEvent ) { break; } } bool success( false ); if ( generatedEvent ) { // Store the daughters for this particle from the Pythia decay tree. // This is a recursive function that will continue looping through // all available daughters until the first set of non-quark and non-gluon // particles are encountered in the Pythia Event structure. // First, clear up the internal vectors storing the daughter // EvtId types and 4-momenta. this->clearDaughterVectors(); // Now store the daughter info. Since this is a recursive function // to loop through the full Pythia decay tree, we do not want to create // EvtParticles here but in the next step. this->storeDaughterInfo( theParticle, 1 ); // Now create the EvtParticle daughters of the (parent) particle. // We need to use the EvtParticle::makeDaughters function // owing to the way EvtParticle stores parent-daughter information. this->createDaughterEvtParticles( theParticle ); //theParticle->printTree(); //theEvent.list(true, true); success = true; } return success; } void EvtPythiaEngine::storeDaughterInfo( EvtParticle* theParticle, int startInt ) { Pythia8::Event& theEvent = _thePythiaGenerator->event; std::vector daugList = theEvent.daughterList( startInt ); std::vector::iterator daugIter; for ( daugIter = daugList.begin(); daugIter != daugList.end(); ++daugIter ) { int daugInt = *daugIter; // Ask if the daughter is a quark or gluon. If so, recursively search again. Pythia8::Particle& daugParticle = theEvent[daugInt]; if ( daugParticle.isQuark() || daugParticle.isGluon() ) { // Recursively search for correct daughter type this->storeDaughterInfo( theParticle, daugInt ); } else { // We have a daughter that is not a quark nor gluon particle. // Make sure we are not double counting particles, since several quarks // and gluons make one particle. // Set the status flag for any "new" particle to say that we have stored it. // Use status flag = 1000 (within the user allowed range for Pythia codes). // Check that the status flag for the particle is not equal to 1000 int statusCode = daugParticle.status(); if ( statusCode != 1000 ) { int daugPDGInt = daugParticle.id(); double px = daugParticle.px(); double py = daugParticle.py(); double pz = daugParticle.pz(); double E = daugParticle.e(); EvtVector4R daughterP4( E, px, py, pz ); // Now store the EvtId and 4-momentum in the internal vectors _daugPDGVector.push_back( daugPDGInt ); _daugP4Vector.push_back( daughterP4 ); // Set the status flag for the Pythia particle to let us know // that we have already considered it to avoid double counting. daugParticle.status( 1000 ); } // Status code != 1000 } } } void EvtPythiaEngine::createDaughterEvtParticles( EvtParticle* theParent ) { if ( theParent == 0 ) { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Error in EvtPythiaEngine::createDaughterEvtParticles. The parent is null" << endl; return; } // Get the list of Pythia decay modes defined for this particle id alias. // It would be easier to just use the decay channel number that Pythia chose to use // for the particle decay, but this is not accessible from the Pythia interface at present. int nDaughters = _daugPDGVector.size(); std::vector daugAliasIdVect( 0 ); EvtId particleId = theParent->getId(); // Check to see if we have an anti-particle. If we do, charge conjugate the particle id to get the // Pythia "alias" we can compare with the defined (particle) Pythia modes. int PDGId = EvtPDL::getStdHep( particleId ); int aliasInt = particleId.getAlias(); int pythiaAliasInt( aliasInt ); if ( PDGId < 0 ) { // We have an anti-particle. EvtId conjPartId = EvtPDL::chargeConj( particleId ); pythiaAliasInt = conjPartId.getAlias(); } std::vector pythiaModes = _pythiaModeMap[pythiaAliasInt]; // Loop over all available Pythia decay modes and find the channel that matches // the daughter ids. Set each daughter id to also use the alias integer. // This will then convert the Pythia generated channel to the EvtGen alias defined one. std::vector::iterator modeIter; bool gotMode( false ); for ( modeIter = pythiaModes.begin(); modeIter != pythiaModes.end(); ++modeIter ) { // Stop the loop if we have the right decay mode channel if ( gotMode ) { break; } int pythiaModeInt = *modeIter; EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel( aliasInt, pythiaModeInt ); if ( decayModel != 0 ) { int nModeDaug = decayModel->getNDaug(); // We need to make sure that the number of daughters match if ( nDaughters == nModeDaug ) { int iModeDaug( 0 ); for ( iModeDaug = 0; iModeDaug < nModeDaug; iModeDaug++ ) { EvtId daugId = decayModel->getDaug( iModeDaug ); int daugPDGId = EvtPDL::getStdHep( daugId ); // Pythia has used the right PDG codes for this decay mode, even for conjugate modes int pythiaPDGId = _daugPDGVector[iModeDaug]; if ( daugPDGId == pythiaPDGId ) { daugAliasIdVect.push_back( daugId ); } } // Loop over EvtGen mode daughters int daugAliasSize = daugAliasIdVect.size(); if ( daugAliasSize == nDaughters ) { // All daughter Id codes are accounted for. Set the flag to stop the loop. gotMode = true; } else { // We do not have the correct daughter ordering. Clear the id vector // and try another mode. daugAliasIdVect.clear(); } } // Same number of daughters } // decayModel != 0 } // Loop over available Pythia modes if ( gotMode == false ) { // We did not find a match for the daughter aliases. Just use the normal PDG codes // from the Pythia decay result int iPyDaug( 0 ); for ( iPyDaug = 0; iPyDaug < nDaughters; iPyDaug++ ) { int daugPDGCode = _daugPDGVector[iPyDaug]; EvtId daugPyId = EvtPDL::evtIdFromStdHep( daugPDGCode ); daugAliasIdVect.push_back( daugPyId ); } } // Make the EvtParticle daughters (with correct alias id's). Their 4-momenta are uninitialised. theParent->makeDaughters( nDaughters, daugAliasIdVect ); // Now set the 4-momenta of the daughters. int iDaug( 0 ); // Can use an iterator here, but we already had to use the vector size... for ( iDaug = 0; iDaug < nDaughters; iDaug++ ) { EvtParticle* theDaughter = theParent->getDaug( iDaug ); // Set the correct 4-momentum for each daughter particle. if ( theDaughter != 0 ) { EvtId theDaugId = daugAliasIdVect[iDaug]; const EvtVector4R theDaugP4 = _daugP4Vector[iDaug]; theDaughter->init( theDaugId, theDaugP4 ); } } } void EvtPythiaEngine::updateParticleLists() { // Use the EvtGen decay table (decay/user.dec) to update Pythia particle // decay modes. Also, make sure the generic and alias Pythia generators // use the same particle data entries as defined by EvtGen (evt.pdl). // Loop over all entries in the EvtPDL particle data table. // Aliases are added at the end with id numbers equal to the // original particle, but with alias integer = lastPDLEntry+1 etc.. int iPDL; int nPDL = EvtPDL::entries(); // Reset the _addedPDGCodes map that keeps track // of any new particles added to the Pythia input data stream _addedPDGCodes.clear(); for ( iPDL = 0; iPDL < nPDL; iPDL++ ) { EvtId particleId = EvtPDL::getEntry( iPDL ); int aliasInt = particleId.getAlias(); // Pythia and EvtGen should use the same PDG codes for particles int PDGCode = EvtPDL::getStdHep( particleId ); // Update pole mass, width, lifetime and mass range double mass = EvtPDL::getMeanMass( particleId ); double width = EvtPDL::getWidth( particleId ); double lifetime = EvtPDL::getctau( particleId ); double mmin = EvtPDL::getMinMass( particleId ); double mmax = EvtPDL::getMaxMass( particleId ); // Particle data pointers. The generic and alias Pythia generator pointers have // their own Pythia8::ParticleData particleData public data members which contain // the particle properties table. We can extract and change individual particle // entries using the particleDataEntryPtr() function within ParticleData. // However, we must be careful when accessing the particleData table. We must do // this directly, since assigning it to another Pythia8::ParticleData object via copying // or assignment will give it a different memory address and it will no longer refer to // the original particleData information from the generator pointer. Pythia8::ParticleDataEntry* entry_generic = _genericPythiaGen->particleData.particleDataEntryPtr( PDGCode ); Pythia8::ParticleDataEntry* entry_alias = _aliasPythiaGen->particleData.particleDataEntryPtr( PDGCode ); // Check that the PDG code is not zero/null and exclude other // special cases, e.g. those reserved for internal generator use if ( entry_generic != 0 && this->validPDGCode( PDGCode ) ) { entry_generic->setM0( mass ); entry_generic->setMWidth( width ); entry_generic->setTau0( lifetime ); if ( std::fabs( width ) > 0.0 ) { entry_generic->setMMin( mmin ); entry_generic->setMMax( mmax ); } } // Check that the PDG code is not zero/null and exclude other // special cases, e.g. those reserved for internal generator use if ( entry_alias != 0 && this->validPDGCode( PDGCode ) ) { entry_alias->setM0( mass ); entry_alias->setMWidth( width ); entry_alias->setTau0( lifetime ); if ( std::fabs( width ) > 0.0 ) { entry_alias->setMMin( mmin ); entry_alias->setMMax( mmax ); } } // Check which particles have a Pythia decay defined. // Get the list of all possible decays for the particle, using the alias integer. // If the particle is not actually an alias, aliasInt = idInt. bool hasPythiaDecays = EvtDecayTable::getInstance()->hasPythia( aliasInt ); if ( hasPythiaDecays ) { int isAlias = particleId.isAlias(); // Decide what generator to use depending on whether we have // an aliased particle or not _thePythiaGenerator = ( isAlias == 1 ? _aliasPythiaGen.get() : _genericPythiaGen.get() ); // Find the Pythia particle name given the standard PDG code integer std::string dataName = _thePythiaGenerator->particleData.name( PDGCode ); bool alreadyStored = ( _addedPDGCodes.find( abs( PDGCode ) ) != _addedPDGCodes.end() ); if ( dataName == " " && !alreadyStored ) { // Particle and its antiparticle do not exist in the Pythia database. // Create a new particle, then create the new decay modes. this->createPythiaParticle( particleId, PDGCode ); } // For the particle, create the Pythia decay modes. // Update Pythia data tables. this->updatePythiaDecayTable( particleId, aliasInt, PDGCode ); } // Loop over Pythia decays } // Loop over EvtPDL entries //EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Writing out changed generic Pythia decay list"<particleData.listChanged(); //EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Writing out changed alias Pythia decay list"<particleData.listChanged(); } bool EvtPythiaEngine::validPDGCode( int PDGCode ) { // Exclude certain PDG codes: void = 0 and special values = 81 to 100, which are reserved // for internal generator use (pseudoparticles) according to PDG guidelines. Also exclude // nu'_tau (nu_L) = 18, which has different masses: Pythia8 = 400 GeV, EvtGen = 0 GeV. bool isValid( true ); int absPDGCode = abs( PDGCode ); if ( absPDGCode == 0 || absPDGCode == 18 ) { // Void and nu_L or nu'_tau isValid = false; } else if ( absPDGCode >= 81 && absPDGCode <= 100 ) { // Pseudoparticles isValid = false; } return isValid; } void EvtPythiaEngine::updatePythiaDecayTable( EvtId& particleId, int aliasInt, int PDGCode ) { // Update the particle data table in Pythia. // The tables store information about the allowed decay modes // where the PDGId for all particles must be positive; anti-particles are stored // with the corresponding particle entry. // Since we do not want to implement CP violation here, just use the same branching // fractions for particle and anti-particle modes. int nModes = EvtDecayTable::getInstance()->getNModes( aliasInt ); int iMode( 0 ); bool firstMode( true ); // Only process positive PDG codes. if ( PDGCode < 0 ) { return; } // Keep track of which decay modes are Pythia decays for each aliasInt std::vector pythiaModes( 0 ); // Loop over the decay modes for this particle for ( iMode = 0; iMode < nModes; iMode++ ) { EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel( aliasInt, iMode ); if ( decayModel != 0 ) { int nDaug = decayModel->getNDaug(); // If the decay mode has no daughters, then that means that there will be // no entries for any submode re-definitions for Pythia. // This sometimes occurs for any mode using non-standard Pythia 6 codes. // Do not refine the decay mode, i.e. accept the Pythia 8 default (if it exists). if ( nDaug > 0 ) { // Check to see if we have a Pythia decay mode std::string modelName = decayModel->getModelName(); if ( modelName == "PYTHIA" ) { // Keep track which decay mode is a Pythia one. We need this in order to // reassign alias Id values for particles generated in the decay. pythiaModes.push_back( iMode ); std::ostringstream oss; oss.setf( std::ios::scientific ); // Write out the absolute value of the PDG code, since // particles and anti-particles occupy the same part of the Pythia table. oss << PDGCode; if ( firstMode ) { // Create a new channel oss << ":oneChannel = "; firstMode = false; } else { // Add the channel oss << ":addChannel = "; } // Select all channels (particle and anti-particle). // For CP violation, or different BFs for particle and anti-particle, // use options 2 or 3 (not here). int onMode( 1 ); oss << onMode << " "; double BF = decayModel->getBranchingFraction(); oss << BF << " "; // Need to convert the old Pythia physics mode integers with the new ones // To do this, get the model argument and write a conversion method. int modeInt = this->getModeInt( decayModel ); oss << modeInt; int iDaug( 0 ); for ( iDaug = 0; iDaug < nDaug; iDaug++ ) { EvtId daugId = decayModel->getDaug( iDaug ); int daugPDG = EvtPDL::getStdHep( daugId ); oss << " " << daugPDG; } // Daughter list _thePythiaGenerator->readString( oss.str() ); } // is Pythia } else { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Warning in EvtPythiaEngine. Trying to redefine Pythia table for " << EvtPDL::name( particleId ) << " for a decay channel that has no daughters." << " Keeping Pythia default (if available)." << endl; } } else { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Error in EvtPythiaEngine. decayModel is null for particle " << EvtPDL::name( particleId ) << " mode number " << iMode << endl; } } // Loop over modes _pythiaModeMap[aliasInt] = pythiaModes; // Now, renormalise the decay branching fractions to sum to 1.0 std::ostringstream rescaleStr; rescaleStr.setf( std::ios::scientific ); rescaleStr << PDGCode << ":rescaleBR = 1.0"; _thePythiaGenerator->readString( rescaleStr.str() ); } int EvtPythiaEngine::getModeInt( EvtDecayBase* decayModel ) { int tmpModeInt( 0 ), modeInt( 0 ); if ( decayModel != 0 ) { int nVars = decayModel->getNArg(); // Just read the first integer, which specifies the Pythia decay model. // Ignore any other values. if ( nVars > 0 ) { tmpModeInt = static_cast( decayModel->getArg( 0 ) ); } } if ( _convertPhysCodes ) { // Extra code to convert the old Pythia decay model integer MDME(ICC,2) to the new one. // This should be removed eventually after updating decay.dec files to use // the new convention. if ( tmpModeInt == 0 ) { modeInt = 0; // phase-space } else if ( tmpModeInt == 1 ) { modeInt = 1; // omega or phi -> 3pi } else if ( tmpModeInt == 2 ) { modeInt = 11; // Dalitz decay } else if ( tmpModeInt == 3 ) { modeInt = 2; // V -> PS PS } else if ( tmpModeInt == 4 ) { modeInt = 92; // onium -> ggg or gg gamma } else if ( tmpModeInt == 11 ) { modeInt = 42; // phase-space of hadrons from available quarks } else if ( tmpModeInt == 12 ) { modeInt = 42; // phase-space for onia resonances } else if ( tmpModeInt == 13 ) { modeInt = 43; // phase-space of at least 3 hadrons } else if ( tmpModeInt == 14 ) { modeInt = 44; // phase-space of at least 4 hadrons } else if ( tmpModeInt == 15 ) { modeInt = 45; // phase-space of at least 5 hadrons } else if ( tmpModeInt >= 22 && tmpModeInt <= 30 ) { modeInt = tmpModeInt + 40; // phase space of hadrons with fixed multiplicity (modeInt - 60) } else if ( tmpModeInt == 31 ) { modeInt = 42; // two or more quarks phase-space; one spectactor quark } else if ( tmpModeInt == 32 ) { modeInt = 91; // qqbar or gg pair } else if ( tmpModeInt == 33 ) { modeInt = 0; // triplet q X qbar, where X = gluon or colour singlet (superfluous, since g's are created anyway) } else if ( tmpModeInt == 41 ) { modeInt = 21; // weak decay phase space, weighting nu_tau spectrum } else if ( tmpModeInt == 42 ) { modeInt = 22; // weak decay V-A matrix element } else if ( tmpModeInt == 43 ) { modeInt = 22; // weak decay V-A matrix element, quarks as jets (superfluous) } else if ( tmpModeInt == 44 ) { modeInt = 22; // weak decay V-A matrix element, parton showers (superfluous) } else if ( tmpModeInt == 48 ) { modeInt = 23; // weak decay V-A matrix element, at least 3 decay products } else if ( tmpModeInt == 50 ) { modeInt = 0; // default behaviour } else if ( tmpModeInt == 51 ) { modeInt = 0; // step threshold (channel switched off when mass daughters > mother mass } else if ( tmpModeInt == 52 || tmpModeInt == 53 ) { modeInt = 0; // beta-factor threshold } else if ( tmpModeInt == 84 ) { modeInt = 42; // unknown physics process - just use phase-space } else if ( tmpModeInt == 101 ) { modeInt = 0; // continuation line } else if ( tmpModeInt == 102 ) { modeInt = 0; // off mass shell particles. } else { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Pythia mode integer " << tmpModeInt << " is not recognised. Using phase-space model" << endl; modeInt = 0; // Use phase-space for anything else } } else { // No need to convert the physics mode integer code modeInt = tmpModeInt; } return modeInt; } void EvtPythiaEngine::createPythiaParticle( EvtId& particleId, int PDGCode ) { // Use the EvtGen name, PDGId and other variables to define the new Pythia particle. EvtId antiPartId = EvtPDL::chargeConj( particleId ); std::string aliasName = EvtPDL::name( particleId ); // If not an alias, aliasName = normal name std::string antiName = EvtPDL::name( antiPartId ); EvtSpinType::spintype spinType = EvtPDL::getSpinType( particleId ); int spin = EvtSpinType::getSpin2( spinType ); int charge = EvtPDL::chg3( particleId ); // Must set the correct colour type manually here, since the evt.pdl file // does not store this information. This is required for quarks otherwise // Pythia cannot generate the decay properly. int PDGId = EvtPDL::getStdHep( particleId ); int colour( 0 ); if ( PDGId == 21 ) { colour = 2; // gluons } else if ( PDGId <= 8 && PDGId > 0 ) { colour = 1; // single quarks } double m0 = EvtPDL::getMeanMass( particleId ); double mWidth = EvtPDL::getWidth( particleId ); double mMin = EvtPDL::getMinMass( particleId ); double mMax = EvtPDL::getMaxMass( particleId ); double tau0 = EvtPDL::getctau( particleId ); std::ostringstream oss; oss.setf( std::ios::scientific ); int absPDGCode = abs( PDGCode ); oss << absPDGCode << ":new = " << aliasName << " " << antiName << " " << spin << " " << charge << " " << colour << " " << m0 << " " << mWidth << " " << mMin << " " << mMax << " " << tau0; // Pass this information to Pythia _thePythiaGenerator->readString( oss.str() ); // Also store the absolute value of the PDG entry // to keep track of which new particles have been added, // which also automatically includes the anti-particle. // We need to avoid creating new anti-particles when // they already exist when the particle was added. _addedPDGCodes[absPDGCode] = 1; } void EvtPythiaEngine::updatePhysicsParameters() { // Update any more Pythia physics (or special particle) requirements/cuts etc.. // This should be used if any of the Pythia 6 parameters like JetSetPar MSTJ(i) = x // are needed. Such commands will need to be implemented using the new interface // pythiaGenerator->readString(cmd); Here cmd is a string telling Pythia 8 // what physics parameters to change. This will need to be done for the generic and // alias generator pointers, as appropriate. // Set the multiplicity level for hadronic weak decays std::string multiWeakCut( "ParticleDecays:multIncreaseWeak = 2.0" ); _genericPythiaGen->readString( multiWeakCut ); _aliasPythiaGen->readString( multiWeakCut ); // Set the multiplicity level for all other decays std::string multiCut( "ParticleDecays:multIncrease = 4.5" ); _genericPythiaGen->readString( multiCut ); _aliasPythiaGen->readString( multiCut ); //Now read in any custom configuration entered in the XML GeneratorCommands commands = EvtExtGeneratorCommandsTable::getInstance()->getCommands( "PYTHIA" ); GeneratorCommands::iterator it = commands.begin(); for ( ; it != commands.end(); it++ ) { Command command = *it; std::vector commandStrings; if ( command["VERSION"] == "PYTHIA6" ) { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Converting Pythia 6 command: " << command["MODULE"] << "(" << command["PARAM"] << ")=" << command["VALUE"] << "..." << endl; commandStrings = convertPythia6Command( command ); } else if ( command["VERSION"] == "PYTHIA8" ) { commandStrings.push_back( command["MODULE"] + ":" + command["PARAM"] + " = " + command["VALUE"] ); } else { EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Pythia command received by EvtPythiaEngine has bad version:" << endl; EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Received " << command["VERSION"] << " but expected PYTHIA6 or PYTHIA8." << endl; EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "The error is likely to be in EvtDecayTable.cpp" << endl; EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "EvtGen will now abort." << endl; ::abort(); } std::string generator = command["GENERATOR"]; if ( generator == "GENERIC" || generator == "Generic" || generator == "generic" || generator == "BOTH" || generator == "Both" || generator == "both" ) { std::vector::iterator it2 = commandStrings.begin(); for ( ; it2 != commandStrings.end(); it2++ ) { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Configuring generic Pythia generator: " << ( *it2 ) << endl; _genericPythiaGen->readString( *it2 ); } } if ( generator == "ALIAS" || generator == "Alias" || generator == "alias" || generator == "BOTH" || generator == "Both" || generator == "both" ) { std::vector::iterator it2 = commandStrings.begin(); for ( ; it2 != commandStrings.end(); it2++ ) { EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Configuring alias Pythia generator: " << ( *it2 ) << endl; _aliasPythiaGen->readString( *it2 ); } } } } double EvtPythiaEngine::hmeAmplitudes( int hmeMode, EvtParticle* mom, EvtDecayAmp* amp ) { + if ( _initialised == false ) { + this->initialise(); + } + // Get the helicity matrix element (return 0 if no matrix element) std::unique_ptr& hme = _hmes[hmeMode]; if ( hme == nullptr ) { return 0.0; } // Return the maximum amplitude if no particle is given std::vector prts; if ( mom == nullptr ) { - prts.push_back( Pythia8::HelicityParticle( mom->getPDGId() ) ); - prts[0].rho[0][0] = 1; + // Use tau particle + int PDGCode = EvtPDL::getStdHep( EvtPDL::getId( "tau-" ) ); + const Pythia8::Particle particle( PDGCode ); + prts.push_back( Pythia8::HelicityParticle( + particle, &_genericPythiaGen->particleData ) ); + + // Set first phase space rho element to unity + prts[0].rho[0][0] = 1.0; + + // Initialise channel constants, including the maximum decay weight hard limit. + // This usually requires the tau mass via Particle::m(), which is initialised + // in the Particle constructor. However this mass is set to zero since we + // only pass the PDG integer to the constructor, and m() is not automatically + // updated by the particle data. So we need to set it (to the pole mass): + double m0 = prts[0].m0(); + prts[0].m( m0 ); + hme->initChannel( prts ); return hme->decayWeightMax( prts ); } + int ih( 0 ); std::vector ihs( 3, -1 ), nhs( 3, 0 ), hs( mom->getNDaug() + 1, 0 ); // Translate the mother to HelicityParticle EvtVector4R p = mom->getP4Lab(); prts.push_back( Pythia8::HelicityParticle( mom->getPDGId(), 0, 0, 0, 0, 0, 0, 0, p.get( 1 ), p.get( 2 ), p.get( 3 ), p.get( 0 ), p.mass(), p.mass(), &_thePythiaGenerator->particleData ) ); prts.back().direction = -1; if ( mom->getSpinStates() > 1 ) { ihs[ih] = 0; nhs[ih] = mom->getSpinStates(); ++ih; } // Translate the children to HelicityParticle for ( int iDtr = 0; iDtr < (int)mom->getNDaug(); ++iDtr ) { EvtParticle* dtr = mom->getDaug( iDtr ); p = dtr->getP4Lab(); prts.push_back( Pythia8::HelicityParticle( dtr->getPDGId(), 0, 0, 0, 0, 0, 0, 0, p.get( 1 ), p.get( 2 ), p.get( 3 ), p.get( 0 ), p.mass(), p.mass(), &_thePythiaGenerator->particleData ) ); prts.back().direction = 1; if ( dtr->getSpinStates() > 1 ) { if ( ih < 3 ) { ihs[ih] = 0; nhs[ih] = dtr->getSpinStates(); ++ih; } else { EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "More than three particles with multiple spins." << endl; } } if ( dtr->getSpinType() == EvtSpinType::NEUTRINO && dtr->getPDGId() < 0 ) hs[iDtr + 1] = 1; } // Calculate and set the amplitudes (maximum three helicities) hme->initChannel( prts ); // The initWaves function in pythia8/include/Pythia8/HelicityMatrixElements.h // needs to be made public (from protected) before we can use it: // hme->initWaves( prts ); if ( ihs[0] < 0 ) { amp->vertex( hme->calculateME( hs ) ); } else { for ( hs[ihs[0]] = 0; hs[ihs[0]] < nhs[0]; ++hs[ihs[0]] ) { if ( ihs[1] < 0 ) { amp->vertex( hs[ihs[0]], hme->calculateME( hs ) ); } else { for ( hs[ihs[1]] = 0; hs[ihs[1]] < nhs[1]; ++hs[ihs[1]] ) { if ( ihs[2] < 0 ) { amp->vertex( hs[ihs[0]], hs[ihs[1]], hme->calculateME( hs ) ); } else { for ( hs[ihs[2]] = 0; hs[ihs[2]] < nhs[2]; ++hs[ihs[2]] ) { amp->vertex( hs[ihs[0]], hs[ihs[1]], hs[ihs[2]], hme->calculateME( hs ) ); } } // End third helicity } } // End second helicity } } // End first helicity // Return the decay weight return hme->decayWeight( prts ); } #endif diff --git a/test/jsonFiles/HmePythiaTau1.json b/test/jsonFiles/HmePythiaTau1.json new file mode 100644 index 0000000..fb4c71d --- /dev/null +++ b/test/jsonFiles/HmePythiaTau1.json @@ -0,0 +1,15 @@ +{ + "parent" : "tau-", + "daughters" : ["e-", "anti-nu_e", "nu_tau"], + "model" : "HMEPYTHIA", + "parameters" : [1531], + "events" : 10000, + "histograms" : [ + {"variable" : "mass", "title" : "mENuE" , "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "psumsq", "title" : "qSq", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "coshel", "title" : "cosE", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "prob", "title" : "Prob", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0} + ], + "outfile" : "HmePythiaTau1.root", + "reference" : "RefHmePythiaTau1.root" +} diff --git a/test/testDecayModel.cc b/test/testDecayModel.cc index 817e8c0..e968c8c 100644 --- a/test/testDecayModel.cc +++ b/test/testDecayModel.cc @@ -1,357 +1,348 @@ /*********************************************************************** * Copyright 1998-2020 CERN for the benefit of the EvtGen authors * * * * This file is part of EvtGen. * * * * EvtGen is free software: you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation, either version 3 of the License, or * * (at your option) any later version. * * * * EvtGen is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with EvtGen. If not, see . * ***********************************************************************/ #include "testDecayModel.hh" #include "EvtGen/EvtGen.hh" #include "EvtGenBase/EvtAbsRadCorr.hh" #include "EvtGenBase/EvtDecayBase.hh" #include "EvtGenBase/EvtId.hh" #include "EvtGenBase/EvtMTRandomEngine.hh" +#include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtParticleFactory.hh" -#include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtVector4R.hh" #ifdef EVTGEN_EXTERNAL #include "EvtGenExternal/EvtExternalGenList.hh" #endif #include #include -#include #include +#include #include testDecayModel::testDecayModel( const std::string& jsonFile ) : - m_jsonFile( jsonFile ), - m_infoVect(), - m_histVect(), - m_nHistos( 0 ) + m_jsonFile( jsonFile ), m_infoVect(), m_histVect(), m_nHistos( 0 ) { } -void testDecayModel::run() { - +void testDecayModel::run() +{ Json::Value config; std::ifstream inputStr( m_jsonFile.c_str() ); inputStr >> config; inputStr.close(); const std::string parentName = config["parent"].asString(); Json::Value jDaughters = config["daughters"]; std::vector daugNames; for ( auto daugN : jDaughters ) { - daugNames.push_back( daugN.asString() ); + daugNames.push_back( daugN.asString() ); } const std::string modelName = config["model"].asString(); Json::Value jParameters = config["parameters"]; std::vector modelPars; for ( auto par : jParameters ) { - modelPars.push_back( par.asDouble() ); + modelPars.push_back( par.asDouble() ); } // Create decay file based on json file input - const std::string decFile = createDecFile( parentName, daugNames, - modelName, modelPars ); + const std::string decFile = createDecFile( parentName, daugNames, modelName, + modelPars ); int nEvents = config["events"].asInt(); std::cout << "Number of events = " << nEvents << std::endl; // ROOT output file const std::string outFileName = config["outfile"].asString(); TFile* outFile = TFile::Open( outFileName.c_str(), "recreate" ); // Define histograms defineHistos( config, outFile ); // Generate events and fill histograms generateEvents( decFile, parentName, nEvents ); // Normalize histograms for ( auto hist : m_histVect ) { - double area = hist->Integral(); - if ( area > 0.0 ) { hist->Scale( 1.0/area ); } + double area = hist->Integral(); + if ( area > 0.0 ) { + hist->Scale( 1.0 / area ); + } } - + // Compare with reference histograms const std::string refFileName = config["reference"].asString(); compareHistos( refFileName ); - + // Write output outFile->cd(); for ( auto hist : m_histVect ) { - hist->Write(); + hist->Write(); } outFile->Close(); - } -std::string testDecayModel::createDecFile( const std::string& parent, - const std::vector& daugNames, - const std::string& modelName, - const std::vector& parameters ) const +std::string testDecayModel::createDecFile( + const std::string& parent, const std::vector& daugNames, + const std::string& modelName, const std::vector& parameters ) const { - // Create (or overwrite) the decay file - std::string decName( modelName ); decName += ".dec"; + std::string decName( modelName ); + decName += ".dec"; std::ofstream decFile; decFile.open( decName.c_str() ); decFile << "Decay " << parent << std::endl; decFile << "1.0"; for ( auto daug : daugNames ) { - decFile << " " << daug; + decFile << " " << daug; } decFile << " " << modelName; for ( auto par : parameters ) { - decFile << " " << par; + decFile << " " << par; } decFile << ";" << std::endl; decFile << "Enddecay" << std::endl; decFile << "End" << std::endl; decFile.close(); - + return decName; - } -void testDecayModel::defineHistos( Json::Value& config, TFile* outFile ) { - +void testDecayModel::defineHistos( Json::Value& config, TFile* outFile ) +{ // Histogram information Json::Value jHistos = config["histograms"]; size_t nHistos = jHistos.size(); - + m_infoVect.reserve( nHistos ); m_histVect.reserve( nHistos ); - - for ( auto hInfo : jHistos ) { - const std::string varName = hInfo["variable"].asString(); - // Integer values that define what particles need to be used - // for invariant mass combinations or helicity angles etc - int d1 = hInfo["d1"].asInt(); - int d2 = hInfo["d2"].asInt(); - - testInfo info( varName, d1, d2 ); - m_infoVect.push_back( info ); - - const std::string varTitle = hInfo["title"].asString(); - int nBins = hInfo["nbins"].asInt(); - double xmin = hInfo["xmin"].asDouble(); - double xmax = hInfo["xmax"].asDouble(); - - TString histName( varName.c_str() ); - if ( d1 != 0 ) { histName += "_"; histName += d1; } - if ( d2 != 0 ) { histName += "_"; histName += d2; } - TH1D* hist = new TH1D( histName.Data(), varTitle.c_str(), nBins, xmin, xmax ); - hist->SetDirectory( outFile ); - m_histVect.push_back( hist ); - - } + for ( auto hInfo : jHistos ) { + const std::string varName = hInfo["variable"].asString(); + // Integer values that define what particles need to be used + // for invariant mass combinations or helicity angles etc + int d1 = hInfo["d1"].asInt(); + int d2 = hInfo["d2"].asInt(); + + testInfo info( varName, d1, d2 ); + m_infoVect.push_back( info ); + + const std::string varTitle = hInfo["title"].asString(); + int nBins = hInfo["nbins"].asInt(); + double xmin = hInfo["xmin"].asDouble(); + double xmax = hInfo["xmax"].asDouble(); + + TString histName( varName.c_str() ); + if ( d1 != 0 ) { + histName += "_"; + histName += d1; + } + if ( d2 != 0 ) { + histName += "_"; + histName += d2; + } + TH1D* hist = new TH1D( histName.Data(), varTitle.c_str(), nBins, xmin, + xmax ); + hist->SetDirectory( outFile ); + m_histVect.push_back( hist ); + } m_nHistos = m_histVect.size(); - } -void testDecayModel::generateEvents( const std::string& decFile, const std::string& parentName, - int nEvents ) { - +void testDecayModel::generateEvents( const std::string& decFile, + const std::string& parentName, int nEvents ) +{ // Define the random number generator auto randomEngine = std::make_unique(); EvtAbsRadCorr* radCorrEngine = nullptr; std::list extraModels; #ifdef EVTGEN_EXTERNAL bool convertPythiaCodes( false ); bool useEvtGenRandom( true ); EvtExternalGenList genList( convertPythiaCodes, "", "gamma", useEvtGenRandom ); radCorrEngine = genList.getPhotosModel(); extraModels = genList.getListOfModels(); #endif - + EvtGen theGen( "../DECAY.DEC", "../evt.pdl", randomEngine.get(), - radCorrEngine, &extraModels ); + radCorrEngine, &extraModels ); theGen.readUDecay( decFile.c_str() ); // Generate the decays EvtId parId = EvtPDL::getId( parentName.c_str() ); for ( int i = 0; i < nEvents; i++ ) { + if ( i % 1000 == 0 ) { + std::cout << "Event " << nEvents - i << std::endl; + } + + // Initial 4-momentum and particle + EvtVector4R pInit( EvtPDL::getMass( parId ), 0.0, 0.0, 0.0 ); + EvtParticle* parent = EvtParticleFactory::particleFactory( parId, pInit ); + + theGen.generateDecay( parent ); - if ( i%1000 == 0 ) { std::cout << "Event " << nEvents - i << std::endl; } - - // Initial 4-momentum and particle - EvtVector4R pInit( EvtPDL::getMass( parId ), 0.0, 0.0, 0.0 ); - EvtParticle* parent = EvtParticleFactory::particleFactory( parId, pInit ); - - theGen.generateDecay( parent ); - - // Store information - for ( size_t i = 0; i < m_nHistos; i++ ) { - - testInfo info = m_infoVect[i]; - double value = getValue( parent, info.getName(), info.getd1(), info.getd2() ); - - TH1D* hist = m_histVect[i]; - if ( hist ) { - hist->Fill( value ); - } - - } - - //parent->printTree(); - - // Cleanup - parent->deleteTree(); - + // Store information + for ( size_t i = 0; i < m_nHistos; i++ ) { + testInfo info = m_infoVect[i]; + double value = getValue( parent, info.getName(), info.getd1(), + info.getd2() ); + + TH1D* hist = m_histVect[i]; + if ( hist ) { + hist->Fill( value ); + } + } + + //parent->printTree(); + + // Cleanup + parent->deleteTree(); } - } double testDecayModel::getValue( EvtParticle* parent, const std::string& varName, - const int d1, const int d2 ) const { - + const int d1, const int d2 ) const +{ double value( 0.0 ); - if ( !parent ) { return value; } + if ( !parent ) { + return value; + } // Daughters - const EvtParticle* par1 = d1 > 0 ? parent->getDaug( d1 - 1 ) : nullptr; - const EvtParticle* par2 = d2 > 0 ? parent->getDaug( d2 - 1 ) : nullptr; + const int NDaugMax( parent->getNDaug() ); + const EvtParticle* par1 = d1 > 0 && d1 <= NDaugMax ? parent->getDaug( d1 - 1 ) + : nullptr; + const EvtParticle* par2 = d2 > 0 && d2 <= NDaugMax ? parent->getDaug( d2 - 1 ) + : nullptr; // 4-momenta (in parent rest frame) const EvtVector4R p1 = par1 != nullptr ? par1->getP4() : EvtVector4R(); const EvtVector4R p2 = par2 != nullptr ? par2->getP4() : EvtVector4R(); if ( !varName.compare( "mass" ) ) { + // Invariant mass + if ( d2 != 0 ) { + // Invariant 4-mass combination of particles d1 and d2 + value = ( p1 + p2 ).mass(); + } else { + // Invariant mass of particle d1 only + value = p1.mass(); + } - // Invariant mass - if ( d2 != 0 ) { - // Invariant 4-mass combination of particles d1 and d2 - value = ( p1 + p2 ).mass(); - } else { - // Invariant mass of particle d1 only - value = p1.mass(); - } - } else if ( !varName.compare( "psumsq" ) ) { - - // Invariant momentum sum squared - value = ( p1 + p2 ).mass2(); + // Invariant momentum sum squared + value = ( p1 + p2 ).mass2(); } else if ( !varName.compare( "pdiffsq" ) ) { - - // Invariant momentum difference squared - value = ( p1 - p2 ).mass2(); + // Invariant momentum difference squared + value = ( p1 - p2 ).mass2(); } else if ( !varName.compare( "mtm" ) ) { + // Momentum of particle d1 + value = p1.d3mag(); - // Momentum of particle d1 - value = p1.d3mag(); - } else if ( !varName.compare( "coshel" ) ) { - - // Cosine of helicity angle - - // Resonance center-of-mass system (d1 and d2) - const EvtVector4R p12 = p1 + p2; - // Boost vector - const EvtVector4R boost( p12.get( 0 ), -p12.get( 1 ), -p12.get( 2 ), -p12.get( 3 ) ); - // Momentum of particle d1 in resonance frame, p1Res - const EvtVector4R p1Res = boostTo( p1, boost ); - // Cosine of angle between p1Res and momentum of resonance in parent frame - double p1ResMag = p1Res.d3mag(); - double p12Mag = p12.d3mag(); - if ( p1ResMag > 0.0 && p12Mag > 0.0 ) { - value = p1Res.dot( p12 ) / ( p1ResMag * p12Mag ); - } + // Cosine of helicity angle + + // Resonance center-of-mass system (d1 and d2) + const EvtVector4R p12 = p1 + p2; + // Boost vector + const EvtVector4R boost( p12.get( 0 ), -p12.get( 1 ), -p12.get( 2 ), + -p12.get( 3 ) ); + // Momentum of particle d1 in resonance frame, p1Res + const EvtVector4R p1Res = boostTo( p1, boost ); + // Cosine of angle between p1Res and momentum of resonance in parent frame + double p1ResMag = p1Res.d3mag(); + double p12Mag = p12.d3mag(); + if ( p1ResMag > 0.0 && p12Mag > 0.0 ) { + value = p1Res.dot( p12 ) / ( p1ResMag * p12Mag ); + } } else if ( !varName.compare( "prob" ) ) { - - // Decay probability - double* dProb = parent->decayProb(); + // Decay probability + double* dProb = parent->decayProb(); if ( dProb ) { - value = *dProb; + value = *dProb; } - } - + return value; - } -void testDecayModel::compareHistos( const std::string& refFileName ) const { - +void testDecayModel::compareHistos( const std::string& refFileName ) const +{ // Compare histograms with the same name, calculating the chi-squared TFile* refFile = TFile::Open( refFileName.c_str(), "read" ); if ( !refFile ) { - std::cout << "Could not open reference file " << refFileName << std::endl; - return; + std::cout << "Could not open reference file " << refFileName << std::endl; + return; } - - for ( auto hist : m_histVect ) { - - const std::string histName = hist->GetName(); - // Get equivalent reference histogram - const TH1* refHist = dynamic_cast( refFile->Get ( histName.c_str() ) ); - - if ( refHist ) { - - double chiSq( 0.0 ); - int nDof( 0 ); - int iGood( 0 ); - double pValue = refHist->Chi2TestX( hist, chiSq, nDof, iGood, "WW" ); - std::cout << "Histogram " << histName << " chiSq/nDof = " - << chiSq << "/" << nDof << ", pValue = " << pValue << std::endl; - } else { - - std::cout << "Could not find reference histogram " << histName << std::endl; - - } - + for ( auto hist : m_histVect ) { + const std::string histName = hist->GetName(); + // Get equivalent reference histogram + const TH1* refHist = dynamic_cast( + refFile->Get( histName.c_str() ) ); + + if ( refHist ) { + double chiSq( 0.0 ); + int nDof( 0 ); + int iGood( 0 ); + double pValue = refHist->Chi2TestX( hist, chiSq, nDof, iGood, "WW" ); + std::cout << "Histogram " << histName << " chiSq/nDof = " << chiSq + << "/" << nDof << ", pValue = " << pValue << std::endl; + + } else { + std::cout << "Could not find reference histogram " << histName + << std::endl; + } } refFile->Close(); - } -int main( int argc, char* argv[] ) { - +int main( int argc, char* argv[] ) +{ if ( argc != 2 ) { - std::cout << "Expecting one argument: json input file" << std::endl; - return -1; + std::cout << "Expecting one argument: json input file" << std::endl; + return -1; } const std::string jsonFile = argv[1]; testDecayModel test( jsonFile ); test.run(); return 0; - } diff --git a/test/testDecayModel.hh b/test/testDecayModel.hh index 10b2196..5fddd8c 100644 --- a/test/testDecayModel.hh +++ b/test/testDecayModel.hh @@ -1,87 +1,80 @@ /*********************************************************************** * Copyright 1998-2020 CERN for the benefit of the EvtGen authors * * * * This file is part of EvtGen. * * * * EvtGen is free software: you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation, either version 3 of the License, or * * (at your option) any later version. * * * * EvtGen is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with EvtGen. If not, see . * ***********************************************************************/ #ifndef TEST_DECAY_MODEL_HH #define TEST_DECAY_MODEL_HH -#include "json/json.h" - #include "TFile.h" #include "TH1D.h" +#include "json/json.h" + #include #include class EvtParticle; class testInfo { - public: testInfo( const std::string& name, int d1, int d2 ) : - m_name( name ), - m_d1( d1 ), - m_d2( d2 ) - {;} - + m_name( name ), m_d1( d1 ), m_d2( d2 ) + { + ; + } + const std::string getName() const { return m_name; } int getd1() const { return m_d1; } int getd2() const { return m_d2; } - + private: std::string m_name; int m_d1; int m_d2; - }; class testDecayModel { - public: - testDecayModel( const std::string& jsonFile ); - void run(); + void run(); private: - std::string createDecFile( const std::string& parent, - const std::vector& daugNames, - const std::string& model, - const std::vector& parameters ) const; + const std::vector& daugNames, + const std::string& model, + const std::vector& parameters ) const; void defineHistos( Json::Value& config, TFile* theFile ); - void generateEvents( const std::string& decFile, const std::string& parentName, - int nEvents ); + void generateEvents( const std::string& decFile, + const std::string& parentName, int nEvents ); double getValue( EvtParticle* rootPart, const std::string& varName, - const int d1, const int d2 ) const; - + const int d1, const int d2 ) const; + void compareHistos( const std::string& refFileName ) const; - + std::string m_jsonFile; std::vector m_infoVect; std::vector m_histVect; size_t m_nHistos; - - }; #endif