diff --git a/EvtGenExternal/EvtPythiaEngine.hh b/EvtGenExternal/EvtPythiaEngine.hh index 47640ba..d0e0dd8 100644 --- a/EvtGenExternal/EvtPythiaEngine.hh +++ b/EvtGenExternal/EvtPythiaEngine.hh @@ -1,98 +1,114 @@ /*********************************************************************** * 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 #ifndef EVTPYTHIAENGINE_HH #define EVTPYTHIAENGINE_HH #include "EvtGenBase/EvtDecayAmp.hh" #include "EvtGenBase/EvtDecayBase.hh" #include "EvtGenBase/EvtId.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtVector4R.hh" #include "EvtGenModels/EvtAbsExternalGen.hh" #include "EvtGenExternal/EvtPythiaRandom.hh" #include "Pythia8/ParticleData.h" #include "Pythia8/Pythia.h" #include #include #include #include // Description: Interface to the Pythia 8 external generator class EvtPythiaEngine : public EvtAbsExternalGen { public: EvtPythiaEngine( std::string xmlDir = "./xmldoc", bool convertPhysCodes = false, bool useEvtGenRandom = true ); virtual ~EvtPythiaEngine(); bool doDecay( EvtParticle* theMother ) override; void initialise() override; + void initialiseHmeAmplitudes( int nA, int nB, int nC ); + double hmeAmplitudes( int hmeMode, EvtParticle* mom = nullptr, EvtDecayAmp* amp = nullptr ); protected: private: void updateParticleLists(); void updatePhysicsParameters(); void createPythiaParticle( EvtId& particleId, int PDGCode ); bool validPDGCode( int PDGCode ); void updatePythiaDecayTable( EvtId& particleId, int aliasInt, int PDGCode ); void storeDaughterInfo( EvtParticle* theParticle, int startInt ); void clearDaughterVectors(); void clearPythiaModeMap(); void createDaughterEvtParticles( EvtParticle* theParent ); int getModeInt( EvtDecayBase* decayModel ); std::unique_ptr _genericPythiaGen; std::unique_ptr _aliasPythiaGen; Pythia8::Pythia* _thePythiaGenerator; std::vector _daugPDGVector; std::vector _daugP4Vector; typedef std::map> PythiaModeMap; PythiaModeMap _pythiaModeMap; - bool _convertPhysCodes, _initialised, _useEvtGenRandom; + bool _convertPhysCodes, _initialised, _useEvtGenRandom, + _initialisedHmeAmplitudes; std::unique_ptr _evtgenRandom; std::map _addedPDGCodes; std::map> _hmes; + + //Spin states available for particle A, B, and C. (for HmeAmplitudes) + int _nA, _nB, _nC; + + //temporary array for amplitudes + EvtComplexPtrPtrPtr _amp, _amp1, _amp3; + + //Rotation matrices + EvtComplexPtrPtr _RA, _RB, _RC; + void setUpRotationMatrices( std::vector helParts ); + void applyRotationMatrices(); + + EvtSpinDensity rotateToHelicityBasis( EvtParticle* mom ) const; }; #endif #endif diff --git a/src/EvtGenExternal/EvtHmePythia.cpp b/src/EvtGenExternal/EvtHmePythia.cpp index 002d53c..2888a82 100644 --- a/src/EvtGenExternal/EvtHmePythia.cpp +++ b/src/EvtGenExternal/EvtHmePythia.cpp @@ -1,91 +1,123 @@ /*********************************************************************** * 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/EvtPDL.hh" #include "EvtGenBase/EvtParticle.hh" +#include "EvtGenBase/EvtSpinType.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 ); + // 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 ); + + // Initialize dimensions, amplitudes and rotation matrices + int ih( 0 ); + std::vector nhs( 3, 1 ); + + int momSpinStates = EvtSpinType::getSpinStates( + EvtPDL::getSpinType( getParentId() ) ); + + if ( momSpinStates > 1 ) { + nhs[ih] = momSpinStates; + ++ih; + } + + for ( int iDaug = 0; iDaug < getNDaug(); ++iDaug ) { + int iDaugSpinState = EvtSpinType::getSpinStates( + EvtPDL::getSpinType( getDaug( iDaug ) ) ); + if ( iDaugSpinState > 1 ) { + if ( ih < 3 ) { + nhs[ih] = iDaugSpinState; + ++ih; + } else { + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "More than three particles with multiple spins." + << std::endl; + } + } + } + + _pythiaEngine->initialiseHmeAmplitudes( nhs[0], nhs[1], nhs[2] ); + 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 ); + double new_prob_max = _pythiaEngine->hmeAmplitudes( _hmeMode, p, this ); + // setProbMax(new_prob_max); } } diff --git a/src/EvtGenExternal/EvtPythiaEngine.cpp b/src/EvtGenExternal/EvtPythiaEngine.cpp index 60261b0..ad375d3 100644 --- a/src/EvtGenExternal/EvtPythiaEngine.cpp +++ b/src/EvtGenExternal/EvtPythiaEngine.cpp @@ -1,972 +1,1651 @@ /*********************************************************************** * 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/EvtDiracParticle.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; + + _initialisedHmeAmplitudes = false; } EvtPythiaEngine::~EvtPythiaEngine() { _thePythiaGenerator = nullptr; this->clearDaughterVectors(); this->clearPythiaModeMap(); + + // Clear up memory used for Hme amplitudes and rotation matrices + + if ( _initialisedHmeAmplitudes ) { + int ia, ib, ic; + for ( ia = 0; ia < _nA; ia++ ) { + delete[] _RA[ia]; + } + delete[] _RA; + + for ( ib = 0; ib < _nB; ib++ ) { + delete[] _RB[ib]; + } + delete[] _RB; + + for ( ic = 0; ic < _nC; ic++ ) { + delete[] _RC[ic]; + } + delete[] _RC; + + for ( ia = 0; ia < _nA; ia++ ) { + for ( ib = 0; ib < _nB; ib++ ) { + delete[] _amp[ia][ib]; + delete[] _amp1[ia][ib]; + delete[] _amp3[ia][ib]; + } + delete[] _amp[ia]; + delete[] _amp1[ia]; + delete[] _amp3[ia]; + } + + delete[] _amp; + delete[] _amp1; + delete[] _amp3; + } } 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]; + + // std::cout << "hmeMode = " << hmeMode << std::endl; + if ( hme == nullptr ) { return 0.0; } + // Hme Amplitudes need to be initialized + if ( _initialisedHmeAmplitudes == false ) { + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << " The hmeAmplitudes are not initialized. Use initialiseHmeAmplitudes." + << endl; + } + + std::cout << " " << std::endl; + std::cout << " ----------- New event --------- " << std::endl; + // Return the maximum amplitude if no particle is given std::vector prts; if ( mom == nullptr ) { // 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 ); + // Factor determined empirically by running. Otherwise prob becomes larger than probmax. + std::cout << "decayWeightMax: " << hme->decayWeightMax( prts ) + << std::endl; + return 5 * hme->decayWeightMax( prts ); // hme->decayWeightMax( prts ); } + // Counter for particles with more than one helicity state int ih( 0 ); - std::vector ihs( 3, -1 ), nhs( 3, 0 ), hs( mom->getNDaug() + 1, 0 ); + + // Vectors containing the number/index of helicity states + // We support only 3 particles with more than one helicity state. + // Initial index of helicity states ihs (for loops), + // indices for all helicity states hs. + std::vector ihs( 3, -1 ), hs( mom->getNDaug() + 1, 0 ); + + // Vector containing pointers to particles with more than one helicity state. + std::vector helParts( 3, nullptr ); // 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; + // The mom is always a tau so we will go always inside this loop. if ( mom->getSpinStates() > 1 ) { ihs[ih] = 0; - nhs[ih] = mom->getSpinStates(); + helParts[ih] = mom; ++ih; } + std::cout << "Momenta in lab frame " << std::endl; + std::cout << "Tau: p = " << p.get( 0 ) << ", " << p.get( 1 ) << ", " + << p.get( 2 ) << ", " << p.get( 3 ) << std::endl; + // 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; + + std::cout << "Daughter " << iDtr << ": PDG code " << dtr->getPDGId() + << ": p = " << p.get( 0 ) << ", " << p.get( 1 ) << ", " + << p.get( 2 ) << ", " << p.get( 3 ) << std::endl; + if ( dtr->getSpinStates() > 1 ) { if ( ih < 3 ) { - ihs[ih] = 0; - nhs[ih] = dtr->getSpinStates(); + ihs[ih] = iDtr + 1; + helParts[ih] = dtr; ++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; } + // std::cout << "Before" << std::endl; + // std::cout << "ihs: " << ihs[0] << ", " << ihs[1] << ", " << ihs[2] << std::endl; + // std::cout << "hs: " << hs[0] << ", " << hs[1] << ", " << hs[2] << std::endl; + // std::cout << "hs[ihs] : " << hs[ihs[0]] << ", " << hs[ihs[1]] << ", " << hs[ihs[2]] << std::endl; + // std::cout << " n: " << _nA << ", " << _nB << ", " << _nC << std::endl; + // 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 ); + hme->initWaves( prts ); + + // std::cout << "Begin " << std::endl; + // std::cout << "Test HME output i: " << std::endl; + int ia, ib, ic; + + // for ( ia = 0; ia < _nA+1; ia++ ) { + // for ( ib = 0; ib < _nB+1; ib++ ) { + // for ( ic = 0; ic < _nC; ic++ ) { + // std::cout << "Test HME output i: " << ia << ", " << ic << ", " << ib << std::endl; + // // std::cout << "ihs: " << ihs[0] << ", " << ihs[1] << ", " << ihs[2] << std::endl; + // // std::vector myHS = {ia, 0, ib, 1}; + // std::vector myHS = {ia, 1, 0, ib}; + // // std::vector myHS = {ia, 1, ib}; + // // std::vector myHS = {ia, 0, ib}; + // std::cout << hme->calculateME(myHS) << std::endl; + // // _amp[ia][ib][ic] = hme->calculateME(myHS); + // } + // } + // } + // std::cout << " dims: " << _nA << ", " << _nB << ", " << _nC << "amp[0] = " << _amp[0] << std::endl; + + // Fill the temporary array with the amplitudes + double prob1 = 0.0; + + // if ( ihs[0] < 0 ) { + // // Only for mother with 1 spin state. So never applies for tau decays. + // _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] = hme->calculateME( hs ); + // prob1 += real( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + // conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ); + // } else { + // if ( ihs[1] < 0 ) { + // // Most simple case: all daughters (except neutrino) are scalar. + // for ( hs[ihs[0]] = 0; hs[ihs[0]] < _nA; ++hs[ihs[0]] ) { + // _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] = hme->calculateME( hs ); + // prob1 += real( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + // conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ); + // // std::cout << "hs:" << hs[ihs[0]] << " " << hs[ihs[1]] << " " << hs[ihs[2]] + // // << "prob1: " << prob1 << ", amp: " << real( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + // // conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ) << std::endl; + // // std::cout << _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] << " vs. " << hme->calculateME( hs ) << std::endl; + // // std::cout << "hs[ihs]: " << hs[ihs[0]] << ", " << hs[ihs[1]] << ", " << hs[ihs[2]] << + // // ": " << _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] << std::endl; + // // std::cout << "ihs: " << ihs[0] << ", " << ihs[1] << ", " << ihs[2] << std::endl; + // // std::cout << hs[ihs[0]] << hs[ihs[1]] << hs[ihs[2]] << " vs. " << hme->calculateME( hs ) << std::endl; + // } + // } else { + // if ( ihs[2] < 0 ) { + // for ( hs[ihs[1]] = 0; hs[ihs[1]] < _nA; ++hs[ihs[1]] ) { + // for ( hs[ihs[0]] = 0; hs[ihs[0]] < _nB; ++hs[ihs[0]] ) { + // // Next case: one daughter is non-scalar. + // if ( _nB <= _nA ) { // First daughter has two spin states. + // _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] = + // hme->calculateME( hs ); + // prob1 += real( + // _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + // conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ); + // } else { // First daughter has more than two spin states. + // // When the first daughter has more spin states than the mother (example tau -> rhonu), + // // somehow Pythia permutes the order of the indices + // _amp[hs[ihs[1]]][hs[ihs[0]]][hs[ihs[2]]] = + // hme->calculateME( hs ); + // // std::cout << hs[ihs[1]] << hs[ihs[0]] << hs[ihs[2]] << " vs. " << hme->calculateME( hs ) << std::endl; + // prob1 += real( + // _amp[hs[ihs[1]]][hs[ihs[0]]][hs[ihs[2]]] * + // conj( _amp[hs[ihs[1]]][hs[ihs[0]]][hs[ihs[2]]] ) ); + // } + // } + // } + // } else { // Not tested yet. Requires tau decaying into final states with two non-scalar particles + // for ( hs[ihs[1]] = 0; hs[ihs[1]] < _nA; ++hs[ihs[1]] ) { + // for ( hs[ihs[0]] = 0; hs[ihs[0]] < _nB; ++hs[ihs[0]] ) { + // for ( hs[ihs[2]] = 0; hs[ihs[2]] < _nC; ++hs[ihs[2]] ) { + // _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] = + // hme->calculateME( hs ); + // prob1 += real( + // _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + // conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ); + // } + // } + // } + // } + // } + // } + + std::cout << "HME elements as obtained from pythia: " << std::endl; if ( ihs[0] < 0 ) { - amp->vertex( hme->calculateME( hs ) ); + // Only for mother with 1 spin state. So never applies for tau decays. + _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] = hme->calculateME( hs ); + prob1 += real( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ); } else { - for ( hs[ihs[0]] = 0; hs[ihs[0]] < nhs[0]; ++hs[ihs[0]] ) { + for ( hs[ihs[0]] = 0; hs[ihs[0]] < _nA; ++hs[ihs[0]] ) { if ( ihs[1] < 0 ) { - amp->vertex( hs[ihs[0]], hme->calculateME( hs ) ); + // Most simple case: all daughters (except neutrino) are scalar. + _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] = hme->calculateME( hs ); + prob1 += real( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ); + // std::cout << "hs:" << hs[ihs[0]] << " " << hs[ihs[1]] << " " << hs[ihs[2]] + // << "prob1: " << prob1 << ", amp: " << real( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + // conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ) << std::endl; + // std::cout << _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] << " vs. " << hme->calculateME( hs ) << std::endl; + // std::cout << "hs[ihs]: " << hs[ihs[0]] << ", " << hs[ihs[1]] << ", " << hs[ihs[2]] << + // ": " << _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] << std::endl; + // std::cout << "ihs: " << ihs[0] << ", " << ihs[1] << ", " << ihs[2] << std::endl; + std::cout << "State index : " << hs[ihs[0]] << hs[ihs[1]] + << hs[ihs[2]] << " : " << hme->calculateME( hs ) + << std::endl; } else { - for ( hs[ihs[1]] = 0; hs[ihs[1]] < nhs[1]; ++hs[ihs[1]] ) { + for ( hs[ihs[1]] = 0; hs[ihs[1]] < _nB; ++hs[ihs[0]] ) { if ( ihs[2] < 0 ) { - amp->vertex( hs[ihs[0]], hs[ihs[1]], - hme->calculateME( hs ) ); + // Next case: one daughter is non-scalar. + _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] = + hme->calculateME( hs ); + prob1 += real( + _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ); + } else { // Not tested yet. Requires tau decaying into final states with two non-scalar particles + for ( hs[ihs[2]] = 0; hs[ihs[2]] < _nC; ++hs[ihs[2]] ) { + _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] = + hme->calculateME( hs ); + prob1 += real( + _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ); + } + } + } + } + } + } + + // Here we apply a rotation from the Helicity basis (Pythia) to the spin basis (EvtGen) + // in the same way as it is done in EvtEvalHelAmp.cpp + setUpRotationMatrices( helParts ); + + applyRotationMatrices(); + + // Here we propagate the rotated amplitudes + double prob2 = 0.0; + + std::cout << "Propagated HME elements: " << std::endl; + + //int ia, ib, ic; + for ( ia = 0; ia < _nA; ia++ ) { + for ( ib = 0; ib < _nB; ib++ ) { + for ( ic = 0; ic < _nC; ic++ ) { + prob2 += real( _amp[ia][ib][ic] * conj( _amp[ia][ib][ic] ) ); + // std::cout << "i:" << ia << " " << ib << " " << ic + // << "prob2: " << prob2 << ", amp: " << real( _amp[ia][ib][ic] * conj( _amp[ia][ib][ic] ) ) << std::endl; + if ( _nA == 1 ) { + if ( _nB == 1 ) { + if ( _nC == 1 ) { + amp->vertex( _amp[ia][ib][ic] ); + } else { + amp->vertex( ic, _amp[ia][ib][ic] ); + } } 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 ) ); + if ( _nC == 1 ) { + amp->vertex( ib, _amp[ia][ib][ic] ); + } else { + amp->vertex( ib, ic, _amp[ia][ib][ic] ); } - } // End third helicity + } + } else { + if ( _nB == 1 ) { + if ( _nC == 1 ) { + amp->vertex( ia, _amp[ia][ib][ic] ); + std::cout << "State index " << ia << ", " << ib + << " : " << _amp[ia][ib][ic] << std::endl; + } else { + amp->vertex( ia, ic, _amp[ia][ib][ic] ); + } + } else { + if ( _nC == 1 ) { + amp->vertex( ia, ib, _amp[ia][ib][ic] ); + } else { + amp->vertex( ia, ib, ic, _amp[ia][ib][ic] ); + } + } } - } // End second helicity + } } - } // End first helicity + } + + // Modified original code to pass rotated elements + // if ( fabs( prob1 - prob2 ) > 0.000001 * prob1 ) { + // EvtGenReport( EVTGEN_INFO, "EvtGen" ) + // << "prob1,prob2:" << prob1 << " " << prob2 << endl; + // ::abort(); + // } + + // if ( ihs[0] < 0 ) { + // // Only for mother with 1 spin state. So never applies for tau decays. + // amp->vertex(_amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]]); + // } else { + // for ( hs[ihs[0]] = 0; hs[ihs[0]] < _nA; ++hs[ihs[0]] ) { + // if ( ihs[1] < 0 ) { + // amp->vertex( hs[ihs[0]], _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ); + // // std::cout << "Case 2. " << std::endl; + // // std::cout << _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] << " vs. " << hme->calculateME( hs ) << std::endl; + // // std::cout << "hs[ihs]: " << hs[ihs[0]] << ", " << hs[ihs[1]] << ", " << hs[ihs[2]] << ": " << _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] << std::endl; + // std::cout << "hs[ihs]: " << hs[ihs[0]] << ", " << hs[ihs[1]] << ", " << hs[ihs[2]] << ":, amp: " << + // real( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + // conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ) << std::endl; + // // std::cout << "ihs: " << ihs[0] << ", " << ihs[1] << ", " << ihs[2] << std::endl; + // } else { + // for ( hs[ihs[1]] = 0; hs[ihs[1]] < _nB; ++hs[ihs[1]] ) { + // if ( ihs[2] < 0 ) { + // amp->vertex( hs[ihs[0]], hs[ihs[1]],_amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ); + // // std::cout << "Case 3. " << std::endl; + // } else { + // for ( hs[ihs[2]] = 0; hs[ihs[2]] < _nC; ++hs[ihs[2]] ) { + // amp->vertex( hs[ihs[0]], hs[ihs[1]], hs[ihs[2]], _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ); + // // std::cout << "Case 4. " << std::endl; + // } + // } // End third helicity + // } + // } // End second helicity + // } + // } // End first helicity + + // // Modified original code to pass the right elements + // + // if ( ihs[0] < 0 ) { + // // Only for mother with 1 spin state. So never applies for tau decays. + // amp->vertex(hme->calculateME( hs )); + // } else { + // if ( ihs[1] < 0 ) { + // // Most simple case: all daughters (except neutrino) are scalar. + // for ( hs[ihs[0]] = 0; hs[ihs[0]] < _nA; ++hs[ihs[0]] ) { + // amp->vertex(hs[ihs[0]], hme->calculateME( hs )); + // // std::cout << "hs:" << hs[ihs[0]] << " " << hs[ihs[1]] << " " << hs[ihs[2]] + // // << "prob1: " << prob1 << ", amp: " << real( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] * + // // conj( _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] ) ) << std::endl; + // // std::cout << _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] << " vs. " << hme->calculateME( hs ) << std::endl; + // // std::cout << "hs[ihs]: " << hs[ihs[0]] << ", " << hs[ihs[1]] << ", " << hs[ihs[2]] << + // // ": " << _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] << std::endl; + // // std::cout << "ihs: " << ihs[0] << ", " << ihs[1] << ", " << ihs[2] << std::endl; + // hs[ihs[0]]++; + // } + // } else { + // if ( ihs[2] < 0 ) { + // for ( hs[ihs[1]] = 0; hs[ihs[1]] < _nA; ++hs[ihs[1]] ) { + // for ( hs[ihs[0]] = 0; hs[ihs[0]] < _nB; ++hs[ihs[0]] ) { + // // Next case: one daughter is non-scalar. + // // if ( _nB <= _nA ) { // First daughter has two spin states. + // // amp->vertex(hs[ihs[0]], hs[ihs[1]], hme->calculateME( hs )); + // // } else { // First daughter has more than two spin states. + // // When the first daughter has more spin states than the mother (example tau -> rhonu), + // // somehow Pythia permutes the order of the indices + // amp->vertex(hs[ihs[1]], hs[ihs[0]], hme->calculateME( hs )); + // hs[ihs[0]]++; + // // } + // } + // } + // } else { // Not tested yet. Requires tau decaying into final states with two non-scalar particles + // for ( hs[ihs[1]] = 0; hs[ihs[1]] < _nA; ++hs[ihs[1]] ) { + // for ( hs[ihs[0]] = 0; hs[ihs[0]] < _nB; ++hs[ihs[0]] ) { + // for ( hs[ihs[2]] = 0; hs[ihs[2]] < _nC; ++hs[ihs[2]] ) { + // amp->vertex( hs[ihs[0]], hs[ihs[1]], hs[ihs[2]], hme->calculateME( hs )); + // } + // } + // } + // } + // } + // } + + // //#### Original code + // if ( ihs[0] < 0 ) { + // // Only for mother with 1 spin state. So never applies for tau decays. + // amp->vertex(hme->calculateME( hs )); + // } else { + // for ( hs[ihs[0]] = 0; hs[ihs[0]] < _nA; ++hs[ihs[0]] ) { + // if ( ihs[1] < 0 ) { + // amp->vertex( hs[ihs[0]], hme->calculateME( hs )); + // // std::cout << "Case 2. " << std::endl; + // // std::cout << _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] << " vs. " << hme->calculateME( hs ) << std::endl; + // // std::cout << "hs[ihs]: " << hs[ihs[0]] << ", " << hs[ihs[1]] << ", " << hs[ihs[2]] << ": " << _amp[hs[ihs[0]]][hs[ihs[1]]][hs[ihs[2]]] << std::endl; + // // std::cout << "ihs: " << ihs[0] << ", " << ihs[1] << ", " << ihs[2] << std::endl; + // } else { + // for ( hs[ihs[1]] = 0; hs[ihs[1]] < _nB; ++hs[ihs[1]] ) { + // if ( ihs[2] < 0 ) { + // amp->vertex( hs[ihs[0]], hs[ihs[1]], hme->calculateME( hs )); + // // std::cout << "Case 3. " << std::endl; + // } else { + // for ( hs[ihs[2]] = 0; hs[ihs[2]] < _nC; ++hs[ihs[2]] ) { + // amp->vertex( hs[ihs[0]], hs[ihs[1]], hs[ihs[2]], hme->calculateME( hs )); + // // std::cout << "Case 4. " << std::endl; + // } + // } // End third helicity + // } + // } // End second helicity + // } + // } // End first helicity + + // std::cout << "prob1: " << prob1 << std::endl; + // std::cout << "prob2: " << prob2 << std::endl; + // std::cout << "Decay weight: " << hme->decayWeight( prts ) << std::endl; + // std::cout << "Decay weight max: " << hme->decayWeightMax( prts ) << std::endl; // Return the decay weight return hme->decayWeight( prts ); + // return hme->decayWeightMax( prts ); + // return prob1; +} + +/* Initialize HME amplitudes as private variables to rotate them. */ +void EvtPythiaEngine::initialiseHmeAmplitudes( int nA, int nB, int nC ) +{ + if ( _initialised == false ) { + this->initialise(); + } + + if ( _initialisedHmeAmplitudes ) { + return; + } + + //Find out how many states each particle has + _nA = nA; + _nB = nB; + _nC = nC; + + // Initialize temporary array for the amplitudes + int ia, ib, ic; + + _amp = new EvtComplexPtrPtr[_nA]; + _amp1 = new EvtComplexPtrPtr[_nA]; + _amp3 = new EvtComplexPtrPtr[_nA]; + + for ( ia = 0; ia < _nA; ia++ ) { + _amp[ia] = new EvtComplexPtr[_nB]; + _amp1[ia] = new EvtComplexPtr[_nB]; + _amp3[ia] = new EvtComplexPtr[_nB]; + for ( ib = 0; ib < _nB; ib++ ) { + _amp[ia][ib] = new EvtComplex[_nC]; + _amp1[ia][ib] = new EvtComplex[_nC]; + _amp3[ia][ib] = new EvtComplex[_nC]; + } + } + + for ( ia = 0; ia < _nA; ia++ ) { + for ( ib = 0; ib < _nB; ib++ ) { + for ( ic = 0; ic < _nC; ic++ ) { + _amp[ia][ib][ic] = 0.0; + } + } + } + + _RA = new EvtComplexPtr[_nA]; + for ( int ia = 0; ia < _nA; ia++ ) { + _RA[ia] = new EvtComplex[_nA]; + } + _RB = new EvtComplexPtr[_nB]; + for ( int ib = 0; ib < _nB; ib++ ) { + _RB[ib] = new EvtComplex[_nB]; + } + _RC = new EvtComplexPtr[_nC]; + for ( int ic = 0; ic < _nC; ic++ ) { + _RC[ic] = new EvtComplex[_nC]; + } + + _initialisedHmeAmplitudes = true; +} + +// Copy pasted from src/EvtGenBase/EvtDiracParticle.cpp just to do local checks. +EvtSpinDensity EvtPythiaEngine::rotateToHelicityBasis( EvtParticle* mom ) const +{ + EvtDiracSpinor spplus; + EvtDiracSpinor spminus; + + double sqmt2 = sqrt( 2. * ( mom->getP4().mass() ) ); + + if ( EvtPDL::getStdHep( mom->getId() ) > 0 ) { + spplus.set( 1.0, 0.0, 0.0, 0.0 ); + spminus.set( 0.0, 1.0, 0.0, 0.0 ); + } else { + spplus.set( 0.0, 0.0, 1.0, 0.0 ); + spminus.set( 0.0, 0.0, 0.0, 1.0 ); + } + + EvtSpinDensity R; + R.setDim( 2 ); + + for ( int i = 0; i < 2; i++ ) { + R.set( 0, i, ( spplus * mom->spParent( i ) ) / sqmt2 ); + R.set( 1, i, ( spminus * mom->spParent( i ) ) / sqmt2 ); + } + + return R; +} + +/* Setup rotation matrices based on the helicity states of the particles + * and the momentum of the daughters. Modified from + * src/EvtGenBase/EvtEvalHelAmp.cpp */ + +void EvtPythiaEngine::setUpRotationMatrices( std::vector helParts ) +{ + int JA2 = helParts[0] != nullptr + ? EvtSpinType::getSpin2( helParts[0]->getSpinType() ) + : 0; + int JB2 = helParts[1] != nullptr + ? EvtSpinType::getSpin2( helParts[1]->getSpinType() ) + : 0; + int JC2 = helParts[2] != nullptr + ? EvtSpinType::getSpin2( helParts[2]->getSpinType() ) + : 0; + + double theta( 0 ), phi( 0 ); + + // Tried to transform the helicity aplitudes of the tau itself, but not needed + // because spin quantization axis should be congruent between EvtGen and Pythia. + // if ( helParts[0] != nullptr ) { + // EvtVector4R pMom = helParts[0]->getP4(); + // theta = acos( pMom.get( 3 ) / pMom.d3mag() ); + // phi = atan2( pMom.get( 2 ), pMom.get( 1 ) ); + // } + + // First find theta and phi of the first daughter in the mother's frame + if ( helParts[1] != nullptr ) { + EvtVector4R pDaug = helParts[1]->getP4(); + theta = acos( pDaug.get( 3 ) / pDaug.d3mag() ); + phi = atan2( pDaug.get( 2 ), pDaug.get( 1 ) ); + } + + switch ( JA2 ) { + case 0: + _RA[0][0] = 1.0; + break; + case 1: + case 2: + case 3: + case 4: + case 5: + case 6: + case 7: + case 8: + + { + // EvtSpinDensity R = helParts[0] ->rotateToHelicityBasis(); + // EvtConst::pi, EvtConst::pi, EvtConst::pi ); + // EvtSpinDensity R = helParts[0] -> rotateToHelicityBasis(0, EvtConst::pi/2, 0); + + // If no angles are given as arguments, the rotation matrix is unity + EvtSpinDensity R = helParts[0]->rotateToHelicityBasis(); + + int i, j, n; + + n = R.getDim(); + + assert( n == _nA ); + + std::cout << "Rotation matrix: tau ->rotateToHelicityBasis() = " + << std::endl; + + // std::cout << phi << theta << -phi << std::endl; + + // std::cout << "Evt ID: " << EvtPDL::getStdHep( helParts[0] -> getId() ) << std::endl; + + for ( i = 0; i < n; i++ ) { + for ( j = 0; j < n; j++ ) { + // for ( j = 1; j < n; j++ ) { + _RA[i][j] = R.get( i, j ); + // if (i==0) _RA[i][j] = 0; + std::cout << "R(" << i << ", " << j << ") = " << _RA[i][j] + << std::endl; + } + } + + EvtSpinDensity F = helParts[0]->getSpinDensityForward(); + + std::cout << "Spin density matrix: tau->getSpinDensityForward() = " + << std::endl; + + n = F.getDim(); + for ( i = 0; i < n; i++ ) { + for ( j = 0; j < n; j++ ) { + std::cout << "i: " << i << " j: " << j << " " + << F.get( i, j ) << std::endl; + } + } + + // Print out Spin density matrix rotated using angles phi, -EvtConst::pi, -phi (just for technical check) + // helParts[0] -> setSpinDensityForwardHelicityBasis(helParts[0]->getSpinDensityForward(), phi, -EvtConst::pi, -phi); + // F = helParts[0]->getSpinDensityForward(); + // // // EvtSpinDensity F = mom->getSpinDensityBackward(); + + // // n = F.getDim(); + + // std::cout << "tau -> setSpinDensityForwardHelicityBasis(tau->getSpinDensityForward()) " << std::endl; + // for ( i = 0; i < n; i++ ) { + // for ( j = 0; j < n; j++ ) { + // std::cout << "i: " << i << " j: " << j << " " << F.get(i,j) + // << std::endl; + // } + // } + // std::cout << " n = " << n << std::endl; + } + + break; + + default: + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "Spin2(JA2)=" << JA2 << " not supported!" << endl; + ::abort(); + } + + switch ( JB2 ) { + case 0: + _RB[0][0] = 1.0; + break; + case 1: + case 2: + case 3: + case 4: + case 5: + case 6: + case 7: + case 8: + + { + int i, j, n; + + EvtSpinDensity R = helParts[1]->rotateToHelicityBasis( phi, theta, + -phi ); + + n = R.getDim(); + + assert( n == _nB ); + + for ( i = 0; i < n; i++ ) { + for ( j = 0; j < n; j++ ) { + _RB[i][j] = conj( R.get( i, j ) ); + } + } + } + + break; + + default: + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "Spin2(JB2)=" << JB2 << " not supported!" << endl; + ::abort(); + } + + switch ( JC2 ) { + case 0: + _RC[0][0] = 1.0; + break; + case 1: + case 2: + case 3: + case 4: + case 5: + case 6: + case 7: + case 8: + + { + int i, j, n; + + // Not sure if this rotation is correct because this is taken from Helamp + // which covers only 2-body decays. + EvtSpinDensity R = helParts[2]->rotateToHelicityBasis( + phi, EvtConst::pi + theta, phi - EvtConst::pi ); + + n = R.getDim(); + + assert( n == _nC ); + + for ( i = 0; i < n; i++ ) { + for ( j = 0; j < n; j++ ) { + _RC[i][j] = conj( R.get( i, j ) ); + } + } + + } + + break; + + default: + EvtGenReport( EVTGEN_ERROR, "EvtGen" ) + << "Spin2(JC2)=" << JC2 << " not supported!" << endl; + ::abort(); + } +} + +/* Function to apply the rotation on the matrix of amplitudes. + * Modified from src/EvtGenBase/EvtEvalHelAmp.cpp */ +void EvtPythiaEngine::applyRotationMatrices() +{ + int ia, ib, ic, i; + + EvtComplex temp; + + for ( ia = 0; ia < _nA; ia++ ) { + for ( ib = 0; ib < _nB; ib++ ) { + for ( ic = 0; ic < _nC; ic++ ) { + temp = 0; + for ( i = 0; i < _nC; i++ ) { + temp += _RC[i][ic] * _amp[ia][ib][i]; + } + _amp1[ia][ib][ic] = temp; + } + } + } + + for ( ia = 0; ia < _nA; ia++ ) { + for ( ic = 0; ic < _nC; ic++ ) { + for ( ib = 0; ib < _nB; ib++ ) { + temp = 0; + for ( i = 0; i < _nB; i++ ) { + temp += _RB[i][ib] * _amp1[ia][i][ic]; + } + _amp3[ia][ib][ic] = temp; + } + } + } + + for ( ib = 0; ib < _nB; ib++ ) { + for ( ic = 0; ic < _nC; ic++ ) { + for ( ia = 0; ia < _nA; ia++ ) { + temp = 0; + for ( i = 0; i < _nA; i++ ) { + temp += _RA[i][ia] * _amp3[i][ib][ic]; + } + _amp[ia][ib][ic] = temp; + } + } + } } #endif diff --git a/test/jsonFiles/HmePythiaTau_BtoTaunu_pi.json b/test/jsonFiles/HmePythiaTau_BtoTaunu_pi.json new file mode 100644 index 0000000..57f471d --- /dev/null +++ b/test/jsonFiles/HmePythiaTau_BtoTaunu_pi.json @@ -0,0 +1,42 @@ +{ + "parent" : "B+", + "daughters" : ["tau+", "nu_tau"], + "grand_daughters" : [["anti-nu_tau", "pi+"], []], + "models" : ["SLN", "HMEPYTHIA", ""], + "parameters" : [[], [1521], []], + "extras" : ["Define dm_incohMix_B0 0", "noPhotos"], + "events" : 100000, + "histograms" : [ + {"variable" : "mass", "title" : "mENuE" , "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "psumsq", "title" : "qSq", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "mtmLab", "title" : "mtmLab", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "mtm2", "title" : "mtm_2", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "pz", "title" : "mtm_2", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "coshel", "title" : "cosHel", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "abscoshel", "title" : "absCosHel", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosTheta", "title" : "cosTheta", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "prob", "title" : "Prob", "d1" : 0, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "cosHelTau", "" : "cosHelTau", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosHelTau", "" : "cosHelTau", "d1" : 1, "d2" : 1, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "E", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "E_over_Eparent", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "p", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "pz", "title" : "mtm_z", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "cosTheta", "" : "cosTheta", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "phi", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 180.0, "xmax" : 180}, + {"variable" : "E_daug1", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "E_over_Eparent_daug1", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "p_daug1", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "pz_daug1", "title" : "mtm_z", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "cosTheta_daug1", "" : "cosTheta", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "phi_daug1", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 180.0, "xmax" : 180}, + {"variable" : "E_daug1", "title" : "", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "E_over_Eparent_daug1", "title" : "", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "p_daug1", "title" : "", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "pz_daug1", "title" : "mtm_z", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "cosTheta_daug1", "" : "cosTheta", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "phi_daug1", "title" : "", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 180.0, "xmax" : 180}, + ], + "outfile" : "HmePythiaTau_BtoTaunu_pi.root", + "reference" : "RefHmePythiaTau_BtoTaunu_pi.root" +} diff --git a/test/jsonFiles/HmePythiaTau_enu.json b/test/jsonFiles/HmePythiaTau_enu.json index 6533038..29859b1 100644 --- a/test/jsonFiles/HmePythiaTau_enu.json +++ b/test/jsonFiles/HmePythiaTau_enu.json @@ -1,15 +1,26 @@ { "parent" : "tau-", "daughters" : ["nu_tau", "e-", "anti-nu_e"], "models" : ["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} + {"variable" : "mass", "title" : "mENuE" , "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "psumsq", "title" : "qSq", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "coshel", "title" : "cosE", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "prob", "title" : "Prob", "d1" : 0, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "E", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "E_over_Eparent", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "p", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "pz", "title" : "mtm_z", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "cosTheta", "" : "cosTheta", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "phi", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 180.0, "xmax" : 180}, + {"variable" : "E", "title" : "", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "p", "title" : "", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "pz", "title" : "mtm_z", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "cosTheta", "" : "cosTheta", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "phi", "title" : "", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : -180.0, "xmax" : 180}, ], "outfile" : "HmePythiaTau_enu.root", "reference" : "RefHmePythiaTau_enu.root" } diff --git a/test/jsonFiles/TAUOLATau_BtoTaunu_pi.json b/test/jsonFiles/TAUOLATau_BtoTaunu_pi.json new file mode 100644 index 0000000..e1f3994 --- /dev/null +++ b/test/jsonFiles/TAUOLATau_BtoTaunu_pi.json @@ -0,0 +1,42 @@ +{ + "parent" : "B+", + "daughters" : ["tau+", "nu_tau"], + "grand_daughters" : [["pi+", "anti-nu_tau"], []], + "models" : ["SLN", "TAUSCALARNU", ""], + "parameters" : [[], [], []], + "extras" : ["Define dm_incohMix_B0 0", "noPhotos"], + "events" : 100000, + "histograms" : [ + {"variable" : "mass", "title" : "mENuE" , "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "psumsq", "title" : "qSq", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "mtmLab", "title" : "mtmLab", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "mtm2", "title" : "mtm_2", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "pz", "title" : "mtm_2", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "coshel", "title" : "cosHel", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "abscoshel", "title" : "absCosHel", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosTheta", "title" : "cosTheta", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "prob", "title" : "Prob", "d1" : 0, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "cosHelTau", "" : "cosHelTau", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "cosHelTau", "" : "cosHelTau", "d1" : 1, "d2" : 1, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "E", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "E_over_Eparent", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "p", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "pz", "title" : "mtm_z", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "cosTheta", "" : "cosTheta", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "phi", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 180.0, "xmax" : 180}, + {"variable" : "E_daug1", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "E_over_Eparent_daug1", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "p_daug1", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "pz_daug1", "title" : "mtm_z", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "cosTheta_daug1", "" : "cosTheta", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "phi_daug1", "title" : "", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 180.0, "xmax" : 180}, + {"variable" : "E_daug1", "title" : "", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "E_over_Eparent_daug1", "title" : "", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 1.0}, + {"variable" : "p_daug1", "title" : "", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "pz_daug1", "title" : "mtm_z", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3}, + {"variable" : "cosTheta_daug1", "" : "cosTheta", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}, + {"variable" : "phi_daug1", "title" : "", "d1" : 2, "d2" : 2, "nbins" : 50, "xmin" : 180.0, "xmax" : 180}, + ], + "outfile" : "TAUOLATau_BtoTaunu_pi.root", + "reference" : "RefTAUOLATau_BtoTaunu_pi.root" +}