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