diff --git a/test/jsonFiles/HmePythiaTau_BtoTauTau.json b/test/jsonFiles/HmePythiaTau_BtoTauTau.json
new file mode 100644
index 0000000..7517179
--- /dev/null
+++ b/test/jsonFiles/HmePythiaTau_BtoTauTau.json
@@ -0,0 +1,21 @@
+{
+ "parent" : "B0",
+ "daughters" : ["tau+", "tau-"],
+ "grand_daughters" : [["anti-nu_tau", "pi+", "pi+", "pi-"], ["nu_tau", "pi-", "pi-", "pi+"]],
+ "models" : ["PHSP", "HMEPYTHIA", "HMEPYTHIA"],
+ "parameters" : [[], [1541], [1541]],
+ "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" : "mtmLab", "title" : "mtmLab", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.0, "xmax" : 3},
+ {"variable" : "mtm2", "title" : "mtm_2", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.0, "xmax" : 3},
+ {"variable" : "pz", "title" : "mtm_2", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.0, "xmax" : 3},
+ {"variable" : "coshel", "title" : "cosHel", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+ {"variable" : "abscoshel", "title" : "absCosHel", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+ {"variable" : "cosTheta", "title" : "cosTheta", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+ {"variable" : "prob", "title" : "Prob", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0}
+ ],
+ "outfile" : "HmePythiaTau_BtoTauTau.root",
+ "reference" : "RefHmePythiaTau_BtoTauTau.root"
+}
diff --git a/test/jsonFiles/HmePythiaTau1.json b/test/jsonFiles/HmePythiaTau_enu.json
similarity index 72%
rename from test/jsonFiles/HmePythiaTau1.json
rename to test/jsonFiles/HmePythiaTau_enu.json
index fb4c71d..6533038 100644
--- a/test/jsonFiles/HmePythiaTau1.json
+++ b/test/jsonFiles/HmePythiaTau_enu.json
@@ -1,15 +1,15 @@
{
"parent" : "tau-",
- "daughters" : ["e-", "anti-nu_e", "nu_tau"],
- "model" : "HMEPYTHIA",
- "parameters" : [1531],
+ "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}
],
- "outfile" : "HmePythiaTau1.root",
- "reference" : "RefHmePythiaTau1.root"
+ "outfile" : "HmePythiaTau_enu.root",
+ "reference" : "RefHmePythiaTau_enu.root"
}
diff --git a/test/testDecayModel.cc b/test/testDecayModel.cc
index e968c8c..a1f75cb 100644
--- a/test/testDecayModel.cc
+++ b/test/testDecayModel.cc
@@ -1,348 +1,745 @@
/***********************************************************************
* 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/EvtVector4R.hh"
#ifdef EVTGEN_EXTERNAL
#include "EvtGenExternal/EvtExternalGenList.hh"
#endif
#include
#include
#include
#include
#include
testDecayModel::testDecayModel( const std::string& jsonFile ) :
m_jsonFile( jsonFile ), m_infoVect(), m_histVect(), m_nHistos( 0 )
{
}
void testDecayModel::run()
{
+ /*! Load input file in json format. */
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() );
+ Json::Value daughters = config["daughters"];
+ std::vector daughterNames;
+ for ( auto daugN : daughters ) {
+ daughterNames.push_back( daugN.asString() );
}
- const std::string modelName = config["model"].asString();
+ Json::Value grandDaughters = config["grand_daughters"];
+ std::vector> grandDaughterNames;
+ for ( auto iGrandDaughters : grandDaughters ) {
+ std::vector iGrandDaughterNames;
+ for ( auto ijGrandDaughter : iGrandDaughters ) {
+ iGrandDaughterNames.push_back( ijGrandDaughter.asString() );
+ }
+ grandDaughterNames.push_back( iGrandDaughterNames );
+ }
+
+ // std::cout << "grandDaughterNames: " << grandDaughterNames[0][0] << std::endl;
- Json::Value jParameters = config["parameters"];
- std::vector modelPars;
- for ( auto par : jParameters ) {
- modelPars.push_back( par.asDouble() );
+ Json::Value models = config["models"];
+ std::vector modelNames;
+ for ( auto iModel : models ) {
+ modelNames.push_back( iModel.asString() );
}
- // Create decay file based on json file input
- const std::string decFile = createDecFile( parentName, daugNames, modelName,
- modelPars );
+ Json::Value parameters = config["parameters"];
+ std::vector> modelParameters;
+ for ( auto iParameters : parameters ) {
+ std::vector iModelParameters;
+ for ( auto ijParameter : iParameters ) {
+ iModelParameters.push_back( ijParameter.asDouble() );
+ }
+ modelParameters.push_back( iModelParameters );
+ }
+
+ Json::Value extras = config["extras"];
+ std::vector extraCommands;
+ for ( auto iExtra : extras ) {
+ extraCommands.push_back( iExtra.asString() );
+ }
+
+ /*! Creates a decay file based on json file input. */
+ const std::string decFile = createDecFile( parentName, daughterNames,
+ grandDaughterNames, modelNames,
+ modelParameters, extraCommands );
int nEvents = config["events"].asInt();
std::cout << "Number of events = " << nEvents << std::endl;
- // ROOT output file
+ bool debugFlag = config["debug_flag"].asBool();
+
+ /*! Define the root output file and histograms to be saved. */
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 );
+ /*! Generate events and fill histograms. */
+ generateEvents( decFile, parentName, nEvents, debugFlag );
- // Normalize histograms
+ // Normalize histograms.
for ( auto hist : m_histVect ) {
double area = hist->Integral();
if ( area > 0.0 ) {
hist->Scale( 1.0 / area );
}
}
- // Compare with reference histograms
+ /*! Compare with reference histograms. */
const std::string refFileName = config["reference"].asString();
compareHistos( refFileName );
- // Write output
+ // Write output.
outFile->cd();
for ( auto hist : m_histVect ) {
hist->Write();
}
outFile->Close();
+ std::cout << "Created output file: " << outFileName.c_str() << std::endl;
}
std::string testDecayModel::createDecFile(
- const std::string& parent, const std::vector& daugNames,
- const std::string& modelName, const std::vector& parameters ) const
+ const std::string& parent, const std::vector& daughterNames,
+ const std::vector>& grandDaughterNames,
+ const std::vector& modelNames,
+ const std::vector>& parameters,
+ const std::vector& extras ) const
{
// Create (or overwrite) the decay file
- std::string decName( modelName );
+ std::string decName( "test" );
+ for ( auto iModelName : modelNames )
+ decName += "_" + iModelName;
decName += ".dec";
std::ofstream decFile;
decFile.open( decName.c_str() );
+ // Create aliases if needed
+ std::vector aliasPrefix;
+ for ( long unsigned int daughter_index = 0;
+ daughter_index < daughterNames.size(); daughter_index++ ) {
+ if ( !grandDaughterNames.empty() &&
+ !grandDaughterNames[daughter_index].empty() ) {
+ decFile << "Alias My" << daughterNames[daughter_index] << " "
+ << daughterNames[daughter_index] << std::endl;
+ aliasPrefix.push_back( "My" );
+ } else {
+ aliasPrefix.push_back( "" );
+ }
+ }
+
+ for ( auto iExtra : extras ) {
+ decFile << iExtra << std::endl;
+ }
+
+ // Parent decay
decFile << "Decay " << parent << std::endl;
decFile << "1.0";
- for ( auto daug : daugNames ) {
- decFile << " " << daug;
+
+ for ( long unsigned int daughter_index = 0;
+ daughter_index < daughterNames.size(); daughter_index++ ) {
+ decFile << " " << aliasPrefix[daughter_index]
+ << daughterNames[daughter_index];
}
- decFile << " " << modelName;
- for ( auto par : parameters ) {
+
+ decFile << " " << modelNames[0];
+
+ for ( auto par : parameters[0] ) {
decFile << " " << par;
}
+
decFile << ";" << std::endl;
decFile << "Enddecay" << std::endl;
+
+ // Daughter decays into granddaughters
+ for ( long unsigned int daughter_index = 0;
+ daughter_index < grandDaughterNames.size(); daughter_index++ ) {
+ if ( grandDaughterNames[daughter_index].empty() )
+ continue;
+ decFile << "Decay " << aliasPrefix[daughter_index]
+ << daughterNames[daughter_index] << std::endl;
+ decFile << "1.0";
+ for ( long unsigned int grandDaughter_index = 0;
+ grandDaughter_index < grandDaughterNames[daughter_index].size();
+ grandDaughter_index++ ) {
+ decFile << " "
+ << grandDaughterNames[daughter_index][grandDaughter_index];
+ }
+ decFile << " " << modelNames[daughter_index + 1];
+ for ( auto par : parameters[daughter_index + 1] ) {
+ 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 )
{
// 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 );
}
m_nHistos = m_histVect.size();
}
void testDecayModel::generateEvents( const std::string& decFile,
- const std::string& parentName, int nEvents )
+ const std::string& parentName, int nEvents,
+ bool debug_flag )
{
// 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 );
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 );
+ // To debug
+
+ if ( debug_flag ) {
+ std::cout << "Parent PDG code: " << parent->getPDGId()
+ << " has daughters " << parent->getNDaug() << std::endl;
+ for ( size_t iDaughter = 0; iDaughter < parent->getNDaug();
+ iDaughter++ ) {
+ std::cout << "Parent PDG code of daughter " << iDaughter
+ << " : " << parent->getDaug( iDaughter )->getPDGId()
+ << " has daughters "
+ << parent->getDaug( iDaughter )->getNDaug()
+ << std::endl;
+
+ for ( size_t iGrandDaughter = 0;
+ iGrandDaughter < parent->getDaug( iDaughter )->getNDaug();
+ iGrandDaughter++ ) {
+ std::cout << "Parent PDG code of grand daughter "
+ << iGrandDaughter << " : "
+ << parent->getDaug( iDaughter )
+ ->getDaug( iGrandDaughter )
+ ->getPDGId()
+ << " has daughters "
+ << parent->getDaug( iDaughter )
+ ->getDaug( iGrandDaughter )
+ ->getNDaug()
+ << std::endl;
+ }
+ }
+ }
+
// 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();
+ // parent->printTree();
// Cleanup
parent->deleteTree();
}
}
double testDecayModel::getValue( EvtParticle* parent, const std::string& varName,
const int d1, const int d2 ) const
{
- double value( 0.0 );
+ double value = std::numeric_limits::quiet_NaN();
if ( !parent ) {
return value;
}
- // Daughters
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;
+ // If variable name contains "_daugX", we are interested in daugthers of daughter X
+ // Else we are interested in daughters
+ EvtParticle* selectedParent = parent;
+ std::string selectedVarName = varName;
+ if ( varName.find( "_daug" ) != std::string::npos ) {
+ // Get daughter index from last charachter in string
+ const int iDaughter = varName.back() - '0';
+ selectedVarName = varName.substr( 0, varName.size() - 6 );
+ if ( iDaughter > 0 && iDaughter <= NDaugMax )
+ selectedParent = parent->getDaug( iDaughter - 1 );
+ else
+ return value;
+ }
+
+ const int sel_NDaugMax( selectedParent->getNDaug() );
+ const EvtParticle* par1 = d1 > 0 && d1 <= sel_NDaugMax
+ ? selectedParent->getDaug( d1 - 1 )
+ : nullptr;
+ const EvtParticle* par2 = d2 > 0 && d2 <= sel_NDaugMax
+ ? selectedParent->getDaug( d2 - 1 )
+ : nullptr;
+ const EvtParticle* par3 = sel_NDaugMax > 0
+ ? selectedParent->getDaug( sel_NDaugMax - 1 )
+ : nullptr;
// 4-momenta (in parent rest frame)
const EvtVector4R p1 = par1 != nullptr ? par1->getP4() : EvtVector4R();
const EvtVector4R p2 = par2 != nullptr ? par2->getP4() : EvtVector4R();
+ const EvtVector4R p3 = par3 != nullptr ? par3->getP4() : EvtVector4R();
+
+ // 4-momenta (in parent lab frame)
+ const EvtVector4R p1_lab = par1 != nullptr ? par1->getP4Lab() : EvtVector4R();
+ const EvtVector4R p3_lab = par3 != nullptr ? par3->getP4Lab() : EvtVector4R();
- if ( !varName.compare( "mass" ) ) {
+ if ( !selectedVarName.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();
}
- } else if ( !varName.compare( "psumsq" ) ) {
+ } else if ( !selectedVarName.compare( "psumsq" ) ) {
// Invariant momentum sum squared
value = ( p1 + p2 ).mass2();
- } else if ( !varName.compare( "pdiffsq" ) ) {
+ } else if ( !selectedVarName.compare( "pdiffsq" ) ) {
// Invariant momentum difference squared
value = ( p1 - p2 ).mass2();
- } else if ( !varName.compare( "mtm" ) ) {
+ } else if ( !selectedVarName.compare( "mass3" ) ) {
+ // Invariant mass of 3 daughters
+ value = ( p1 + p2 + p3 ).mass();
+
+ } else if ( !selectedVarName.compare( "cosTheta3" ) ) {
+ // Cosines of polar angle of first daughter
+ const EvtVector4R p123 = p1 + p2 + p3;
+ if ( p123.d3mag() > 0.0 ) {
+ value = p123.get( 3 ) / p123.d3mag();
+ }
+
+ } else if ( !selectedVarName.compare( "mtmLab" ) ) {
+ // Momentum of particle d1
+ value = p1_lab.d3mag();
+
+ } else if ( !selectedVarName.compare( "mtm" ) ) {
// Momentum of particle d1
value = p1.d3mag();
- } else if ( !varName.compare( "coshel" ) ) {
+ } else if ( !selectedVarName.compare( "mtm2" ) ) {
+ // Momentum of particle d1
+ double p1_lab_x = p1_lab.get( 1 );
+ double p1_lab_y = p1_lab.get( 2 );
+ double p1_lab_z = p1_lab.get( 3 );
+ double p1_lab_mag = p1_lab_x * p1_lab_x + p1_lab_y * p1_lab_y +
+ p1_lab_z * p1_lab_z;
+
+ value = p1_lab_mag;
+
+ } else if ( !selectedVarName.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 );
}
- } else if ( !varName.compare( "prob" ) ) {
+ } else if ( !selectedVarName.compare( "abscoshel" ) ) {
+ // Absolute 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 = std::abs( p1Res.dot( p12 ) / ( p1ResMag * p12Mag ) );
+ }
+
+ } else if ( !selectedVarName.compare( "cosHelTau" ) ) {
+ // Works only for B -> X nu_1 tau -> pi nu_2.
+ // Cosine of helicity angle between pi momentum and oposite W momentum -(nu_1 + tau) in tau rest frame
+ // Index d1 must match with tau
+ // p3 (momentum of last daughter), is taken as neutrino momentum
+
+ // W momentum
+ const EvtVector4R p_W = p1_lab + p3_lab;
+
+ // Index d2 must match the index of the pion (daughter of tau)
+ const EvtParticle* pion = selectedParent->getDaug( d1 - 1 )->getDaug(
+ d2 - 1 );
+ const EvtVector4R p_pion = pion != nullptr ? pion->getP4Lab()
+ : EvtVector4R();
+
+ // Boost vector to tau frame
+ const EvtVector4R boost( p1_lab.get( 0 ), -p1_lab.get( 1 ),
+ -p1_lab.get( 2 ), -p1_lab.get( 3 ) );
+
+ // Boost both momenta to tau frame
+ const EvtVector4R p_W_boosted = boostTo( p_W, boost );
+ const EvtVector4R p_pion_boosted = boostTo( p_pion, boost );
+ // Cosine of angle between oposite W momentum and pion momentum in tau frame
+ double p_W_boostedMag = p_W_boosted.d3mag();
+ double p_pion_boostedMag = p_pion_boosted.d3mag();
+
+ if ( p_W_boostedMag > 0.0 && p_pion_boostedMag > 0.0 ) {
+ value = -p_W_boosted.dot( p_pion_boosted ) /
+ ( p_W_boostedMag * p_pion_boostedMag );
+ }
+
+ } else if ( !selectedVarName.compare( "cosHelDiTau" ) ) {
+ // Works only for B -> tau (pi nu) tau -> (pi nu).
+ // Cosine of helicity angle between pi momentum and oposite B momentum in tau rest frame
+ // Index d1 must match with tau
+
+ // B momentum
+ const EvtVector4R p_B = parent->getP4Lab();
+
+ // Index d2 must match the index of the pion (daughter of tau)
+ const EvtParticle* pion = selectedParent->getDaug( d1 - 1 )->getDaug(
+ d2 - 1 );
+ const EvtVector4R p_pion = pion != nullptr ? pion->getP4Lab()
+ : EvtVector4R();
+
+ // Boost vector to tau frame
+ const EvtVector4R boost( p1_lab.get( 0 ), -p1_lab.get( 1 ),
+ -p1_lab.get( 2 ), -p1_lab.get( 3 ) );
+
+ // Boost both momenta to tau frame
+ const EvtVector4R p_B_boosted = boostTo( p_B, boost );
+ const EvtVector4R p_pion_boosted = boostTo( p_pion, boost );
+ // Cosine of angle between oposite W momentum and pion momentum in tau frame
+ double p_B_boostedMag = p_B_boosted.d3mag();
+ double p_pion_boostedMag = p_pion_boosted.d3mag();
+
+ if ( p_B_boostedMag > 0.0 && p_pion_boostedMag > 0.0 ) {
+ value = -p_B_boosted.dot( p_pion_boosted ) /
+ ( p_B_boostedMag * p_pion_boostedMag );
+ }
+
+ } else if ( !selectedVarName.compare( "cosHelDiTau_over05" ) ) {
+ // Works only for B -> tau (pi nu) tau -> (pi nu).
+ // Cosine of helicity angle between pi momentum and oposite B momentum in tau rest frame
+ // Index d1 must match with tau
+
+ // B momentum
+ const EvtVector4R p_B = parent->getP4Lab();
+
+ // Index d2 must match the index of the pion (daughter of tau)
+ const EvtParticle* pion_1 = selectedParent->getDaug( 0 )->getDaug( d2 -
+ 1 );
+ const EvtVector4R p_pion_1 = pion_1 != nullptr ? pion_1->getP4Lab()
+ : EvtVector4R();
+
+ const EvtVector4R p_first_daughter =
+ selectedParent->getDaug( 0 )->getP4Lab();
+
+ // Boost vector to tau frame
+ const EvtVector4R boost_1( p_first_daughter.get( 0 ),
+ -p_first_daughter.get( 1 ),
+ -p_first_daughter.get( 2 ),
+ -p_first_daughter.get( 3 ) );
+
+ // Boost both momenta to tau frame
+ const EvtVector4R p_B_boosted_1 = boostTo( p_B, boost_1 );
+ const EvtVector4R p_pion_boosted_1 = boostTo( p_pion_1, boost_1 );
+ // Cosine of angle between oposite W momentum and pion momentum in tau frame
+
+ double hel1 = -p_B_boosted_1.dot( p_pion_boosted_1 ) /
+ ( p_B_boosted_1.d3mag() * p_pion_boosted_1.d3mag() );
+
+ if ( hel1 > 0.5 ) {
+ // Works only for B -> tau (pi nu) tau -> (pi nu).
+ // Cosine of helicity angle between pi momentum and oposite B momentum in tau rest frame
+ // Index d1 must match with tau
+
+ // Index d2 must match the index of the pion (daughter of tau)
+ const EvtParticle* pion = selectedParent->getDaug( d1 - 1 )->getDaug(
+ d2 - 1 );
+ const EvtVector4R p_pion = pion != nullptr ? pion->getP4Lab()
+ : EvtVector4R();
+
+ // Boost vector to tau frame
+ const EvtVector4R boost( p1_lab.get( 0 ), -p1_lab.get( 1 ),
+ -p1_lab.get( 2 ), -p1_lab.get( 3 ) );
+
+ // Boost both momenta to tau frame
+ const EvtVector4R p_B_boosted = boostTo( p_B, boost );
+ const EvtVector4R p_pion_boosted = boostTo( p_pion, boost );
+ // Cosine of angle between oposite W momentum and pion momentum in tau frame
+ double p_B_boostedMag = p_B_boosted.d3mag();
+ double p_pion_boostedMag = p_pion_boosted.d3mag();
+
+ if ( p_B_boostedMag > 0.0 && p_pion_boostedMag > 0.0 ) {
+ value = -p_B_boosted.dot( p_pion_boosted ) /
+ ( p_B_boostedMag * p_pion_boostedMag );
+ }
+ }
+ } else if ( !selectedVarName.compare( "cosHelDiTau_under05" ) ) {
+ // Works only for B -> tau (pi nu) tau -> (pi nu).
+ // Cosine of helicity angle between pi momentum and oposite B momentum in tau rest frame
+ // Index d1 must match with tau
+
+ // B momentum
+ const EvtVector4R p_B = parent->getP4Lab();
+
+ // Index d2 must match the index of the pion (daughter of tau)
+ const EvtParticle* pion_1 = selectedParent->getDaug( 0 )->getDaug( d2 -
+ 1 );
+ const EvtVector4R p_pion_1 = pion_1 != nullptr ? pion_1->getP4Lab()
+ : EvtVector4R();
+
+ const EvtVector4R p_first_daughter =
+ selectedParent->getDaug( 0 )->getP4Lab();
+
+ // Boost vector to tau frame
+ const EvtVector4R boost_1( p_first_daughter.get( 0 ),
+ -p_first_daughter.get( 1 ),
+ -p_first_daughter.get( 2 ),
+ -p_first_daughter.get( 3 ) );
+
+ // Boost both momenta to tau frame
+ const EvtVector4R p_B_boosted_1 = boostTo( p_B, boost_1 );
+ const EvtVector4R p_pion_boosted_1 = boostTo( p_pion_1, boost_1 );
+ // Cosine of angle between oposite W momentum and pion momentum in tau frame
+
+ double hel1 = -p_B_boosted_1.dot( p_pion_boosted_1 ) /
+ ( p_B_boosted_1.d3mag() * p_pion_boosted_1.d3mag() );
+
+ if ( hel1 < -0.5 ) {
+ // Works only for B -> tau (pi nu) tau -> (pi nu).
+ // Cosine of helicity angle between pi momentum and oposite B momentum in tau rest frame
+ // Index d1 must match with tau
+
+ // Index d2 must match the index of the pion (daughter of tau)
+ const EvtParticle* pion = selectedParent->getDaug( d1 - 1 )->getDaug(
+ d2 - 1 );
+ const EvtVector4R p_pion = pion != nullptr ? pion->getP4Lab()
+ : EvtVector4R();
+
+ // Boost vector to tau frame
+ const EvtVector4R boost( p1_lab.get( 0 ), -p1_lab.get( 1 ),
+ -p1_lab.get( 2 ), -p1_lab.get( 3 ) );
+
+ // Boost both momenta to tau frame
+ const EvtVector4R p_B_boosted = boostTo( p_B, boost );
+ const EvtVector4R p_pion_boosted = boostTo( p_pion, boost );
+ // Cosine of angle between oposite W momentum and pion momentum in tau frame
+ double p_B_boostedMag = p_B_boosted.d3mag();
+ double p_pion_boostedMag = p_pion_boosted.d3mag();
+
+ if ( p_B_boostedMag > 0.0 && p_pion_boostedMag > 0.0 ) {
+ value = -p_B_boosted.dot( p_pion_boosted ) /
+ ( p_B_boostedMag * p_pion_boostedMag );
+ }
+ }
+ } else if ( !selectedVarName.compare( "cosTheta" ) ) {
+ // Cosines of polar angle of first daughter
+ double p1_lab_mag = p1_lab.d3mag();
+ if ( p1_lab_mag > 0.0 ) {
+ value = p1_lab.get( 3 ) / p1_lab_mag;
+ }
+
+ } else if ( !selectedVarName.compare( "phi" ) ) {
+ // Azimuthal angle of first daughter
+ double p1_lab_mag = p1_lab.d3mag();
+ if ( p1_lab_mag > 0.0 ) {
+ value = atan2( p1_lab.get( 1 ), p1_lab.get( 2 ) ) * 180 / M_PI;
+ }
+
+ } else if ( !selectedVarName.compare( "E" ) ) {
+ // Energy of first daughter
+ value = p1_lab.get( 0 );
+
+ } else if ( !selectedVarName.compare( "E_over_Eparent" ) ) {
+ // Energy of first daughter
+ value = p1_lab.get( 0 ) / selectedParent->getP4Lab().get( 0 );
+
+ } else if ( !selectedVarName.compare( "E_over_Eparent_over05" ) ) {
+ // First daughter E_over_Eparent
+ double E_over_Eparent_2 =
+ parent->getDaug( 0 )->getDaug( d2 - 1 )->getP4Lab().get( 0 ) /
+ parent->getDaug( 0 )->getP4Lab().get( 0 );
+
+ if ( E_over_Eparent_2 > 0.5 ) {
+ value = p1_lab.get( 0 ) / selectedParent->getP4Lab().get( 0 );
+ }
+
+ } else if ( !selectedVarName.compare( "p" ) ) {
+ // Energy of first daughter
+ value = p1_lab.d3mag();
+
+ } else if ( !selectedVarName.compare( "pz" ) ) {
+ // Momentum of particle d1
+ value = p1_lab.get( 3 );
+
+ } else if ( !selectedVarName.compare( "prob" ) ) {
// Decay probability
- double* dProb = parent->decayProb();
+ double* dProb = selectedParent->decayProb();
if ( dProb ) {
value = *dProb;
}
}
return value;
}
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;
}
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[] )
{
if ( argc != 2 ) {
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 5fddd8c..b8910f0 100644
--- a/test/testDecayModel.hh
+++ b/test/testDecayModel.hh
@@ -1,80 +1,84 @@
/***********************************************************************
* 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 "TFile.h"
#include "TH1D.h"
#include "json/json.h"
+#include "math.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 )
{
;
}
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();
private:
- std::string createDecFile( const std::string& parent,
- const std::vector& daugNames,
- const std::string& model,
- const std::vector& parameters ) const;
+ std::string createDecFile(
+ const std::string& parent, const std::vector& daughterNames,
+ const std::vector>& grandDaughterNames,
+ const std::vector& models,
+ const std::vector>& parameters,
+ const std::vector& extras ) const;
void defineHistos( Json::Value& config, TFile* theFile );
void generateEvents( const std::string& decFile,
- const std::string& parentName, int nEvents );
+ const std::string& parentName, int nEvents,
+ bool debugFlag );
double getValue( EvtParticle* rootPart, const std::string& varName,
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