Page MenuHomeHEPForge
Files F19244605

No OneTemporary
Actions

Size
19 KB
Referenced Files
None
Subscribers
None

View Options

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 <https://www.gnu.org/licenses/>. *
***********************************************************************/
#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 <fstream>
#include <iostream>
-#include <memory>
#include <list>
+#include <memory>
#include <vector>
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<std::string> 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<double> 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<std::string>& daugNames,
- const std::string& modelName,
- const std::vector<double>& parameters ) const
+std::string testDecayModel::createDecFile(
+ const std::string& parent, const std::vector<std::string>& daugNames,
+ const std::string& modelName, const std::vector<double>& 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<EvtMTRandomEngine>();
EvtAbsRadCorr* radCorrEngine = nullptr;
std::list<EvtDecayBase*> 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<TH1*>( 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<TH1*>(
+ 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 <https://www.gnu.org/licenses/>. *
***********************************************************************/
#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 <string>
#include <vector>
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<std::string>& daugNames,
- const std::string& model,
- const std::vector<double>& parameters ) const;
+ const std::vector<std::string>& daugNames,
+ const std::string& model,
+ const std::vector<double>& 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<testInfo> m_infoVect;
std::vector<TH1D*> m_histVect;
size_t m_nHistos;
-
-
};
#endif

File Metadata

Mime Type
text/x-diff
Expires
Tue, Sep 30, 4:44 AM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6564794
Default Alt Text
(19 KB)

Event Timeline

HEPForge.org · .