diff --git a/History.txt b/History.txt
index b9d7f41..68a0917 100644
--- a/History.txt
+++ b/History.txt
@@ -1,486 +1,496 @@
//==========================================================================
//
// History file for EvtGen
//
//===========================================================================
+14th July 2017 John Back
+ Only create external generator objects if they don't already exist in
+ EvtExternalGenFactory. Modified configure script to work with Pythia 8.2x
+
+5th July 2017 Michal Kreps
+ Register the VTOSLL model.
+
+14th June 2017 John Back
+ Add isNeutralKaon() boolean function and corrected comments in EvtDDalitz.
+
8th May 2017 Michal Kreps
Fix bug in EvtbTosllVectorAmp to recognise Bs --> K*bar mu mu decay as
b --> d ll transition.
8th May 2017 Michal Kreps
Significantly simplify way how we decide on decay mode and daughters
ordering in DDalitz model. With new code by definition all orderings of
daughters in the decay file will yield same output.
4th May 2017 John Back
Further fixes to DDalitz particle ordering (including charge-conjugates):
Mode 5: D0 -> K- K0bar K+ and K+ K- K0bar
Mode 12: D0 -> pi0 pi- pi+ and pi+ pi0 pi-
Removed unneeded index ordering checks for mode 10 (D+ -> pi- pi+ pi+)
and mode 11 (Ds+ -> pi- pi+ pi+)
27th April 2017 John Back
Fixed DDalitz particle ordering for mode 7: D+ -> pi+ K- K+ and K+ pi+ K-
and their charge-conjugates
7th April 2017 John Back
Modified EvtGenExternal/EvtPythiaEngine to ensure that the EvtGen-based
instances of Pythia8 (for generic and alias decays) use the same
particle properties as defined by EvtGen, courtesy Patrick Robbe (LHCb).
5th April 2017 Michal Kreps
Fixed indexing in copy constructor of Evt3Rank3C, which would otherwise
produce an infinite loop; bug report from David Grellscheid.
3rd November 2016 John Back
Modified EvtFlatQ2 model to work for all B -> X lepton lepton modes, as
well as adding an extra phase space factor to correct for the dip at low
q^2, courtesy of Marcin Chrzaszcz & Thomas Blake (LHCb), which is enabled
by using "FLATQ2 1" instead of just "FLATQ2" in the decay file.
13th October 2016 John Back
Added the TauolaCurrentOption decfile keyword to select the hadronic
current in Tauola; default is the BaBar-tuned current option (int = 1).
EvtParticles can store double attributes using the functions
setAttributeDouble(name, double) and getAttributeDouble(name), which can
be useful for storing and retrieving amplitude weights, for example.
The analogous EvtParticle integer attribute interface remains unchanged:
setAttribute(name, int) and getAttribute(name).
13th September 2016 John Back
Modified EvtTauolaEngine to use internal Tauola spin matrices for
tau pair events by temporarily setting the PDG id of the mother as a
boson, keeping the same 4-momentum. The BaBar hadronic currents are now
used by default. Also added the ability to change some Tauola parameters
using the "Define" keyword in decay files. Added an example decay file
illustrating the new features: validation/TauolaFiles/Btautau.dec
9th September 2016 Michal Kreps
Reimplement code in EvtBTo3pi.F, EvtBTo3piMPP.F, EvtBTo3piP00.F and
EvtBToKpipi.F in C++ in order to remove dependence on Fortran compiler.
With this, there is no internal Fortran code in EvtGen.
R01-06-00-------------------------------------------------------------------
1st June 2016 John Back
New tag incorporating all changes below. Recommended external packages
are HepMC 2.06.09, pythia 8.186, Photos++ 3.61 and Tauola++ 1.1.5.
28th April 2016 Michal Kreps
For Ds+ --> 2pi+ pi- there was double counting of branching fraction
resulting in total branching fraction being 1.5 times larger than measured
one. Fix by revisiting submodes, which now fill total Ds --> 3pi.
25th April 2016 Michal Kreps
Added DECAY.DEC/XML, which contain updated semileptonic charm and beauty
branching fractions using the 2014 PDG, tuned to minimize disagreements
between measurements and EvtGen for both inclusive and exclusive decays.
Updated the evt.pdl particle properties file to the PDG 2014 edition.
Implemented new LQCD form factors for Lb --> L mu mu from arXiv paper
1602.01399 (EvtRareLbToLllFFlQCD); old LQCD form factors are removed.
18th March 2016 John Back
Fixed incorrect spinor algebra used in S1 -> 1/2 S2, 1/2 -> S3 S4 decays
in EvtDiracParticle::rotateToHelicityBasis() functions, courtesy of
Luis Miguel Garcia Martin and the IFIC Valencia LHCb group.
19th Feburary 2016 John Back
Fixed bug in the definition of the initial spinor term Sinit in
EvtRareLbToLll::HadronicAmpRS(), from Tom Blake (LHCb).
12th February 2016 John Back
From LHCb, added extensions to the EvtHQET2(FF) model for semitauonic
decays from Brian Hamilton, which needs a patch to EvtSemiLeptonicAmp
from Jack Wimberley to ensure that the q^2 range is physical when
finding the maximum amplitude probability.
2nd December 2015 John Back
From LHCb, added EvtKStopizmumu model for KS -> pi0 mu mu decays based on
JHEP08(1998)004, courtesy of Veronika Chobanova, Diego Martinez Santos
and Jeremy Dalseno. Added EvtConst::Fermi for Fermi coupling constant.
R01-05-00-------------------------------------------------------------------
21st October 2015 John Back
New tag incorporating all changes below. Recommended external packages
are HepMC 2.06.09, pythia 8.186, Photos++ 3.61 and Tauola++ 1.1.5.
Added the EvtB2MuMuMuNu model for simulating the very rare four-leptonic
decays B- -> mu+ mu- anti-nu_mu mu-, courtesy Nikolai Nikitin.
16th October 2015 John Back
Updated the configure script to automatically select the library names
for PHOTOS++; version 3.56 and below uses Fortran, version 3.61 and above
uses C++ only (default). Avoid using v3.60, since it does not work.
This needs the PHOTOS libraries built before EvtGen is configured.
Modified setupEvtGen.sh to use Photos++ v3.61.
7th October 2015 John Back
Updated EvtGenExternal/EvtPhotosEngine to check that additional particles
from the outgoing vertex are indeed (FSR) photons, since later versions of
PHOTOS introduce pair emission, where particles may not always be photons.
Added the genRootDecayChain.cc validation program to create ROOT files
containing information about the complete decay tree. Two example test
decay files BKstarGamma.dec and BuDst0rhop.dec can be used with this; the
first tests PHOTOS, the second looks at sequential decay chain storage.
The plotBKstarGamma.C ROOT macro can be used for B -> K* gamma plots.
2nd October 2015 John Back
Modified EvtSVPHelAmp and added a new EvtSVPHelCPMix model, implementing
the complete mixing phenomenology of Bs to vector gamma decays, courtesy
of Clara Remon (LHCb).
EvtD0mixDalitz code: cleanup, inverted q/p for decays of D0bar (simplifies
user decay files) and fixed y parameter bug, courtesy of Jordi Tico (LHCb).
Changed the initialisation order of the infrared cut-off in EvtPhotosEngine.
This actually has no effect, since the exponentiation function sets it to the
same 1e-7 value, but it's now in the correct order if we need to update it.
Removed all remaining obsolete pragma (Win32) warnings from some classes.
23rd September 2015 Michal Kreps
Reimplement the real Spence function in C++ and removed its fortran
implementation.
15th September 2015 Michal Kreps
Fixed accessed uninitialised memory in EvtPDL.cpp, line 213.
Modified the configure and setupEvtGen.sh scripts to work on Mac; needed
Mac compilation patch files added to the new "platform" subdirectory.
10th September 2015 John Back
Updated setupEvtGen.sh to use the recommended external packages:
HepMC 2.06.09, pythia 8.186, Photos++ 3.56 and Tauola++ 1.1.5.
Fixed form-factor calculations for the BTOSLLBALL model 6 used to generate
b -> sll decays, courtesy of Christoph Langenbruch and David Loh (LHCb).
Affects B->K*ll, B->rholl and B->omegall, particularly the electron modes.
In the validation directory, added runPhotosTest.sh for testing FSR in
Upsilon(4S) -> e+ e- decays, and changed the plot comparison scripts to
use the 2nd directory "oldRootFiles" (which could be a soft-link) for
including ROOT histograms made from a previous version of EvtGen.
27th August 2015 John Back
Added Mersenne-Twister random number generator (RNG) EvtMTRandomEngine.
It requires c++11 compiler features (>= gcc 4.7), which should
automatically be enabled by the configure script. Introduced the
preprocessor environment variable EVTGEN_CPP11 for c++11 features.
EvtMTRandomEngine is the default RNG for the validation and test examples
if c++11 features are enabled.
Added a phase-space test validation/genPHSP.sh and PhaseSpacePlots.C to
visually check the flatness of Dalitz plots in order to ensure that the
RNG is not producing biased results that depend on particle ordering.
Added the models EvtbsToLLLLAmp and EvtbsToLLLLHyperCP for
B0_q -> l+ l- l+ l- decays (SM and one supersymmetric scenario),
courtesy of Nikolai Nikitin and Konstantin Toms. Documentation provided
in doc/evt_BQTOLLLL_model.pdf and doc/evt_BQTOLLLLHYPERCP_model.pdf.
Changed the installation and set-up script name to be just setupEvtGen.sh;
it uses the VERSION variable to specify the required tag. List of tags
are available using either "svn ls -v http://svn.cern.ch/guest/evtgen/tags"
or by going to http://svn.cern.ch/guest/evtgen/tags in a web browser.
12th June 2015 John Back
Changed the width of chi_b1 in evt.pdl from 9.8928 GeV (!) to zero.
1st May 2015 John Back
Added Bc -> scalar ell nu (EvtBcSMuNu) and Bc -> tensor ell nu
(EvtBcTMuNu) decays, courtesy of Jack Wimberley, LHCb. Also included the
chi_c1 mode for EvtBcVMuNu.
R01-04-00-------------------------------------------------------------------
2nd April 2015 John Back
Removed the EvtStdlibRandomEngine class since this can produce biases
to kinematic distributions when one or more of the daughters is a
resonance, such as B0 -> K pi psi (thanks to Antonio Augusto Alves Jr
who discovered this issue). EvtSimpleRandomEngine is now the default
random number generator in the validation and test examples.
Incorporated several additions and modifications from LHCb:
a) From Michal Kreps, Tom Blake & Christoph Langenbruch,
added rare Lb --> Lambda^(*) ell ell models described in
arXiv:1108.6129, with various form factors from Gutsche et al.
(arXiv:1301.3737) and lattice QCD (arXiv:1212.4827)
b) From Andrew Crocombe, implemented Bs --> K* form factors
from Ball-Zwicky and z-parametrization form factors from
arXiv:1006.4945 for EvtbTosllBallFF
c) Christoph Langenbruch fixed the Bs -> phi ll form factors in
EvtbTosllBallFF; T3 showed a non-physical pole at very low
q2 which significantly affected the electron mode
d) From Michal Kreps, removed semicolons from wrong places to
clear warnings when compiled with the "-pedantic" option.
9th October 2014 John Back
Change svnweb.cern.ch to svn.cern.ch in the setup script.
1st April 2014 John Back
In EvtReport, modified the logging output severity status flags
to have the "EVTGEN_" prefix, e.g. INFO becomes EVTGEN_INFO.
The global report() function has been renamed to EvtGenReport().
31st March 2014 John Back
Added the ability to store named attributes for EvtParticles in the
form of a map<string, int>. The setAttribute(name, value) stores the
required value, while getAttribute(name) retrieves the integer value.
This is used in EvtPhotosEngine to specify the final-state radiation
"FSR" attribute to 1 for any additional photons (EvtPhotonParticles)
created by Photos++. It also stores the "ISR" attribute, but this
is always set to zero, since only FSR photons are created.
If the named attribute doesn't exist, then getAttribute() returns zero.
29th January 2014 Daniel Craik
Removed mass assertion on GS shape in EvtDalitzReso to allow it to also
be used for charged rho resonances.
27th January 2014 John Back
Minor corrections to Vub models to remove further gcc 4.8 warnings.
Updated configure script to work for MacOS clang (from Genser team).
R01-03-00-------------------------------------------------------------------
9th January 2014 John Back
New tag version "1.3.0", incorporating all changes below.
Replaced auto-install script to work with this version as well as
the latest versions of all external generator packages.
Updated README to mention the new CERN-based web pages for Photos++
and Tauola++.
8th January 2014 John Back
Fix gcc 4.6 and 4.8 compilation warnings, courtesy of
Patrick Robbe (LHCb); main changes are removal of unused variables.
Changed the EvtPythiaEngine class and configure script to use new
Pythia 8 header locations; Pythia 8.180 or above is now required.
7th January 2014 John Back
Modified EvtBCVFF to correct the Kiselev form factors from
Jack Wimberley (LHCb).
9th October 2013 Daniel Craik
Added Gounaris-Sakurai and Gaussian shapes to EvtGenericDalitz
and set sensible defaults for LASS parameters.
19th September 2013 John Back
Modified EvtGenExternal/EvtPythiaEngine to keep track of any new
particles that are added to the default Pythia database to avoid
duplicating particle/anti-particle entries that could override
previously defined Pythia decay chains.
18th September 2013 John Back
Added Mac OS flags for the configure script and src/Makefile.
15th July 2013 Daniel Craik
Added flag to turn on scaling of LASS amplitude by M/q in EvtDalitzReso
15th July 2013 Daniel Craik
EvtParserXML now accepts file names containing environment variables,
exponential non-resonant shape in EvtDalitzReso now defined as exp(-alpha*m^2),
LASS shape in EvtDalitzReso now takes a cutoff parameter
4th July 2013 Daniel Craik
Added LASS, exponential non-resonant and linear non-resonant shapes to EvtGenericDalitz.
3rd July 2013 Daniel Craik
Fixed auto-install script for R01-02-00.
1st July 2013 Daniel Craik
Added auto-install script for R01-02-00.
R01-02-00-------------------------------------------------------------------
15th May 2013 John Back
New tag, version "1.2.0", incorporating all changes below.
14th May 2013 Michal Kreps
Added Blatt-Weisskopf barrier factors up to L=5 in
EvtGenBase/EvtBlattWeisskopf::compute().
14th May 2013 John Back
Added additional entries (appended at the end) to the evt.pdl particle
data file courtesy of Romulus Godang and Belle II colleagues.
14th March 2013 John Back
Added the method EvtParticle::getPDGId() to get the PDG integer for a
particle directly (which just calls EvtPDL::getStdHep()).
Added a check in EvtPhotosEngine::doDecay to skip Photos if a given
particle has too many daughters (>= 10) to avoid a problem with a
hard coded upper limit in the Photos PHOENE subroutine.
2nd February 2013 Daniel Craik
Updated EvtDalitzTable to estimate probMax when it is missing from a
Dalitz model.
1st February 2013 John Back
Added the ability to read in Pythia 6 commands in ascii decay files in
EvtDecayTable::readDecayFile (this was already possible in xml files).
Modified the Photos++ engine default settings to be more suited to B
decays (from LHCb defaults).
31st January 2013 John Back
Added the ability to read in Pythia 8 commands in ascii decay files
in EvtDecayTable::readDecayFile. They can be set using the syntax:
"PythiaTypeParam module:variable=value", where Type = Generic, Alias or
Both for specifying whether the parameter is for the generic or alias
Pythia decay engines (or both). The 2nd argument must not contain any
spaces. Fixed the list of commands strings being used in the
EvtPythiaEngine class (i.e. Pythia parameters that can be set via decay
files).
31st January 2013 Daniel Craik
Added named parameters to various decay models.
30th January 2013 John Back
Fixed some of the parameter arguments used in the EvtSVSCPiso model.
24th January 2013 John Back
Set the Photos++ and Tauola++ models to use the EvtGen random number
engine when useEvtGenRandom is set to true in the EvtExternalGenList
constructor.
23rd January 2013 John Back
Added EvtGenExternal/EvtPythiaRandom to allow the use of the EvtGen
random number engine to also be used for the random engine for Pythia 8.
Added a boolean (useEvtGenRandom, default = true) within the
EvtExternalGenList constructor to use this feature.
18th December 2012 John Back
Corrected some wrong daughter ordering assignments for decay modes 5 and
12 in EvtDDalitz. Updated validation/DalitzDecays.xml to also contain
D decay mode 12, as well as various modes that may use K_S0 and K_L0.
Added validation/genDDalitzModes.sh and updated
validation/compareDalitz.C to do a complete comparison of the D Dalitz
modes with the xml versions.
11th December 2012 Daniel Craik
Updated the Xml parser to support named model parameters. Updated the
generic Dalitz model to use named model parameters as an example.
15th October 2012 John Back
Make EvtSimpleRandomEngine inherit from EvtRandomEngine to avoid
crash in EvtGen.cpp when no random engine is defined
(from Bjoern Spruck).
R01-01-00-------------------------------------------------------------------
4th October 2012 John Back
New tag, version "1.1.0", incorporating all changes below.
Provide proper default constructors for EvtVector4R and
EvtPhotonParticle. Modified the validation and test code to also
compile/link in the case of no external generators being included.
3rd October 2012 John Back
Corrected the t3 vector form factor values for the Ball-Zwicky '05
model (modelId = 6) in EvtbTosllBallFF::getVectorFF(), which
were set to t3tilde instead.
18th September 2012 John Back
Moved the external generator engines to a new sub-directory
EvtGenExternal. Building the code now creates 2 libraries:
libEvtGen.so (Base+Models) and libEvtGenExternal.so.
This allows anyone to ignore using the new external generators
if required (by not creating/loading the 2nd library).
Added prefix option to the configure script/Makefile to allow the user
to specify an installation directory for the include files, libraries,
DECAY.DEC and evt.pdl files (for Genser).
14th September 2012 Michal Kreps
Fixed the calculation of the angle between decay planes in the function
EvtKine::EvtDecayAngleChi. Fixed typo in EvtLb2Lll decay model. Only
some NP scenarious could be affected, SM one is definitely unaffected.
13th September 2012 John Back
Added the use of the environment variables EVTGEN_PHOTOS, EVTGEN_PYTHIA
and EVTGEN_TAUOLA to specify if the Photos, Pythia and/or Tauola engine
classes are used or not. These variables are set by the configure script,
depending if the library paths are specified for these generators.
R01-00-01--------------------------------------------------------------------
12th September 2012 John Back
New tag incorporating all changes below, since R01-00-00.
11th September 2012 John Back
Modified the Photos and Tauola engine classes to use the new
Photospp and Tauolapp namespaces that are present in the latest
versions of Photos++(3.5) and Tauola++(1.0.7). Updated the configure file
to get the correct location of the Tauola++ include files.
Added the D0->pi+pi-pi0 decay mode in EvtDDalitz from Marco Gersabeck
and Frederic Dreyer (LHCb).
Added new decay models/classes from Alexey Luchinsky (LHCb):
EvtBcVMuNu, EvtTVP, EvtWnPi, EvtSVP, EvtXPsiGamma, EvtBcVNpi
29th June 2012 John Back
Corrected mass(squared) variables filled in the Dalitz TTree in
validation/genExampleRootFiles.
15th May 2012 Daniel Craik
Updated EvtD0gammaDalitz to deal with D mesons from neutral B->DK
Added save function to validation/compareDalitz.C.
11th May 2012 Daniel Craik
Replaced BaBar specific configuration for BlattWeisskopf birth factors.
Updated XML conversion script to handle new configuration.
Fixed some bugs in the XML conversion script related to particle
modifications.
9th May 2012 Daniel Craik
Added latex documentation for xml decay files.
2nd May 2012 Daniel Craik
Added resDaughters attribute to the Dalitz resonance xml tag to
simplify defining symmetric resonances. Updated validation xml files to
use the new functionality.
27th April 2012 Daniel Craik
Upgraded EvtGenericDalitz to use EvtDalitzReso for resonances.
Added validation to compare EvtGenericDalitz to all 11 EvtDDalitz modes.
Added a root macro to quickly compare two Dalitz decays for validation.
24th April 2012 John Back
Solved two bugs in the EvtD0gammaDalitz model (from Jordi Tico, LHCb):
configuration of the conjugated model, and using only the B charge
to determine the model used, not the D flavour.
17th April 2012 Daniel Craik
Updated the GenericDalitz validation code to use the same probMax
values as DDalitz. Added XML decay file parsing to EvtGen::readUDecay.
Dec files are still the default.
30th March 2012 John Back
Update maximum probability values in EvtDDalitz::initProbMax()
for all DDalitz modes.
23rd March 2012 John Back
Added the EvtEta2MuMuGamma decay model from LHCb.
21st March 2012 John Back
Added EvtD0gammaDalitz decay model from LHCb.
20th March 2012 Daniel Craik
Added backwards compatibility for Pythia 6 commands in the XML
configuration. Updated decay file conversion tool to convert JetSetPar
lines to pythia6Param tags.
19th March 2012 Daniel Craik
Added infrastructure to pass commands to external generators.
XML config now takes Pythia8 configuration commands.
16th March 2012 Daniel Craik
Added the ability to define particles from the PDL for Dalitz decay
resonances instead of defining mass, width and spin seperately. Renamed
the lifetime attribute of Dalitz decay resonaces to width to avoid
confusion. Added further validation code for the generic Dalitz model.
15th March 2012 Daniel Craik
Added validation code for xml decay files and the generic Dalitz model.
R01-00-00 ------------------------------------------------------------------
6th March 2012 John Back
First official version for Genser (evtgen 1.0.0) that includes
support for external generators: Pythia8, Photos++ and Tauola++.
This also includes a preliminary version of creating Dalitz plot
decay models using EvtGenericDalitz.
diff --git a/configure b/configure
index 368e3f2..375e3b7 100755
--- a/configure
+++ b/configure
@@ -1,415 +1,420 @@
#! /bin/bash
#Package-specific:>
GENERATOR=EvtGen
#VERSION=11
# === Default values of input parameters ====
PATH_TO_EXTERNAL=
CLHEP_VERSION=
PREFIX=.
PYTHIADIR=${PYTHIADIR}
PHOTOSDIR=${PHOTOSDIR}
TAUOLADIR=${TAUOLADIR}
HEPMCDIR=${HEPMCDIR}
COMPMODE=OPT # Default: compiler optimization mode
EVTGEN_PHOTOS=
EVTGEN_PYTHIA=
EVTGEN_TAUOLA=
EVTGEN_EXTERNAL=
#<
#--------------------------------------------------------------
check_arguments () {
for arg in $* ; do
if [ "x$arg" = "x--enable-debug" ] ; then
COMPMODE=DBG
elif [ "x$arg" = "x--help" ] ; then
echo -e \
"\nUsage: ./configure [options] , where options are:\n\n"\
"--help : prints this help\n"\
"--prefix=[installation path] : specifies installation path\n"\
"--hepmcdir=[full directory path] : specifies location of HepMC base directory (required)\n"\
"--pythiadir=[full directory path] : specifies location of Pythia 8 base directory (optional)\n"\
"--photosdir=[full directory path] : specifies location of PHOTOS base directory (optional)\n"\
"--tauoladir=[full directory path] : specifies location of TAUOLA base directory (optional)\n"\
"--enable-debug : turns on debugging flags (optional)\n"
exit
else
if [ "x${arg}" = "x${arg/=/}" ] ; then
echo "${arg}: wrong option. Ignored." >&2
else
option=${arg/=*/}
value=${arg/*=/}
if [ "x${option}" = "x--prefix" ] ; then
PREFIX=${value}
elif [ "x${option}" = "x--pythiadir" ] ; then
PYTHIADIR=${value}
elif [ "x${option}" = "x--hepmcdir" ] ; then
HEPMCDIR=${value}
elif [ "x${option}" = "x--photosdir" ] ; then
PHOTOSDIR=${value}
elif [ "x${option}" = "x--tauoladir" ] ; then
TAUOLADIR=${value}
else
echo "${arg}: wrong option. Ignored." >&2
fi
fi
fi
done
return
}
#---------------------------------------
# Check that we can see (some of) the header files for the external generators.
# If not, print a warning message and ask the user to re-run configure
# specifying the base directory correct location(s).
check_dirs() {
echo "PREFIX = $PREFIX"
hepmcFile=${HEPMCDIR}/include/HepMC/GenParticle.h
pythiaFile=${PYTHIADIR}/include/Pythia8/Pythia.h
photosFile=${PHOTOSDIR}/include/Photos/Photos.h
tauolaFile=${TAUOLADIR}/include/Tauola/Tauola.h
stop=0
if [ -f $hepmcFile ]; then
echo "HEPMCDIR set to $HEPMCDIR";
else
echo "Mandatory: please specify the location of the HepMC base directory"\
"using the configure option --hepmcdir=[full directory path]"\
"or set the HEPMCDIR environment variable"
stop=1
fi
if [ -f $pythiaFile ]; then
echo "PYTHIADIR set to $PYTHIADIR";
EVTGEN_PYTHIA=1
EVTGEN_EXTERNAL=1
else
echo "To use PYTHIA: --pythiadir=[full directory path] or set the PYTHIADIR"\
"environment variable"
fi
if [ -f $photosFile ]; then
echo "PHOTOSDIR set to $PHOTOSDIR";
EVTGEN_PHOTOS=1
EVTGEN_EXTERNAL=1
else
echo "To use PHOTOS: --photosdir=[full directory path] or set the PHOTOSDIR"\
" environment variable"
fi
if [ -f $tauolaFile ]; then
echo "TAUOLADIR set to $TAUOLADIR";
EVTGEN_TAUOLA=1
EVTGEN_EXTERNAL=1
else
echo "To use TAUOLA: --tauoladir=[full directory path] or set the TAUOLADIR"\
" environment variable"
fi
if [ $stop = 1 ] ; then
exit 1
fi
}
#--------------------------------------------------------------
check_extlibs () {
# Set default shared libary lists for external packages
HEPMCLIBLIST="-lHepMC"
# For PHOTOS version 3.61 and above: default setting.
# Aside: avoid 3.60 since this does not work properly for EvtGen!
PHOTOSLIBLIST="-lPhotospp -lPhotosppHepMC"
# For PHOTOS version 3.56 and below: check to see
# if we need to modify this library list
if [ -f ${PHOTOSDIR}/lib/libPhotosFortran.so ]; then
PHOTOSLIBLIST="-lPhotosCxxInterface -lPhotosFortran"
fi
TAUOLALIBLIST="-lTauolaCxxInterface -lTauolaFortran"
- PYTHIALIBLIST="-lpythia8 -llhapdfdummy"
+
+ PYTHIALIBLIST="-lpythia8"
+ # For PYTHIA 8.1x and below, check to see if we need lhapdfdummy
+ if [ -f ${PYTHIADIR}/lib/liblhapdfdummy.so ]; then
+ PYTHIALIBLIST="-lpythia8 -llhapdfdummy"
+ fi
}
#-------------------------------------------------------------
check_arguments $*
check_dirs
check_extlibs
echo "Compilation mode is ${COMPMODE}"
#Package-specific:>
#
# User-changeable part, experts -----------------------------
#
CXX=g++
CXXFLAGS_OPT="-O2 -Wall"
CXXFLAGS_DBG="-g -Wall"
CFLAGS_OPT=-O2
CFLAGS_DBG=-g
FC=gfortran
FFLAGS_OPT="-O2 -Wuninitialized"
FFLAGS_DBG=-g
FFLAGSSHARED=-fPIC
CFLAGSSHARED=-fPIC
CXXFLAGSSHARED=-fPIC
LDFLAGSSHARED="${CXXFLAGS_OPT} -pthread -fPIC"
#
# Find platform.
#
ARCH=`uname`
theGcc=`g++ --version | grep '[0-9]\.[0-9]\.[0-9]' -o | head -1 | awk -F . '{print $1}'`
if [ ${theGcc} = 4 ]; then
ARCH=${ARCH}-gcc4
fi
if [ ${theGcc} = 5 ]; then
ARCH=${ARCH}-gcc5
fi
if [ ${theGcc} = 6 ]; then
ARCH=${ARCH}-gcc6
fi
echo "Platform is $ARCH"
#default platform settings:
FFLAGS="${FFLAGS_OPT}"
CFLAGS="${CFLAGS_OPT}"
CXXFLAGS="${CXXFLAGS_OPT}"
FLIBS="-lfrtbegin -lg2c"
SOFLAGS="-soname"
if [ ${COMPMODE} = OPT ]; then
FFLAGS="${FFLAGS_OPT}"
CFLAGS="${CFLAGS_OPT}"
CXXFLAGS="${CXXFLAGS_OPT}"
fi
if [ ${COMPMODE} = DBG ]; then
FFLAGS="${FFLAGS_DBG}"
CFLAGS="${CFLAGS_DBG}"
CXXFLAGS="${CXXFLAGS_DBG}"
fi
if [ $ARCH = Linux ]; then
FFLAGS="${FFLAGS_OPT} -Wno-globals"
CFLAGS="${CFLAGS_OPT}"
CXXFLAGS="${CXXFLAGS_OPT}"
FLIBS="-lfrtbegin -lg2c"
if [ ${COMPMODE} = OPT ]; then
FFLAGS="${FFLAGS_OPT}"
CFLAGS="${CFLAGS_OPT}"
CXXFLAGS="${CXXFLAGS_OPT}"
fi
if [ ${COMPMODE} = DBG ]; then
FFLAGS="${FFLAGS_DBG} -Wno-globals"
CFLAGS="${CFLAGS_DBG}"
CXXFLAGS="${CXXFLAGS_DBG}"
fi
fi
# Linux platform with gcc4: new Fortran90 compiler.
if [ $ARCH = Linux-gcc4 ]; then
FFLAGS="${FFLAGS_OPT}"
CFLAGS="${CFLAGS_OPT}"
CXXFLAGS="${CXXFLAGS_OPT}"
FLIBS="-lgfortran"
if [ ${COMPMODE} = OPT ]; then
FFLAGS="${FFLAGS_OPT}"
CFLAGS="${CFLAGS_OPT}"
CXXFLAGS="${CXXFLAGS_OPT}"
fi
if [ ${COMPMODE} = DBG ]; then
FFLAGS="${FFLAGS_DBG}"
CFLAGS="${CFLAGS_DBG}"
CXXFLAGS="${CXXFLAGS_DBG}"
fi
fi
if [ $ARCH = Linux-gcc5 ]; then
FFLAGS="${FFLAGS_OPT}"
CFLAGS="${CFLAGS_OPT}"
CXXFLAGS="${CXXFLAGS_OPT}"
FLIBS="-lgfortran"
if [ ${COMPMODE} = OPT ]; then
FFLAGS="${FFLAGS_OPT}"
CFLAGS="${CFLAGS_OPT}"
CXXFLAGS="${CXXFLAGS_OPT}"
fi
if [ ${COMPMODE} = DBG ]; then
FFLAGS="${FFLAGS_DBG}"
CFLAGS="${CFLAGS_DBG}"
CXXFLAGS="${CXXFLAGS_DBG}"
fi
fi
if [ $ARCH = Linux-gcc6 ]; then
FFLAGS="${FFLAGS_OPT}"
CFLAGS="${CFLAGS_OPT}"
CXXFLAGS="${CXXFLAGS_OPT}"
FLIBS="-lgfortran"
if [ ${COMPMODE} = OPT ]; then
FFLAGS="${FFLAGS_OPT}"
CFLAGS="${CFLAGS_OPT}"
CXXFLAGS="${CXXFLAGS_OPT}"
fi
if [ ${COMPMODE} = DBG ]; then
FFLAGS="${FFLAGS_DBG}"
CFLAGS="${CFLAGS_DBG}"
CXXFLAGS="${CXXFLAGS_DBG}"
fi
fi
# Add C++11 options if required
theGcc2=`g++ --version | grep '[0-9]\.[0-9]\.[0-9]' -o | head -1 | awk -F . '{print $2}'`
gccVar=`echo $theGcc`
gccVar2=`echo $theGcc2`
# Flag to make sure c++11 is enabled for compiling certain classes, e.g. MT random engine.
CPP11=1
# Check that we have gcc version 4.7 and above
if [ $gccVar -le 4 ]; then
if [ $gccVar2 -le 6 ]; then
echo "Not enabling c++11 features"
CPP11=0
fi
fi
if [ $CPP11 = 1 ]; then
echo "c++11 is enabled"
fi
# Mac OS platform with gcc4
if [[ $ARCH == Darwin* ]]; then
tt=`gfortran -print-search-dirs|grep libraries|cut -d'=' -f2|sed 's/:/ /g'`
LIBGFORTRANLOCATION=''
for i in $tt
do
if [ -e $i/libgfortran.so ]
then
LIBGFORTRANLOCATION=$i
break
elif [ -e $i/libgfortran.dylib ]
then
LIBGFORTRANLOCATION=$i
break
fi
done
FLIBS="-L${LIBGFORTRANLOCATION} -lgfortran"
SOFLAGS="-install_name"
fi
#Platform & opt/dbg - independent flags and variables:
echo -n "Creating config.mk ... "
rm -f config.mk
cat > config.mk << EOF
GENERATOR = ${GENERATOR}
VERSION = ${VERSION}
PREFIX = ${PREFIX}
SHELL = /bin/bash
ARCH = ${ARCH}
LCGPLATFORM = ${LCGPLATFORM}
FC = ${FC}
FFLAGS = ${FFLAGS}
CFLAGS = ${CFLAGS}
CXX = ${CXX}
CXXFLAGS = ${CXXFLAGS}
SOFLAGS = ${SOFLAGS}
EVTGENDIR = $(pwd)
TMPDIR = \$(EVTGENDIR)/tmp
INCLUDEDIR = \$(EVTGENDIR)
HEPMCDIR = ${HEPMCDIR}
HEPMCINCDIR = ${HEPMCDIR}/include
HEPMCLIBDIR = ${HEPMCDIR}/lib
HEPMCLIBLIST = ${HEPMCLIBLIST}
SRCDIR = \$(EVTGENDIR)/src
LIBDIR_SHARED = \$(EVTGENDIR)/lib
LIBDIR_ARCHIVE = \$(EVTGENDIR)/lib/archive
LIBDIRLIST = -lEvtGen
LIB_SHARED = \$(LIBDIR_SHARED)/lib\$(GENERATOR).so
LIB_ARCHIVE = \$(LIBDIR_ARCHIVE)/lib\$(GENERATOR).a
LIBEXT_SHARED = \$(LIBDIR_SHARED)/lib\$(GENERATOR)External.so
LIBEXT_ARCHIVE = \$(LIBDIR_ARCHIVE)/lib\$(GENERATOR)External.a
# Flags:
#
FFLAGSSHARED = ${FFLAGSSHARED}
CFLAGSSHARED = ${CFLAGSSHARED}
CXXFLAGSSHARED = ${CXXFLAGSSHARED}
LDFLAGSSHARED = ${LDFLAGSSHARED}
FLIBS = ${FLIBS}
# External generators:
EVTGEN_EXTERNAL = ${EVTGEN_EXTERNAL}
EOF
if [ "${EVTGEN_PYTHIA}" == "1" ] ; then
echo "EVTGEN_PYTHIA = 1" >> config.mk
echo "EXTRAFLAGS += -D EVTGEN_PYTHIA" >> config.mk
echo "PYTHIADIR = ${PYTHIADIR}" >> config.mk
echo "PYTHIAINCDIR = ${PYTHIADIR}/include" >> config.mk
echo "PYTHIALIBDIR = ${PYTHIADIR}/lib" >> config.mk
echo "PYTHIALIBLIST = ${PYTHIALIBLIST}" >> config.mk
fi
if [ "${EVTGEN_PHOTOS}" == "1" ] ; then
echo "EVTGEN_PHOTOS = 1" >> config.mk
echo "EXTRAFLAGS += -D EVTGEN_PHOTOS" >> config.mk
echo "PHOTOSDIR = ${PHOTOSDIR}" >> config.mk
echo "PHOTOSINCDIR = ${PHOTOSDIR}/include" >> config.mk
echo "PHOTOSLIBDIR = ${PHOTOSDIR}/lib" >> config.mk
echo "PHOTOSLIBLIST = ${PHOTOSLIBLIST}" >> config.mk
fi
if [ "${EVTGEN_TAUOLA}" == "1" ] ; then
echo "EVTGEN_TAUOLA = 1" >> config.mk
echo "EXTRAFLAGS += -D EVTGEN_TAUOLA" >> config.mk
echo "TAUOLADIR = ${TAUOLADIR}" >> config.mk
echo "TAUOLAINCDIR = ${TAUOLADIR}/include" >> config.mk
echo "TAUOLALIBDIR = ${TAUOLADIR}/lib" >> config.mk
echo "TAUOLALIBLIST = ${TAUOLALIBLIST}" >> config.mk
fi
if [ ${CPP11} == "1" ] ; then
echo "" >> config.mk
echo "# c++11 features:" >> config.mk
echo "EVTGEN_CPP11 = 1" >> config.mk
echo "CXXFLAGS += -std=c++11 -D EVTGEN_CPP11" >> config.mk
fi
echo " done"
set -v
mkdir -p tmp/ tmp/EvtGenBase tmp/EvtGenModels tmp_ext/EvtGenExternal lib/ lib/archive
set +v
echo -e "\n\nType 'make' to build everything or 'make help' to list build targets.\n\n"
exit
diff --git a/src/EvtGenExternal/EvtExternalGenFactory.cpp b/src/EvtGenExternal/EvtExternalGenFactory.cpp
index fc192e3..2dbc5ee 100644
--- a/src/EvtGenExternal/EvtExternalGenFactory.cpp
+++ b/src/EvtGenExternal/EvtExternalGenFactory.cpp
@@ -1,157 +1,174 @@
//--------------------------------------------------------------------------
//
// Environment:
// This software is part of the EvtGen package. If you use all or part
// of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
// Copyright (C) 2011 University of Warwick, UK
//
// Module: EvtExternalGenFactory
//
// Description: A factory type method to create engines for external physics
// generators like Pythia.
//
// Modification history:
//
// John Back April 2011 Module created
//
//------------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenExternal/EvtExternalGenFactory.hh"
#ifdef EVTGEN_PYTHIA
#include "EvtGenExternal/EvtPythiaEngine.hh"
#endif
#ifdef EVTGEN_PHOTOS
#include "EvtGenExternal/EvtPhotosEngine.hh"
#endif
#ifdef EVTGEN_TAUOLA
#include "EvtGenExternal/EvtTauolaEngine.hh"
#endif
#include <iostream>
using std::endl;
EvtExternalGenFactory::EvtExternalGenFactory() {
_extGenMap.clear();
}
EvtExternalGenFactory::~EvtExternalGenFactory() {
ExtGenMap::iterator iter;
for (iter = _extGenMap.begin(); iter != _extGenMap.end(); ++iter) {
EvtAbsExternalGen* theGenerator = iter->second;
delete theGenerator;
}
_extGenMap.clear();
}
EvtExternalGenFactory* EvtExternalGenFactory::getInstance() {
static EvtExternalGenFactory* theFactory = 0;
if (theFactory == 0) {
theFactory = new EvtExternalGenFactory();
}
return theFactory;
}
void EvtExternalGenFactory::definePythiaGenerator(std::string xmlDir,
bool convertPhysCodes,
bool useEvtGenRandom) {
// Only define the generator if we have the external ifdef variable set
#ifdef EVTGEN_PYTHIA
int genId = EvtExternalGenFactory::PythiaGenId;
- EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Defining EvtPythiaEngine: data tables defined in "
- <<xmlDir<<endl;
+ EvtAbsExternalGen* pythiaGenerator = _extGenMap[genId];
+ if (!pythiaGenerator) {
- if (convertPhysCodes == true) {
- EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Pythia 6 codes in decay files will be converted to Pythia 8 codes"<<endl;
- } else {
- EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Pythia 8 codes need to be used in decay files"<<endl;
- }
+ EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Defining EvtPythiaEngine: data tables defined in "
+ <<xmlDir<<endl;
- if (useEvtGenRandom == true) {
- EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Using EvtGen random engine for Pythia 8 as well"<<endl;
- }
+ if (convertPhysCodes == true) {
+ EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Pythia 6 codes in decay files will be converted to Pythia 8 codes"<<endl;
+ } else {
+ EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Pythia 8 codes need to be used in decay files"<<endl;
+ }
- EvtAbsExternalGen* pythiaGenerator = new EvtPythiaEngine(xmlDir, convertPhysCodes, useEvtGenRandom);
- _extGenMap[genId] = pythiaGenerator;
+ if (useEvtGenRandom == true) {
+ EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Using EvtGen random engine for Pythia 8 as well"<<endl;
+ }
+
+ pythiaGenerator = new EvtPythiaEngine(xmlDir, convertPhysCodes, useEvtGenRandom);
+ _extGenMap[genId] = pythiaGenerator;
+
+ }
#endif
}
void EvtExternalGenFactory::definePhotosGenerator(std::string photonType, bool useEvtGenRandom) {
#ifdef EVTGEN_PHOTOS
int genId = EvtExternalGenFactory::PhotosGenId;
- EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Defining EvtPhotosEngine using photonType = "<<photonType<<endl;
- EvtAbsExternalGen* photosGenerator = new EvtPhotosEngine(photonType, useEvtGenRandom);
- _extGenMap[genId] = photosGenerator;
+
+ EvtAbsExternalGen* photosGenerator = _extGenMap[genId];
+ if (!photosGenerator) {
+
+ EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Defining EvtPhotosEngine using photonType = "<<photonType<<endl;
+
+ photosGenerator = new EvtPhotosEngine(photonType, useEvtGenRandom);
+ _extGenMap[genId] = photosGenerator;
+ }
#endif
}
void EvtExternalGenFactory::defineTauolaGenerator(bool useEvtGenRandom) {
#ifdef EVTGEN_TAUOLA
int genId = EvtExternalGenFactory::TauolaGenId;
- EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Defining EvtTauolaEngine."<<endl;
- EvtAbsExternalGen* tauolaGenerator = new EvtTauolaEngine(useEvtGenRandom);
- _extGenMap[genId] = tauolaGenerator;
+
+ EvtAbsExternalGen* tauolaGenerator = _extGenMap[genId];
+ if (!tauolaGenerator) {
+
+ EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Defining EvtTauolaEngine."<<endl;
+
+ tauolaGenerator = new EvtTauolaEngine(useEvtGenRandom);
+ _extGenMap[genId] = tauolaGenerator;
+ }
#endif
}
EvtAbsExternalGen* EvtExternalGenFactory::getGenerator(int genId) {
EvtAbsExternalGen* theGenerator(0);
ExtGenMap::iterator iter;
if ((iter = _extGenMap.find(genId)) != _extGenMap.end()) {
// Retrieve the external generator engine
theGenerator = iter->second;
}
return theGenerator;
}
void EvtExternalGenFactory::initialiseAllGenerators() {
ExtGenMap::iterator iter;
for (iter = _extGenMap.begin(); iter != _extGenMap.end(); ++iter) {
EvtAbsExternalGen* theGenerator = iter->second;
if (theGenerator != 0) {
theGenerator->initialise();
}
}
}
diff --git a/test/evtgenlhc_test1.cc b/test/evtgenlhc_test1.cc
index 149acb9..eba0130 100644
--- a/test/evtgenlhc_test1.cc
+++ b/test/evtgenlhc_test1.cc
@@ -1,5714 +1,5714 @@
//#@# Dalitz plot for D0 --> K- pi+ pi0 decay:
//#@# 1: Mass(K-, pi+)
//#@# 2: Mass(pi+,pi0)
//
//--------------------------------------------------------------------------
//
// Environment:
// This software is part of the EvtGen package developed jointly
// for the BaBar and CLEO collaborations. If you use all or part
// of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
// Copyright (C) 1998 Caltech, UCSB
//
// Module: testEvtGen.cc
//
// Description:
//
// This program invokes the EvtGen event generator package
// for testing various decay models that are implemented.
//
// Modification history:
//
// DJL/RYD Sometime long ago Module created
//
//------------------------------------------------------------------------
//
//
//
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtVectorParticle.hh"
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtKine.hh"
#include "EvtGenBase/EvtGammaMatrix.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtRandomEngine.hh"
#include "EvtGenBase/EvtDecayTable.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtStdHep.hh"
#include "EvtGenBase/EvtSecondary.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGen/EvtGen.hh"
#include "EvtGenBase/EvtParticleFactory.hh"
#include "EvtGenBase/EvtSimpleRandomEngine.hh"
#include "EvtGenBase/EvtMTRandomEngine.hh"
#include "EvtGenBase/EvtIdSet.hh"
#include "EvtGenBase/EvtParser.hh"
#include "EvtGenBase/EvtAbsRadCorr.hh"
#include "EvtGenBase/EvtDecayBase.hh"
#ifdef EVTGEN_EXTERNAL
#include "EvtGenExternal/EvtExternalGenList.hh"
#endif
#include <cstdio>
#include <fstream>
#include <sstream>
#include <cmath>
#include <string>
#include <vector>
#include <cstdlib>
#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TApplication.h"
#include "TROOT.h"
using std::vector;
void runFile(int nevent,char* fname,EvtGen& myGenerator);
void runPrint(int nevent,char* fname,EvtGen& myGenerator);
void runFileVpho(int nevent,char* fname,EvtGen& myGenerator);
void runTest1(int nevent,EvtGen& myGenerator);
void runTest2(int nevent,EvtGen& myGenerator);
void runOmega(int nevent,EvtGen& myGenerator);
void runChi1Kstar(int nevent,EvtGen& myGenerator);
void runPi0Dalitz(int nevent,EvtGen& myGenerator);
void runMix(int nevent,EvtGen& myGenerator);
void runBMix(int nevent,EvtGen& myGenerator,std::string userFile,std::string rootFile);
void runDDalitz(int nevent,EvtGen& myGenerator);
void runPiPiCPT(int nevent,EvtGen& myGenerator);
void runPiPiPiPi(int nevent,EvtGen& myGenerator);
void runD2Pi(int nevent,EvtGen& myGenerator);
void runJetsetTab3(int nevent,EvtGen& myGenerator);
void runHelAmp(int nevent,EvtGen& myGenerator,std::string userFile,std::string rootFile);
void runHelAmp2(int nevent,EvtGen& myGenerator);
void runJpsiKs(int nevent,EvtGen& myGenerator);
void runDump(int nevent,EvtGen& myGenerator);
void runD1(int nevent,EvtGen& myGenerator);
void runGenericCont(int nevent,EvtGen& myGenerator);
void runPiPiPi(int nevent,EvtGen& myGenerator);
void runBHadronic(int nevent,EvtGen& myGenerator);
void runSingleB(int nevent,EvtGen& myGenerator);
void runA2Pi(int nevent,EvtGen& myGenerator);
void runAlias();
void runRepeat(int nevent);
void runPhotos(int nevent,EvtGen& myGenerator);
void runTrackMult(int nevent,EvtGen& myGenerator);
void runGeneric(int neventOrig,EvtGen& myGenerator,
std::string listfile);
void runFinalStates(int nevent,EvtGen& myGenerator);
std::vector<std::string> findFinalState(EvtParticle *p);
void runKstarnunu(int nevent,EvtGen& myGenerator);
void runBsmix(int nevent,EvtGen& myGenerator);
void runTauTauPiPi(int nevent,EvtGen& myGenerator);
void runTauTauEE(int nevent,EvtGen& myGenerator);
void runTauTau2Pi2Pi(int nevent,EvtGen& myGenerator);
void runTauTau3Pi3Pi(int nevent,EvtGen& myGenerator);
void runJPsiKstar(int nevent,EvtGen& myGenerator,int modeInt);
void runSVVCPLH(int nevent,EvtGen& myGenerator);
void runSVSCPLH(int nevent,EvtGen& myGenerator);
void runSSDCP(int nevent,EvtGen& myGenerator);
void runKstarstargamma(int nevent,EvtGen& myGenerator);
void runDSTARPI(int nevent,EvtGen& myGenerator);
void runETACPHIPHI(int nevent,EvtGen& myGenerator);
void runVVPiPi(int nevent,EvtGen& myGenerator);
void runSVVHelAmp(int nevent,EvtGen& myGenerator);
void runSVVHelAmp2(int nevent,EvtGen& myGenerator);
void runPartWave(int nevent,EvtGen& myGenerator);
void runPartWave2(int nevent,EvtGen& myGenerator);
void runTwoBody(int nevent,EvtGen& myGenerator,std::string decfile,std::string rootFile);
void runPiPi(int nevent,EvtGen& myGenerator);
void runA1Pi(int nevent,EvtGen& myGenerator);
void runCPTest(int nevent,EvtGen& myGenerator);
void runSemic(int nevent,EvtGen& myGenerator);
void runKstarll(int nevent,EvtGen& myGenerator);
void runKll(int nevent,EvtGen& myGenerator);
void runHll(int nevent,EvtGen& myGenerator, char* mode);
void runVectorIsr(int nevent,EvtGen& myGenerator);
void runBsquark(int nevent,EvtGen& myGenerator);
void runK3gamma(int nevent,EvtGen& myGenerator);
void runLambda(int nevent,EvtGen& myGenerator);
void runBtoXsgamma(int nevent,EvtGen& myGenerator);
void runBtoK1273gamma(int nevent,EvtGen& myGenerator);
void runCheckRotBoost();
void runMassCheck(int nevent, EvtGen& myGenerator, int partnum);
void runJpsiPolarization(int nevent, EvtGen& myGenerator);
void runDDK(int nevent, EvtGen &myGenerator);
int countInclusive(std::string name, EvtParticle *root,TH1F* mom=0,TH1F* mass=0 );
int countInclusiveParent(std::string name, EvtParticle *root,EvtIdSet setIds,
TH1F* mom=0);
int countInclusiveSubTree(std::string name, EvtParticle *root,EvtIdSet setIds,
TH1F* mom=0);
void runBaryonic(int nEvent, EvtGen& myGenerator);
int main(int argc, char* argv[]){
// Define the random number generator
EvtRandomEngine* myRandomEngine = 0;
#ifdef EVTGEN_CPP11
// Use the Mersenne-Twister generator (C++11 only)
myRandomEngine = new EvtMTRandomEngine();
#else
myRandomEngine = new EvtSimpleRandomEngine();
#endif
if (!TROOT::Initialized()) {
static TROOT root("RooTuple", "RooTuple ROOT in EvtGen");
}
if (argc==1){
EvtVector4R p(0.0,1.0,0.0,0.0);
EvtVector4R k(0.0,0.0,1.0,0.0);
EvtTensor4C T=dual(EvtGenFunctions::directProd(p,k));
EvtGenReport(EVTGEN_INFO,"EvtGen") << "p:"<<p<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "k:"<<k<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "T=dual(directProd(p,k)):"<<T<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "T03:"<<T.get(0,3)<<std::endl;
return 1;
}
EvtAbsRadCorr* radCorrEngine = 0;
std::list<EvtDecayBase*> extraModels;
#ifdef EVTGEN_EXTERNAL
EvtExternalGenList genList;
radCorrEngine = genList.getPhotosModel();
extraModels = genList.getListOfModels();
#endif
EvtGen myGenerator("../DECAY_2010.DEC", "../evt.pdl", myRandomEngine,
radCorrEngine, &extraModels);
if (!strcmp(argv[1],"file")) {
int nevent=atoi(argv[2]);
runFile(nevent,argv[3],myGenerator);
}
if (!strcmp(argv[1],"print")) {
int nevent=atoi(argv[2]);
runPrint(nevent,argv[3],myGenerator);
}
if (!strcmp(argv[1],"filevpho")) {
int nevent=atoi(argv[2]);
runFileVpho(nevent,argv[3],myGenerator);
}
if (!strcmp(argv[1],"test1")) {
int nevent=atoi(argv[2]);
runTest1(nevent,myGenerator);
}
if (!strcmp(argv[1],"chi1kstar")) {
int nevent=atoi(argv[2]);
runChi1Kstar(nevent,myGenerator);
}
if (!strcmp(argv[1],"test2")) {
int nevent=atoi(argv[2]);
runTest2(nevent,myGenerator);
}
if (!strcmp(argv[1],"omega")) {
int nevent=atoi(argv[2]);
runOmega(nevent,myGenerator);
}
if (!strcmp(argv[1],"alias")) {
runAlias();
}
if (!strcmp(argv[1],"repeat")) {
int nevent=atoi(argv[2]);
runRepeat(nevent);
}
if (!strcmp(argv[1],"photos")) {
int nevent=atoi(argv[2]);
runPhotos(nevent,myGenerator);
}
if (!strcmp(argv[1],"trackmult")) {
int nevent=atoi(argv[2]);
runTrackMult(nevent,myGenerator);
}
if (!strcmp(argv[1],"generic")) {
int nevent=atoi(argv[2]);
std::string listfile("");
if ( argc==4) listfile=argv[3];
runGeneric(nevent,myGenerator,listfile);
}
if (!strcmp(argv[1],"finalstates")) {
int nevent=atoi(argv[2]);
runFinalStates(nevent,myGenerator);
}
if (!strcmp(argv[1],"kstarnunu")) {
int nevent=atoi(argv[2]);
runKstarnunu(nevent,myGenerator);
}
if (!strcmp(argv[1],"bsmix")) {
int nevent=atoi(argv[2]);
runBsmix(nevent,myGenerator);
}
if (!strcmp(argv[1],"BtoXsgamma")) {
int nevent=atoi(argv[2]);
runBtoXsgamma(nevent,myGenerator);
}
if (!strcmp(argv[1],"BtoK1273gamma")) {
int nevent=atoi(argv[2]);
runBtoK1273gamma(nevent,myGenerator);
}
if (!strcmp(argv[1],"pi0dalitz")) {
int nevent=atoi(argv[2]);
runPi0Dalitz(nevent,myGenerator);
}
if (!strcmp(argv[1],"ddalitz")) {
int nevent=atoi(argv[2]);
runDDalitz(nevent,myGenerator);
}
if (!strcmp(argv[1],"kstarll")) {
int nevent=atoi(argv[2]);
runKstarll(nevent, myGenerator);
}
if (!strcmp(argv[1],"kll")) {
int nevent=atoi(argv[2]);
runKll(nevent, myGenerator);
}
if (!strcmp(argv[1],"hll")) {
int nevent=atoi(argv[2]);
runHll(nevent, myGenerator, argv[3]);
}
if (!strcmp(argv[1],"vectorisr")) {
int nevent=atoi(argv[2]);
runVectorIsr(nevent, myGenerator);
}
if (!strcmp(argv[1],"bsquark")) {
int nevent=atoi(argv[2]);
runBsquark(nevent, myGenerator);
}
if (!strcmp(argv[1],"k3gamma")) {
int nevent=atoi(argv[2]);
runK3gamma(nevent, myGenerator);
}
if (!strcmp(argv[1],"lambda")) {
int nevent=atoi(argv[2]);
runLambda(nevent, myGenerator);
}
if (!strcmp(argv[1],"tautaupipi")) {
int nevent=atoi(argv[2]);
runTauTauPiPi(nevent, myGenerator);
}
if (!strcmp(argv[1],"tautauee")) {
int nevent=atoi(argv[2]);
runTauTauEE(nevent, myGenerator);
}
if (!strcmp(argv[1],"tautau2pi2pi")) {
int nevent=atoi(argv[2]);
runTauTau2Pi2Pi(nevent, myGenerator);
}
if (!strcmp(argv[1],"tautau3pi3pi")) {
int nevent=atoi(argv[2]);
runTauTau3Pi3Pi(nevent, myGenerator);
}
if (!strcmp(argv[1],"jpsikstar")) {
int nevent=atoi(argv[2]);
int modeInt=atoi(argv[3]);
runJPsiKstar(nevent, myGenerator,modeInt);
}
if (!strcmp(argv[1],"svvcplh")) {
int nevent=atoi(argv[2]);
runSVVCPLH(nevent, myGenerator);
}
if (!strcmp(argv[1],"svscplh")) {
int nevent=atoi(argv[2]);
runSVSCPLH(nevent, myGenerator);
}
if (!strcmp(argv[1],"ssdcp")) {
int nevent=atoi(argv[2]);
runSSDCP(nevent, myGenerator);
}
if (!strcmp(argv[1],"kstarstargamma")) {
int nevent=atoi(argv[2]);
runKstarstargamma(nevent, myGenerator);
}
if (!strcmp(argv[1],"dstarpi")) {
int nevent=atoi(argv[2]);
runDSTARPI(nevent, myGenerator);
}
if (!strcmp(argv[1],"etacphiphi")) {
int nevent=atoi(argv[2]);
runETACPHIPHI(nevent, myGenerator);
}
if (!strcmp(argv[1],"vvpipi")) {
int nevent=atoi(argv[2]);
runVVPiPi(nevent, myGenerator);
}
if (!strcmp(argv[1],"svvhelamp")) {
int nevent=atoi(argv[2]);
runSVVHelAmp(nevent, myGenerator);
}
if (!strcmp(argv[1],"partwave")) {
int nevent=atoi(argv[2]);
runPartWave(nevent, myGenerator);
}
if (!strcmp(argv[1],"partwave2")) {
int nevent=atoi(argv[2]);
runPartWave2(nevent, myGenerator);
}
if (!strcmp(argv[1],"twobody")) {
int nevent=atoi(argv[2]);
runTwoBody(nevent, myGenerator, argv[3], argv[4]);
}
if (!strcmp(argv[1],"pipipi")) {
int nevent=atoi(argv[2]);
runPiPiPi(nevent, myGenerator);
}
if (!strcmp(argv[1],"bhadronic")) {
int nevent=atoi(argv[2]);
runBHadronic(nevent, myGenerator);
}
if (!strcmp(argv[1],"singleb")) {
int nevent=atoi(argv[2]);
runSingleB(nevent, myGenerator);
}
if (!strcmp(argv[1],"pipi")) {
int nevent=atoi(argv[2]);
runPiPi(nevent, myGenerator);
}
if (!strcmp(argv[1],"pipipipi")) {
int nevent=atoi(argv[2]);
runPiPiPiPi(nevent, myGenerator);
}
if (!strcmp(argv[1],"a2pi")) {
int nevent=atoi(argv[2]);
runA2Pi(nevent, myGenerator);
}
if (!strcmp(argv[1],"helamp")) {
int nevent=atoi(argv[2]);
runHelAmp(nevent, myGenerator, argv[3], argv[4]);
}
if (!strcmp(argv[1],"helamp2")) {
int nevent=atoi(argv[2]);
runHelAmp2(nevent, myGenerator);
}
if (!strcmp(argv[1],"d2pi")) {
int nevent=atoi(argv[2]);
runD2Pi(nevent, myGenerator);
}
if (!strcmp(argv[1],"a1pi")) {
int nevent=atoi(argv[2]);
runA1Pi(nevent, myGenerator);
}
if (!strcmp(argv[1],"cptest")) {
int nevent=atoi(argv[2]);
runCPTest(nevent, myGenerator);
}
if (!strcmp(argv[1],"pipicpt")) {
int nevent=atoi(argv[2]);
runPiPiCPT(nevent, myGenerator);
}
if (!strcmp(argv[1],"jpsiks")) {
int nevent=atoi(argv[2]);
runJpsiKs(nevent, myGenerator);
}
if (!strcmp(argv[1],"dump")) {
int nevent=atoi(argv[2]);
runDump(nevent, myGenerator);
}
if (!strcmp(argv[1],"genericcont")) {
int nevent=atoi(argv[2]);
runGenericCont(nevent, myGenerator);
}
if (!strcmp(argv[1],"d1")) {
int nevent=atoi(argv[2]);
runD1(nevent, myGenerator);
}
if (!strcmp(argv[1],"mix")) {
int nevent=atoi(argv[2]);
runMix(nevent, myGenerator);
}
if (!strcmp(argv[1],"bmix")) {
int nevent=atoi(argv[2]);
runBMix(nevent, myGenerator, argv[3], argv[4]);
}
if (!strcmp(argv[1],"semic")) {
int nevent=atoi(argv[2]);
runSemic(nevent,myGenerator);
}
if (!strcmp(argv[1],"ddk")) {
int nevent=atoi(argv[2]);
runDDK(nevent,myGenerator);
}
if (!strcmp(argv[1],"checkmass")) {
int nevent=atoi(argv[2]);
int partnum=atoi(argv[3]);
runMassCheck(nevent,myGenerator,partnum);
}
if (!strcmp(argv[1],"jpsipolarization")) {
int nevent=atoi(argv[2]);
runJpsiPolarization(nevent,myGenerator);
}
//*******************************************************
//test of the rotations and boosts performed in EvtGen.
// Added by Lange and Ryd Jan 5,2000.
//*******************************************************
if (!strcmp(argv[1],"checkrotboost")) {
runCheckRotBoost();
}
if ( !strcmp(argv[1], "baryonic") )
{
runBaryonic( atoi(argv[2]), myGenerator );
}
delete myRandomEngine;
return 0;
}
void runFile(int nevent,char* fname, EvtGen &myGenerator) {
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
int count;
char udecay_name[100];
strcpy(udecay_name,fname);
myGenerator.readUDecay(udecay_name);
count=1;
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
root_part->deleteTree();
}while (count++<nevent);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runPrint(int nevent,char* fname, EvtGen &myGenerator) {
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
int count;
char udecay_name[100];
strcpy(udecay_name,fname);
myGenerator.readUDecay(udecay_name);
count=1;
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
root_part->printTree();
root_part->deleteTree();
}while (count++<nevent);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runFileVpho(int nevent,char* fname, EvtGen &myGenerator) {
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
int count;
char udecay_name[100];
strcpy(udecay_name,fname);
myGenerator.readUDecay(udecay_name);
count=1;
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(VPHO,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
root_part->deleteTree();
}while (count++<nevent);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
///////////////
void runJpsiPolarization(int nevent, EvtGen &myGenerator) {
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId JPSI=EvtPDL::getId(std::string("J/psi"));
int count;
myGenerator.readUDecay("exampleFiles/GENERIC.DEC");
myGenerator.readUDecay("exampleFiles/JPSITOLL.DEC");
TFile *file=new TFile("jpsipolar.root", "RECREATE");
TH1F* coshel = new TH1F("h1","cos hel",50,-1.0,1.0);
TH1F* coshelHigh = new TH1F("h2","cos hel pstar gt 1.1",50,-1.0,1.0);
TH1F* coshelLow = new TH1F("h3","cos hel pstar lt 1.1",50,-1.0,1.0);
count=1;
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p = root_part;
do{
EvtId type=p->getId();
if (p->getId()==JPSI) {
EvtVector4R p4psi=p->getP4Lab();
EvtVector4R p4Daug=p->getDaug(0)->getP4Lab();
double dcostheta=EvtDecayAngle(p_init,p4psi,p4Daug);
coshel->Fill(dcostheta);
if ( p4psi.d3mag() > 1.1 ) {
coshelHigh->Fill(dcostheta);
}
else{
coshelLow->Fill(dcostheta);
}
}
p=p->nextIter(root_part);
}while(p!=0);
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runMassCheck(int nevent, EvtGen &myGenerator,
int partnum) {
int count;
static EvtId myPart=EvtPDL::evtIdFromStdHep(partnum);
TFile *file=new TFile("checkmass.root", "RECREATE");
TH1F* mass = new TH1F("h1","Mass",500,0.0,2.5);
count=1;
do{
mass->Fill(EvtPDL::getMass(myPart));
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runPi0Dalitz(int nevent, EvtGen &myGenerator) {
static EvtId PI0=EvtPDL::getId(std::string("pi0"));
TFile *file=new TFile("pi0dalitz.root", "RECREATE");
TH1F* q2 = new TH1F("h1","q2",50,0.0,0.02);
int count;
myGenerator.readUDecay("exampleFiles/PI0DALITZ.DEC");
count=1;
do{
EvtVector4R p_init(EvtPDL::getMass(PI0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(PI0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
EvtVector4R ep=root_part->getDaug(0)->getP4Lab();
EvtVector4R em=root_part->getDaug(1)->getP4Lab();
EvtVector4R gamma=root_part->getDaug(2)->getP4Lab();
q2->Fill( (ep+em).mass2() );
// EvtGenReport(EVTGEN_INFO,"EvtGen") << ep << em << gamma <<std::endl;
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
//*******************************************************************************
void runTest1(int nevent, EvtGen &myGenerator) {
// TFile *file=new TFile("test1.root", "RECREATE");
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
// int first=0;
// char **second;
// TApplication *theApp = new TApplication("App", &first, second);
TFile *file = new TFile("test1.root", "RECREATE", "Example");
TH1F * costhetaB = new TH1F("hcosthetaB","costhetaB",50,-1.0,1.0);
TH1F* phiB = new TH1F("hphiB","phiB",50,-EvtConst::pi,EvtConst::pi);
TH1F* Elep = new TH1F("hEl","E?l!",50,0.0,2.5);
TH1F* q2 = new TH1F("hq2","q^2!",44,0.0,11.0);
TH1F* ctv = new TH1F("hctv","ctv",50,-1.0,1.0);
TH1F* chi_low_ctv = new TH1F("hcostv1","[h] for cos[Q]?V!\"L#0",50,0.0,EvtConst::twoPi);
TH1F* chi_high_ctv = new TH1F("hcostv2","[h] for cos[Q]?V!\"G#0",50,0.0,EvtConst::twoPi);
TH1F* dt = new TH1F("hdt","dt",50,-5.0,5.0);
int count;
EvtVector4R p4b0,p4b0b,p4dstar,p4e,p4nu,p4d,p4pi,p4pip,p4pim;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/TEST1.DEC");
// EvtGen myGenerator(decay_name,pdttable_name,myRandomEngine);
myGenerator.readUDecay(udecay_name);
double costhetaV;
count=1;
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
p4b0=root_part->getDaug(0)->getP4Lab();
p4b0b=root_part->getDaug(1)->getP4Lab();
p4dstar=root_part->getDaug(0)->getDaug(0)->getP4Lab();
p4e=root_part->getDaug(0)->getDaug(1)->getP4Lab();
p4nu=root_part->getDaug(0)->getDaug(2)->getP4Lab();
p4d=root_part->getDaug(0)->getDaug(0)->getDaug(0)->getP4Lab();
p4pi=root_part->getDaug(0)->getDaug(0)->getDaug(1)->getP4Lab();
p4pip=root_part->getDaug(1)->getDaug(0)->getP4Lab();
p4pim=root_part->getDaug(1)->getDaug(1)->getP4Lab();
costhetaB->Fill(p4b0.get(3)/p4b0.d3mag());
phiB->Fill(atan2(p4b0.get(1),p4b0.get(2)));
Elep->Fill(p4b0*p4e/p4b0.mass());
q2->Fill((p4e+p4nu).mass2());
dt->Fill(root_part->getDaug(1)->getLifetime()-root_part->getDaug(0)->getLifetime());
costhetaV=EvtDecayAngle(p4b0,p4d+p4pi,p4d);
ctv->Fill(costhetaV);
if (costhetaV<0.0){
chi_low_ctv->Fill(EvtDecayAngleChi(p4b0,p4d,p4pi,p4e,p4nu));
}
else{
chi_high_ctv->Fill(EvtDecayAngleChi(p4b0,p4d,p4pi,p4e,p4nu));
}
root_part->deleteTree();
}while (count++<nevent);
file->Write();
file->Close();
// delete theApp;
// hfile.write();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
//*******************************************************************************
void runDDK(int nevent, EvtGen &myGenerator) {
// TFile *file=new TFile("test1.root", "RECREATE");
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
int count;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/GENERIC.DEC");
myGenerator.readUDecay(udecay_name);
count=1;
static EvtId kp=EvtPDL::getId(std::string("K+"));
static EvtId km=EvtPDL::getId(std::string("K-"));
static EvtId ks=EvtPDL::getId(std::string("K_S0"));
static EvtId kl=EvtPDL::getId(std::string("K_L0"));
static EvtId k0=EvtPDL::getId(std::string("K0"));
static EvtId kb=EvtPDL::getId(std::string("anti-K0"));
static EvtId d0=EvtPDL::getId(std::string("D0"));
static EvtId dp=EvtPDL::getId(std::string("D+"));
static EvtId dm=EvtPDL::getId(std::string("D-"));
static EvtId db=EvtPDL::getId(std::string("anti-D0"));
static EvtIdSet theKs(kp,km,ks,kl,k0,kb);
static EvtIdSet theDs(d0,dp,dm,db);
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
static EvtId BP=EvtPDL::getId(std::string("B+"));
static EvtId BM=EvtPDL::getId(std::string("B-"));
static EvtIdSet theBs(B0B,B0,BP,BM);
int nDDK=0;
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *theB01=root_part->getDaug(0);
EvtParticle *theB02=root_part->getDaug(1);
int nD=0;
int nK=0;
EvtParticle *p=theB01;
do{
EvtId type=p->getId();
EvtId typePar=p->getParent()->getId();
if (theDs.contains(type)) nD++;
if (theKs.contains(type) && theBs.contains(typePar)) nK++;
p=p->nextIter(theB01);
}while(p!=0);
if ( nD==2 && nK==1 ) nDDK++;
nD=0;
nK=0;
p=theB02;
do{
EvtId type=p->getId();
EvtId typePar=p->getParent()->getId();
if (theDs.contains(type)) nD++;
if (theKs.contains(type) && theBs.contains(typePar)) nK++;
p=p->nextIter(theB02);
}while(p!=0);
if ( nD==2 && nK==1 ) nDDK++;
root_part->deleteTree();
}while (count++<nevent);
EvtGenReport(EVTGEN_INFO,"EvtGen") << nDDK << " " << (count-1) << " " << nDDK/float(2*(count-1)) << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
//*******************************************************************************
void runTest2(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("test2.root", "RECREATE");
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
TH1F* costhetaB = new TH1F("h1","cos[Q]?B!",50,-1.0,1.0);
TH1F* phiB = new TH1F("h2","[f]?B!",50,-EvtConst::pi,
EvtConst::pi);
TH1F* dt = new TH1F("h3","[D]t",100,-5.0,5.0);
TH1F* costhetaJpsi = new TH1F("h4","cos[Q]?J/[y]!",
50,-1.0,1.0);
TH1F* costhetaKstar = new TH1F("h5","cos[Q]?K*!",
50,-1.0,1.0);
TH1F* chi = new TH1F("h6","[h]",50,0.0,EvtConst::twoPi);
TH1F* chi1 = new TH1F("h26","[h] [D]t\"L#0",50,0.0,
EvtConst::twoPi);
TH1F* chi2 = new TH1F("h27","[h] [D]t\"G#0",50,0.0,
EvtConst::twoPi);
TH1F* costhetaomega = new TH1F("h7","costhetaomega",50,-1.0,1.0);
TH1F* costhetaomega1 = new TH1F("h20","costhetaomega1",50,-1.0,1.0);
TH1F* costhetaomega2 = new TH1F("h21","costhetaomega2",50,-1.0,1.0);
TH1F* costhetaomega3 = new TH1F("h22","costhetaomega3",50,-1.0,1.0);
TH1F* omegaamp = new TH1F("h8","omegaamp",50,0.0,0.05);
TH1F* omegaamp1 = new TH1F("h9","omegaamp1",50,0.0,0.05);
TH1F* omegaamp2 = new TH1F("h10","omegaamp2",50,0.0,0.05);
TH1F* omegaamp3 = new TH1F("h11","omegaamp3",50,0.0,0.05);
TH2F* chi1vscoskstarl =
new TH2F("h30","[h] vs. cos[Q]?J/[y]! [D]t\"L#0",
20,0.0,EvtConst::twoPi,20,-1.0,1.0);
TH2F* chi1vscoskstarg =
new TH2F("h31","[h] vs. cos[Q]?J/[y]! [D]t\"G#0",
20,0.0,EvtConst::twoPi,20,-1.0,1.0);
int count=1;
EvtVector4R p4_b0,p4_b0b,p4_psi,p4_kstar,p4_mup,p4_mum;
EvtVector4R p4_kz,p4_pi0,p4_pi1,p4_pi2,p4_pi3,p4_omega;
EvtVector4R p4_pi1_omega,p4_pi2_omega,p4_pi3_omega;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/TEST2.DEC");
//EvtGen myGenerator(decay_name,pdttable_name,myRandomEngine);
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
p4_b0=root_part->getDaug(0)->getP4Lab();
p4_b0b=root_part->getDaug(1)->getP4Lab();
p4_psi=root_part->getDaug(1)->getDaug(0)->getP4Lab();
p4_kstar=root_part->getDaug(1)->getDaug(1)->getP4Lab();
p4_mup=root_part->getDaug(1)->getDaug(0)->getDaug(0)->getP4Lab();
p4_mum=root_part->getDaug(1)->getDaug(0)->getDaug(1)->getP4Lab();
p4_kz=root_part->getDaug(1)->getDaug(1)->getDaug(0)->getP4Lab();
p4_pi0=root_part->getDaug(1)->getDaug(1)->getDaug(1)->getP4Lab();
p4_omega=root_part->getDaug(0)->getDaug(0)->getP4Lab();
p4_pi1=root_part->getDaug(0)->getDaug(0)->getDaug(0)->getP4Lab();
p4_pi2=root_part->getDaug(0)->getDaug(0)->getDaug(1)->getP4Lab();
p4_pi3=root_part->getDaug(0)->getDaug(0)->getDaug(2)->getP4Lab();
//get momentum in the omega restframe
p4_pi1_omega=root_part->getDaug(0)->getDaug(0)->getDaug(0)->getP4();
p4_pi2_omega=root_part->getDaug(0)->getDaug(0)->getDaug(1)->getP4();
p4_pi3_omega=root_part->getDaug(0)->getDaug(0)->getDaug(2)->getP4();
EvtVector3R p3_perp=cross(EvtVector3R(p4_pi2_omega.get(0),
p4_pi2_omega.get(1),
p4_pi2_omega.get(2)),
EvtVector3R(p4_pi3_omega.get(0),
p4_pi3_omega.get(1),
p4_pi3_omega.get(2)));
EvtVector4R p4_perp(p3_perp.d3mag(),p3_perp.get(0),
p3_perp.get(1),p3_perp.get(2));
root_part->getDaug(0)->getDaug(0)->getDaug(0)->setP4(p4_perp);
p4_perp=root_part->getDaug(0)->getDaug(0)->getDaug(0)->getP4Lab();
EvtVector4R p4_perpprime=p4_omega-p4_perp;
double d_omegaamp=EvtVector3R(p4_pi1_omega.get(0),
p4_pi1_omega.get(1),
p4_pi1_omega.get(2))*p3_perp;
d_omegaamp*=d_omegaamp;
d_omegaamp*=20.0;
double d_dt=root_part->getDaug(1)->getLifetime()-
root_part->getDaug(0)->getLifetime();
double d_costhetaJpsi=EvtDecayAngle(p4_b0b,p4_mup+p4_mum,p4_mup);
double d_costhetaKstar=EvtDecayAngle(p4_b0b,p4_pi0+p4_kz,p4_pi0);
double d_chi=EvtDecayAngleChi(p4_b0b,p4_pi0,p4_kz,p4_mup, p4_mum );
costhetaB->Fill(p4_b0.get(3)/p4_b0.d3mag());
phiB->Fill(atan2(p4_b0.get(1),p4_b0.get(2)));
dt->Fill(d_dt);
costhetaJpsi->Fill(d_costhetaJpsi);
costhetaKstar->Fill(d_costhetaKstar);
chi->Fill(d_chi);
if (d_dt<0.0) {
chi1->Fill(d_chi);
chi1vscoskstarl->Fill(d_chi,d_costhetaJpsi,1.0);
}
if (d_dt>0.0) {
chi2->Fill(d_chi);
chi1vscoskstarg->Fill(d_chi,d_costhetaJpsi,1.0);
}
double d_costhetaomega=EvtDecayAngle(p4_b0b,p4_perp+p4_perpprime,p4_perp);
costhetaomega->Fill(d_costhetaomega);
if (d_omegaamp<0.001) costhetaomega1->Fill(d_costhetaomega);
if (d_omegaamp>0.02) costhetaomega2->Fill(d_costhetaomega);
if (std::fabs(d_omegaamp-0.015)<0.005) costhetaomega3->Fill(d_costhetaomega);
omegaamp->Fill(d_omegaamp);
if (d_costhetaomega<-0.5) omegaamp1->Fill(d_omegaamp);
if (d_costhetaomega>0.5) omegaamp2->Fill(d_omegaamp);
if (std::fabs(d_costhetaomega)<0.5) omegaamp3->Fill(d_omegaamp);
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runOmega(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("omega.root", "RECREATE");
static EvtId OMEGA=EvtPDL::getId(std::string("omega"));
TH2F* dalitz = new TH2F("h1","E1 vs E2",50,0.0,0.5,50,0.0,0.5);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/OMEGA.DEC");
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMeanMass(OMEGA),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(OMEGA,
p_init);
//root_part->setDiagonalSpinDensity();
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
dalitz->Fill(root_part->getDaug(0)->getP4().get(0),
root_part->getDaug(1)->getP4().get(0),1.0);
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runChi1Kstar(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("chi1kstar.root", "RECREATE");
static EvtId B0=EvtPDL::getId(std::string("B0"));
TH1F* costhetaChi1 = new TH1F("h1","cos[Q]?J/[x]!",
50,-1.0,1.0);
TH1F* costhetaKstar = new TH1F("h2","cos[Q]?K*!",
50,-1.0,1.0);
TH1F* chi = new TH1F("h3","[x]",50,-EvtConst::pi,EvtConst::pi);
int count=1;
EvtVector4R p4_b,p4_chi,p4_kstar,p4_gamma,p4_psi,p4_k,p4_p;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/CHI1KSTAR.DEC");
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(B0,
p_init);
myGenerator.generateDecay(root_part);
p4_b=root_part->getP4Lab();
p4_chi=root_part->getDaug(0)->getP4Lab();
p4_kstar=root_part->getDaug(1)->getP4Lab();
p4_psi=root_part->getDaug(0)->getDaug(0)->getP4Lab();
p4_gamma=root_part->getDaug(0)->getDaug(1)->getP4Lab();
p4_k=root_part->getDaug(1)->getDaug(0)->getP4Lab();
p4_p=root_part->getDaug(1)->getDaug(1)->getP4Lab();
double d_costhetaChi1=EvtDecayAngle(p4_b,p4_chi,p4_psi);
double d_costhetaKstar=EvtDecayAngle(p4_b,p4_kstar,p4_k);
double d_chi=EvtDecayAngleChi(p4_b,p4_k,p4_p,p4_psi, p4_gamma );
if (d_chi>EvtConst::pi) d_chi-=EvtConst::twoPi;
costhetaChi1->Fill(d_costhetaChi1);
costhetaKstar->Fill(d_costhetaKstar);
chi->Fill(d_chi);
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runAlias() {
EvtId idpip=EvtPDL::getId(std::string("pi+"));
EvtPDL::alias(idpip,std::string("my_pi+"));
EvtId myidpip=EvtPDL::getId(std::string("my_pi+"));
EvtId idpim=EvtPDL::getId(std::string("pi-"));
EvtPDL::alias(idpim,std::string("my_pi-"));
EvtId myidpim=EvtPDL::getId(std::string("my_pi-"));
EvtId idpi0=EvtPDL::getId(std::string("pi0"));
EvtPDL::alias(idpi0,std::string("my_pi0"));
EvtId myidpi0=EvtPDL::getId(std::string("my_pi0"));
EvtPDL::aliasChgConj(myidpip,myidpim);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Id pi+:" << idpip << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Id pi-:" << idpim << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Id pi0:" << idpi0 << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Id my_pi+:" << myidpip << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Id my_pi-:" << myidpim << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Id my_pi0:" << myidpi0 << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Chg conj pi+:" <<
EvtPDL::chargeConj(idpip) << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Chg conj pi-:" <<
EvtPDL::chargeConj(idpim) << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Chg conj pi0:" <<
EvtPDL::chargeConj(idpi0) << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Chg conj my_pi+:" <<
EvtPDL::chargeConj(myidpip) << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Chg conj my_pi-:" <<
EvtPDL::chargeConj(myidpim) << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Chg conj my_pi0:" <<
EvtPDL::chargeConj(myidpi0) << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runRepeat(int nevent) {
int i;
for(i=0;i<nevent;i++){
EvtDecayTable::getInstance()->readDecayFile(std::string("../DECAY_2010.DEC"));
}
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runPhotos(int nevent, EvtGen &myGenerator) {
static EvtId PSI=EvtPDL::getId(std::string("J/psi"));
TFile *file=new TFile("photos.root", "RECREATE");
TH1F* mee = new TH1F("h1","mee",60,3.0,3.12);
int count=1;
EvtVector4R e1,e2;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/PHOTOS.DEC");
//EvtGen myGenerator(decay_name,pdttable_name,myRandomEngine);
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMass(PSI),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(PSI,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
e1=root_part->getDaug(0)->getP4Lab();
e2=root_part->getDaug(1)->getP4Lab();
mee->Fill( (e1+e2).mass() );
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runFinalStates(int nevent, EvtGen &myGenerator) {
//Parse the table of particles to find..
EvtParser parser;
parser.read(std::string("exampleFiles/finalstates.list"));
std::vector< std::string > dList[20];
int dListNum[20];
std::vector< std::string > *dListItem = 0;
std::string dListName[20];
int ik,lk;
std::string tk="";
int tline=-1;
std::string parent;
for(ik=0;ik<parser.getNToken();ik++){
lk=tline;
tline=parser.getLineofToken(ik);
tk=parser.getToken(ik);
if ( lk!=tline&&tline>2) {
if ( tline>1) {
dList[tline-3]=*dListItem;
dListItem=new std::vector<std::string>;
dListNum[tline-3]=0;
dListName[tline-2]=parser.getToken(ik);
}
}
else{
if ( tline==1 ) {
//This is the parent particle name
parent=parser.getToken(ik);
dListItem=new std::vector<std::string>;
}
else{
//This is one of the daughters
if ( tline!=2 || (lk==tline) ) {
dListItem->push_back(parser.getToken(ik));
}
if ( tline==2 && (lk!=tline) ) {
dListName[tline-2]=parser.getToken(ik);
}
}
}
}
dList[tline-2]=*dListItem;
dListNum[tline-2]=0;
static EvtId parId=EvtPDL::getId(parent);
int count=0;
do{
if (count==1000*(count/1000)) {
//if (count==1*(count/1)) {
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Event:"<< count << std::endl;
//EvtGenReport(EVTGEN_INFO,"EvtGen") << HepRandom::getTheSeed()<<std::endl;
}
EvtVector4R p_init(EvtPDL::getMass(parId),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(parId,
p_init);
if ( parent=="Upsilon(4S)") {
root_part->setVectorSpinDensity();
}
else{
root_part->setDiagonalSpinDensity();
}
myGenerator.generateDecay(root_part);
EvtParticle* p = root_part;
std::vector<std::string> fs=findFinalState(p);
int j;
for(j=0; j<(tline-1); j++) {
std::vector<std::string> temp=dList[j];
if ( temp.size() == fs.size() ) {
bool foundIt=true;
unsigned int k, l;
std::vector<bool> alreadyUsed(temp.size());
for(k=0; k<temp.size(); k++) {
bool foundThisOne=false;
for(l=0; l<temp.size(); l++) {
if ( k==0 ) alreadyUsed[l]=false;
if ( foundThisOne||alreadyUsed[l] ) continue;
if ( temp[k]==fs[l] ) {
alreadyUsed[l]=true;
foundThisOne=true;
// EvtGenReport(EVTGEN_INFO,"EvtGen") << "found daughter " << k << " " << l << std::endl;
}
}
if ( !foundThisOne ) foundIt=false;
}
if ( foundIt ) {//EvtGenReport(EVTGEN_INFO,"EvtGen") << "found a cand \n"; (histo1[j])->Fill(0.5);
//EvtGenReport(EVTGEN_INFO,"EvtGen") << "found one " << j << std::endl;
dListNum[j]++;}
}
}
root_part->deleteTree();
count++;
}while(count<nevent);
int j;
for (j=0; j<(tline-1) ; j++) EvtGenReport(EVTGEN_INFO,"EvtGen") << dListName[j].c_str() << " " << j << " " << dListNum[j] << " " << count << " "<< (dListNum[j]/(1.0*count)) << std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
std::vector<std::string> findFinalState( EvtParticle *tree){
EvtParticle *p=tree;
std::vector<std::string> fs;
static EvtId ep=EvtPDL::getId(std::string("e+"));
static EvtId em=EvtPDL::getId(std::string("e-"));
static EvtId kp=EvtPDL::getId(std::string("K+"));
static EvtId km=EvtPDL::getId(std::string("K-"));
static EvtId mup=EvtPDL::getId(std::string("mu+"));
static EvtId mum=EvtPDL::getId(std::string("mu-"));
static EvtId pip=EvtPDL::getId(std::string("pi+"));
static EvtId pim=EvtPDL::getId(std::string("pi-"));
static EvtId pi0=EvtPDL::getId(std::string("pi0"));
static EvtId pr=EvtPDL::getId(std::string("p+"));
static EvtId apr=EvtPDL::getId(std::string("anti-p-"));
static EvtId ne=EvtPDL::getId(std::string("n0"));
static EvtId ane=EvtPDL::getId(std::string("anti-n0"));
do{
EvtId type=p->getId();
if (type==ep) fs.push_back(std::string("e+"));
if (type==em) fs.push_back(std::string("e-"));
if (type==mup) fs.push_back(std::string("mu+"));
if (type==mum) fs.push_back(std::string("mu-"));
if (type==kp) fs.push_back(std::string("K+"));
if (type==km) fs.push_back(std::string("K-"));
if (type==pip) fs.push_back(std::string("pi+"));
if (type==pim) fs.push_back(std::string("pi-"));
if (type==pi0) fs.push_back(std::string("pi0"));
if (type==pr) fs.push_back(std::string("p+"));
if (type==apr) fs.push_back(std::string("anti-p-"));
if (type==ne) fs.push_back(std::string("n0"));
if (type==ane) fs.push_back(std::string("anti-n0"));
p=p->nextIter();
}while(p!=0);
return fs;
}
void runTrackMult(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("trackmult.root","RECREATE");
TH1F* trackAll = new TH1F("trackAll","trackAll",12,1.0,25.0);
TH1F* trackNoSL = new TH1F("trackNoSL","trackNoSL",12,1.0,25.0);
TH1F* track1SL = new TH1F("track1SL","track1SL",12,1.0,25.0);
TH1F* track2SL = new TH1F("track2SL","track2SL",12,1.0,25.0);
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
static EvtId BP=EvtPDL::getId(std::string("B+"));
static EvtId BM=EvtPDL::getId(std::string("B-"));
//look for these tracks in generic events
static EvtId ep=EvtPDL::getId(std::string("e+"));
static EvtId em=EvtPDL::getId(std::string("e-"));
static EvtId mup=EvtPDL::getId(std::string("mu+"));
static EvtId mum=EvtPDL::getId(std::string("mu-"));
static EvtId pip=EvtPDL::getId(std::string("pi+"));
static EvtId pim=EvtPDL::getId(std::string("pi-"));
static EvtId kp=EvtPDL::getId(std::string("K+"));
static EvtId km=EvtPDL::getId(std::string("K-"));
static EvtId pp=EvtPDL::getId(std::string("p+"));
static EvtId pm=EvtPDL::getId(std::string("anti-p-"));
static EvtIdSet theTracks(ep,em,mup,mum,pip,pim,kp,km,pp,pm);
static EvtIdSet theLeptons(ep,em,mup,mum);
static EvtIdSet theBs(B0,B0B,BP,BM);
int count=1;
EvtParticle *p;
myGenerator.readUDecay("exampleFiles/GENERIC.DEC");
int totTracks=0;
int totTracksNoSL=0;
int totTracks1SL=0;
int totTracks2SL=0;
int totNoSL=0;
int tot1SL=0;
int tot2SL=0;
do{
int evTracks=0;
if (count==1000*(count/1000)) {
EvtGenReport(EVTGEN_INFO,"EvtGen") << count << std::endl;
//EvtGenReport(EVTGEN_INFO,"EvtGen") << HepRandom::getTheSeed()<<std::endl;
}
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
p=root_part;
int howManySL=0;
do{
if (theTracks.contains(p->getId())) {
totTracks+=1;
evTracks+=1;
}
if ( theLeptons.contains(p->getId())) {
if ( p->getParent() ) {
if ( theBs.contains(p->getParent()->getId())) howManySL+=1;
}
}
p=p->nextIter(root_part);
}while(p!=0);
//Now need to figure out which histogram to book
trackAll->Fill(evTracks);
if ( howManySL==0 ) {trackNoSL->Fill(evTracks); totNoSL+=1; totTracksNoSL+=evTracks;}
if ( howManySL==1 ) {track1SL->Fill(evTracks); tot1SL+=1; totTracks1SL+=evTracks;}
if ( howManySL==2 ) {track2SL->Fill(evTracks); tot2SL+=1; totTracks2SL+=evTracks;}
root_part->deleteTree();
}while (count++<nevent);
double aveMulti=float(totTracks)/float(nevent);
double aveMultiNoSL=float(totTracksNoSL)/float(totNoSL);
double aveMulti1SL=float(totTracks1SL)/float(tot1SL);
double aveMulti2SL=float(totTracks2SL)/float(tot2SL);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Your average multiplicity="<<aveMulti<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Your average multiplicity for no B->semileptonic events="<<aveMultiNoSL<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Your average multiplicity for 1 B->semileptonic events="<<aveMulti1SL<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Your average multiplicity for 2 B->semileptonic events="<<aveMulti2SL<<std::endl;
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runGeneric(int neventOrig, EvtGen &myGenerator,
std::string listfile) {
int nevent=abs(neventOrig);
//Parse the table of particles to find..
TFile *file=new TFile("generic.root","RECREATE");
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
static EvtId BP=EvtPDL::getId(std::string("B+"));
static EvtId BM=EvtPDL::getId(std::string("B-"));
static EvtIdSet theBs(B0B,B0,BP,BM);
static EvtIdSet theB0B(B0B);
static EvtIdSet theB0(B0);
static EvtIdSet theBP(BP);
static EvtIdSet theBM(BM);
static EvtId D0=EvtPDL::getId(std::string("D0"));
static EvtId D0B=EvtPDL::getId(std::string("anti-D0"));
static EvtId DP=EvtPDL::getId(std::string("D+"));
static EvtId DM=EvtPDL::getId(std::string("D-"));
static EvtIdSet theDs(D0B,D0,DP,DM);
static EvtIdSet theD0B(D0B);
static EvtIdSet theD0(D0);
static EvtIdSet theDP(DP);
static EvtIdSet theDM(DM);
int count;
int stable_list[1];
stable_list[0]=-1;
EvtParticle *p;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/GENERIC.DEC");
myGenerator.readUDecay(udecay_name);
if (listfile!="") {
EvtParser parser;
if (parser.read(listfile)!=0){
EvtGenReport(EVTGEN_ERROR,"EvtGen") << "Will terminate."<<std::endl;
exit(-1);
}
std::vector<TH1F*> histo1(parser.getNToken());
std::vector<TH1F*> histo2(parser.getNToken());
std::vector<TH1F*> massHisto(parser.getNToken());
int ik;
std::string tk,tkname;
for(ik=0;ik<(parser.getNToken()/2);ik++){
tk=parser.getToken(2*ik);
tkname=parser.getToken(1+2*ik);
histo1[ik] = new TH1F(tkname.c_str(),tkname.c_str(), 30,0.0,3.0);
char *directName;
directName=new char[(strlen(tkname.c_str())+8)];
directName = strcpy(directName,tkname.c_str());
directName = strcat(directName,"Direct");
histo2[ik] = new TH1F(directName,directName,
30,0.0,3.0);
delete directName;
char *massName;
massName=new char[(strlen(tkname.c_str())+4)];
massName = strcpy(massName,tkname.c_str());
massName = strcat(massName,"Mass");
massHisto[ik] = new TH1F(massName,massName, 3000,0.0,5.0);
delete massName;
}
count=1;
std::vector<int> temp(parser.getNToken()/2,0);
std::vector<int> tempB(parser.getNToken()/2,0);
std::vector<int> tempB0B(parser.getNToken()/2,0);
std::vector<int> tempB0(parser.getNToken()/2,0);
std::vector<int> tempBP(parser.getNToken()/2,0);
std::vector<int> tempBM(parser.getNToken()/2,0);
std::vector<int> tempD(parser.getNToken()/2,0);
std::vector<int> tempD0B(parser.getNToken()/2,0);
std::vector<int> tempD0(parser.getNToken()/2,0);
std::vector<int> tempDP(parser.getNToken()/2,0);
std::vector<int> tempDM(parser.getNToken()/2,0);
do{
//EvtGenReport(EVTGEN_INFO,"EvtGen") << "new event\n";
if (count==1000*(count/1000)) {
EvtGenReport(EVTGEN_INFO,"EvtGen") << count << std::endl;
//EvtGenReport(EVTGEN_INFO,"EvtGen") << HepRandom::getTheSeed()<<std::endl;
}
EvtVector4R p_init(sqrt(EvtPDL::getMass(UPS4)*EvtPDL::getMass(UPS4)+5.9*5.9),0.0,0.0,5.9);
EvtParticle* root_part=0;
if ( neventOrig > 0 ) {
root_part=EvtParticleFactory::particleFactory(UPS4,p_init);
}
else{
root_part=EvtParticleFactory::particleFactory(VPHO, p_init);
}
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
p=root_part;
//EvtStdHep stdhep;
//stdhep.init();
//root_part->makeStdHep(stdhep);
//EvtGenReport(EVTGEN_INFO,"EvtGen") <<stdhep<<std::endl;
//EvtGenReport(EVTGEN_INFO,"EvtGen") <<secondary<<std::endl;
std::string token;
int itok;
for(itok=0;itok<(parser.getNToken()/2);itok++){
token=parser.getToken(2*itok);
//temp[itok]+=countInclusive(token,root_part);
temp[itok]+=countInclusive(token,root_part,histo1[itok],massHisto[itok]);
tempB[itok]+=countInclusiveParent(token,root_part,theBs,histo2[itok]);
tempB0[itok]+=countInclusiveSubTree(token,root_part,theB0);
tempB0B[itok]+=countInclusiveSubTree(token,root_part,theB0B);
tempBP[itok]+=countInclusiveSubTree(token,root_part,theBP);
tempBM[itok]+=countInclusiveSubTree(token,root_part,theBM);
// tempD[itok]+=countInclusiveParent(token,root_part,theDs);
// tempD0[itok]+=countInclusiveSubTree(token,root_part,theD0);
// tempD0B[itok]+=countInclusiveSubTree(token,root_part,theD0B);
// tempDP[itok]+=countInclusiveSubTree(token,root_part,theDP);
// tempDM[itok]+=countInclusiveSubTree(token,root_part,theDM);
}
// numd0+=countInclusive("D0",root_part);
// numd0b+=countInclusive("anti-D0",root_part);
// numdp+=countInclusive("D+",root_part);
// numdm+=countInclusive("D-",root_part);
// root_part->printTree();
root_part->deleteTree();
}while (count++<nevent);
int itok;
std::string token;
for(itok=0;itok<(parser.getNToken()/2);itok++){
token=parser.getToken(2*itok);
float br=0.5*float(temp[itok])/float(nevent);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Found "<<temp[itok]<<" "<<token.c_str()
<< " in " << nevent << " events. Average number of "
<< token.c_str()<<" per B meson="<<br<<std::endl;
br=0.5*float(tempB[itok])/float(nevent);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Found "<<tempB[itok]<<" "<<token.c_str()
<< " produced directly in decays of B mesons avg. br.fr.="
<<br<<std::endl;
br=2.0*float(tempB0[itok])/float(nevent);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Found "<<tempB0[itok]<<" "<<token.c_str()
<< " in decay tree of B0, br.fr.="<<br<<std::endl;
br=2.0*float(tempB0B[itok])/float(nevent);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Found "<<tempB0B[itok]<<" "<<token.c_str()
<< " in decay tree of anti-B0, br.fr.="<<br<<std::endl;
br=2.0*float(tempBP[itok])/float(nevent);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Found "<<tempBP[itok]<<" "<<token.c_str()
<< " in decay tree of B+, br.fr.="<<br<<std::endl;
br=2.0*float(tempBM[itok])/float(nevent);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Found "<<tempBM[itok]<<" "<<token.c_str()
<< " in decay tree of B-, br.fr.="<<br<<std::endl;
// br=0.5*float(tempD[itok])/float(numd0+numd0b+numdm+numdp);
// EvtGenReport(EVTGEN_INFO,"EvtGen") << "Found "<<tempD[itok]<<" "<<token
// << " produced directly in decays of D mesons avg. br.fr.="
// <<br<<std::endl;
// br=2.0*float(tempD0[itok])/float(numd0);
// EvtGenReport(EVTGEN_INFO,"EvtGen") << "Found "<<tempD0[itok]<<" "<<token
// << " in decay of D0, br.fr.="<<br<<std::endl;
// br=2.0*float(tempD0B[itok])/float(numd0b);
// EvtGenReport(EVTGEN_INFO,"EvtGen") << "Found "<<tempD0B[itok]<<" "<<token
// << " in decay of anti-D0, br.fr.="<<br<<std::endl;
// br=2.0*float(tempDP[itok])/float(numdp);
// EvtGenReport(EVTGEN_INFO,"EvtGen") << "Found "<<tempDP[itok]<<" "<<token
// << " in decay of D+, br.fr.="<<br<<std::endl;
// br=2.0*float(tempDM[itok])/float(numdm);
// EvtGenReport(EVTGEN_INFO,"EvtGen") << "Found "<<tempDM[itok]<<" "<<token
// << " in decay of D-, br.fr.="<<br<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "*******************************************\n";
}
}
else{
count=1;
do{
if (count==1000*(count/1000)) {
EvtGenReport(EVTGEN_INFO,"EvtGen") << count << std::endl;
//EvtGenReport(EVTGEN_INFO,"EvtGen") << HepRandom::getTheSeed()<<std::endl;
}
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=0;
if ( neventOrig > 0 ) {
root_part=EvtParticleFactory::particleFactory(UPS4,p_init);
}
else{
root_part=EvtParticleFactory::particleFactory(VPHO, p_init);
}
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
p=root_part;
root_part->deleteTree();
}while (count++<nevent);
}
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runKstarnunu(int nevent, EvtGen &myGenerator) {
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
TFile *file=new TFile("kstarnunu.root", "RECREATE");
TH1F* q2 = new TH1F("h1","q2",
50,0.0,25.0);
TH1F* enu = new TH1F("h2","Neutrino energy",
50,0.0,5.0);
TH1F* x = new TH1F("h3","Total neutrino energy/B mass",
50,0.5,0.9);
int count;
EvtVector4R kstar,nu,nub;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/KSTARNUNU.DEC");
//EvtGen myGenerator(decay_name,pdttable_name,myRandomEngine);
myGenerator.readUDecay(udecay_name);
count=1;
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(B0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
kstar=root_part->getDaug(0)->getP4Lab();
nu=root_part->getDaug(1)->getP4Lab();
nub=root_part->getDaug(2)->getP4Lab();
q2->Fill( (nu+nub).mass2() );
enu->Fill( nu.get(0) );
x->Fill( (nu.get(0)+nub.get(0))/root_part->mass() );
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runBsmix(int nevent, EvtGen &myGenerator) {
static EvtId BS0=EvtPDL::getId(std::string("B_s0"));
static EvtId BSB=EvtPDL::getId(std::string("anti-B_s0"));
TFile *file=new TFile("bsmix.root", "RECREATE");
TH1F* tmix = new TH1F("h1","tmix (mm)",100,0.0,5.0);
TH1F* tnomix = new TH1F("h2","tnomix (mm)",100,0.0,5.0);
int count;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/BSMIX.DEC");
//EvtGen myGenerator(decay_name,pdttable_name,myRandomEngine);
myGenerator.readUDecay(udecay_name);
count=1;
do{
EvtVector4R p_init(EvtPDL::getMass(BS0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(BS0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
double t=root_part->getLifetime();
int mix=0;
if (root_part->getNDaug()==1){
if (root_part->getDaug(0)->getId()==BSB){
mix=1;
}
}
if (mix==0) tnomix->Fill(t);
if (mix==1) tmix->Fill(t);
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runSemic(int nevent, EvtGen &myGenerator) {
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
static EvtId GAMM=EvtPDL::getId(std::string("gamma"));
static EvtId PSI=EvtPDL::getId(std::string("J/psi"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
static EvtId BS0=EvtPDL::getId(std::string("B_s0"));
static EvtId BSB=EvtPDL::getId(std::string("anti-B_s0"));
static EvtId PI0=EvtPDL::getId(std::string("pi0"));
static EvtId PIP=EvtPDL::getId(std::string("pi+"));
static EvtId PIM=EvtPDL::getId(std::string("pi-"));
static EvtId EP=EvtPDL::getId(std::string("e+"));
static EvtId KP=EvtPDL::getId(std::string("K+"));
static EvtId MUP=EvtPDL::getId(std::string("mu+"));
static EvtId EM=EvtPDL::getId(std::string("e-"));
static EvtId KM=EvtPDL::getId(std::string("K-"));
static EvtId MUM=EvtPDL::getId(std::string("mu-"));
static EvtId DST0=EvtPDL::getId(std::string("D*0"));
static EvtId DSTB=EvtPDL::getId(std::string("anti-D*0"));
static EvtId DSTP=EvtPDL::getId(std::string("D*+"));
static EvtId DSTM=EvtPDL::getId(std::string("D*-"));
static EvtId D0=EvtPDL::getId(std::string("D0"));
static EvtId D0B=EvtPDL::getId(std::string("anti-D0"));
static EvtId DP=EvtPDL::getId(std::string("D+"));
static EvtId DM=EvtPDL::getId(std::string("D-"));
static EvtId K0L=EvtPDL::getId(std::string("K_L0"));
static EvtId D1P1P=EvtPDL::getId(std::string("D_1+"));
static EvtId D1P1N=EvtPDL::getId(std::string("D_1-"));
static EvtId D1P10=EvtPDL::getId(std::string("D_10"));
static EvtId D1P1B=EvtPDL::getId(std::string("anti-D_10"));
static EvtId D3P2P=EvtPDL::getId(std::string("D_2*+"));
static EvtId D3P2N=EvtPDL::getId(std::string("D_2*-"));
static EvtId D3P20=EvtPDL::getId(std::string("D_2*0"));
static EvtId D3P2B=EvtPDL::getId(std::string("anti-D_2*0"));
static EvtId D3P1P=EvtPDL::getId(std::string("D'_1+"));
static EvtId D3P1N=EvtPDL::getId(std::string("D'_1-"));
static EvtId D3P10=EvtPDL::getId(std::string("D'_10"));
static EvtId D3P1B=EvtPDL::getId(std::string("anti-D'_10"));
static EvtId D3P0P=EvtPDL::getId(std::string("D_0*+"));
static EvtId D3P0N=EvtPDL::getId(std::string("D_0*-"));
static EvtId D3P00=EvtPDL::getId(std::string("D_0*0"));
static EvtId D3P0B=EvtPDL::getId(std::string("anti-D_0*0"));
static EvtId D21S0P=EvtPDL::getId(std::string("D(2S)+"));
static EvtId D21S0N=EvtPDL::getId(std::string("D(2S)-"));
static EvtId D21S00=EvtPDL::getId(std::string("D(2S)0"));
static EvtId D21S0B=EvtPDL::getId(std::string("anti-D(2S)0"));
static EvtId D23S1P=EvtPDL::getId(std::string("D*(2S)+"));
static EvtId D23S1N=EvtPDL::getId(std::string("D*(2S)-"));
static EvtId D23S10=EvtPDL::getId(std::string("D*(2S)0"));
static EvtId D23S1B=EvtPDL::getId(std::string("anti-D*(2S)0"));
static EvtIdSet radExitDstar(D23S1P,D23S1N,D23S10,D23S1B);
TFile *file=new TFile("semic.root", "RECREATE");
TH1F* Dpe_q2 = new TH1F("h11","q2 for B0B ->D+ e- nu",50,0.0,12.0);
TH1F* Dpe_elep = new TH1F("h12","Elep for B0B ->D+ e- nu",50,0.0,2.5);
TH1F* Dme_q2 = new TH1F("h13","q2 for B0 ->D- e+ nu",50,0.0,12.0);
TH1F* Dme_elep = new TH1F("h14","Elep for B0 ->D- e+ nu",50,0.0,2.5);
TH1F* D0e_q2 = new TH1F("h15","q2 for B- ->D0 e- nu",50,0.0,12.0);
TH1F* D0e_elep = new TH1F("h16","Elep for B- ->D0 e- nu",50,0.0,2.5);
TH1F* D0Be_q2 = new TH1F("h17","q2 for B+ ->D0B e+ nu",50,0.0,12.0);
TH1F* D0Be_elep = new TH1F("h18","Elep for B+ ->D0B e+ nu",50,0.0,2.5);
TH1F* Dstpe_q2 = new TH1F("h21","q2 for B0B ->D*+ e- nu",50,0.0,12.0);
TH1F* Dstpe_elep = new TH1F("h22","Elep for B0B ->D*+ e- nu",50,0.0,2.5);
TH1F* Dstme_q2 = new TH1F("h23","q2 for B0 ->D*- e+ nu",50,0.0,12.0);
TH1F* Dstme_elep = new TH1F("h24","Elep for B0 ->D*- e+ nu",50,0.0,2.5);
TH1F* Dst0e_q2 = new TH1F("h25","q2 for B- ->D*0 e- nu",50,0.0,12.0);
TH1F* Dst0e_elep = new TH1F("h26","Elep for B*- ->D*0 e- nu",50,0.0,2.5);
TH1F* Dst0Be_q2 = new TH1F("h27","q2 for B+ ->D*0B e+ nu",50,0.0,12.0);
TH1F* Dst0Be_elep = new TH1F("h28","Elep for B+ ->D*0B e+ nu",50,0.0,2.5);
TH1F* D1P1pe_q2 = new TH1F("h31","q2 for B0B ->1P1+ e- nu",50,0.0,12.0);
TH1F* D1P1pe_elep = new TH1F("h32","Elep for B0B ->1P1+ e- nu",50,0.0,2.5);
TH1F* D1P1me_q2 = new TH1F("h33","q2 for B0 ->1P1- e+ nu",50,0.0,12.0);
TH1F* D1P1me_elep = new TH1F("h34","Elep for B0 ->1P1- e+ nu",50,0.0,2.5);
TH1F* D1P10e_q2 = new TH1F("h35","q2 for B- ->1P10 e- nu",50,0.0,12.0);
TH1F* D1P10e_elep = new TH1F("h36","Elep for B*- ->1P10 e- nu",50,0.0,2.5);
TH1F* D1P10Be_q2 = new TH1F("h37","q2 for B+ ->1P1B e+ nu",50,0.0,12.0);
TH1F* D1P10Be_elep = new TH1F("h38","Elep for B+ ->1P1B e+ nu",50,0.0,2.5);
TH1F* D3P0pe_q2 = new TH1F("h41","q2 for B0B ->3P0+ e- nu",50,0.0,12.0);
TH1F* D3P0pe_elep = new TH1F("h42","Elep for B0B ->3P0+ e- nu",50,0.0,2.5);
TH1F* D3P0me_q2 = new TH1F("h43","q2 for B0 ->3P0- e+ nu",50,0.0,12.0);
TH1F* D3P0me_elep = new TH1F("h44","Elep for B0 ->3P0- e+ nu",50,0.0,2.5);
TH1F* D3P00e_q2 = new TH1F("h45","q2 for B- ->3P00 e- nu",50,0.0,12.0);
TH1F* D3P00e_elep = new TH1F("h46","Elep for B*- ->3P00 e- nu",50,0.0,2.5);
TH1F* D3P00Be_q2 = new TH1F("h47","q2 for B+ ->3P0B e+ nu",50,0.0,12.0);
TH1F* D3P00Be_elep = new TH1F("h48","Elep for B+ ->3P0B e+ nu",50,0.0,2.5);
TH1F* D3P1pe_q2 = new TH1F("h51","q2 for B0B ->3P1+ e- nu",50,0.0,12.0);
TH1F* D3P1pe_elep = new TH1F("h52","Elep for B0B ->3P1+ e- nu",50,0.0,2.5);
TH1F* D3P1me_q2 = new TH1F("h53","q2 for B0 ->3P1- e+ nu",50,0.0,12.0);
TH1F* D3P1me_elep = new TH1F("h54","Elep for B0 ->3P1- e+ nu",50,0.0,2.5);
TH1F* D3P10e_q2 = new TH1F("h55","q2 for B- ->3P10 e- nu",50,0.0,12.0);
TH1F* D3P10e_elep = new TH1F("h56","Elep for B*- ->3P10 e- nu",50,0.0,2.5);
TH1F* D3P10Be_q2 = new TH1F("h57","q2 for B+ ->3P1B e+ nu",50,0.0,12.0);
TH1F* D3P10Be_elep = new TH1F("h58","Elep for B+ ->3P1B e+ nu",50,0.0,2.5);
TH1F* D3P2pe_q2 = new TH1F("h61","q2 for B0B ->3P2+ e- nu",50,0.0,12.0);
TH1F* D3P2pe_elep = new TH1F("h62","Elep for B0B ->3P2+ e- nu",50,0.0,2.5);
TH1F* D3P2me_q2 = new TH1F("h63","q2 for B0 ->3P2- e+ nu",50,0.0,12.0);
TH1F* D3P2me_elep = new TH1F("h64","Elep for B0 ->3P2- e+ nu",50,0.0,2.5);
TH1F* D3P20e_q2 = new TH1F("h65","q2 for B- ->3P20 e- nu",50,0.0,12.0);
TH1F* D3P20e_elep = new TH1F("h66","Elep for B*- ->3P20 e- nu",50,0.0,2.5);
TH1F* D3P20Be_q2 = new TH1F("h67","q2 for B+ ->3P2B e+ nu",50,0.0,12.0);
TH1F* D3P20Be_elep = new TH1F("h68","Elep for B+ ->3P2B e+ nu",50,0.0,2.5);
TH1F* phiL = new TH1F("h69","phi",50,-3.1416,3.1416);
int count;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/SEMIC.DEC");
//EvtGen myGenerator(decay_name,pdttable_name,myRandomEngine);
myGenerator.readUDecay(udecay_name);
count=1;
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
int i;
for (i=0;i<2;i++){
EvtId meson=root_part->getDaug(i)->getDaug(0)->getId();
EvtId lepton=root_part->getDaug(i)->getDaug(1)->getId();
EvtVector4R lep=root_part->getDaug(i)->getDaug(1)->getP4Lab();
phiL->Fill(atan2(lep.get(1),lep.get(2)));
EvtVector4R nu=root_part->getDaug(i)->getDaug(2)->getP4Lab();
double q2=(lep+nu).mass2();
double elep=root_part->getDaug(i)->getDaug(1)->getP4().get(0);
if (meson==DP&&lepton==EM){
Dpe_q2->Fill(q2);
Dpe_elep->Fill(elep);
}
if (meson==DM&&lepton==EP){
Dme_q2->Fill(q2);
Dme_elep->Fill(elep);
}
if (meson==D0&&lepton==EM){
D0e_q2->Fill(q2);
D0e_elep->Fill(elep);
}
if (meson==D0B&&lepton==EP){
D0Be_q2->Fill(q2);
D0Be_elep->Fill(elep);
}
if (meson==DSTP&&lepton==EM){
Dstpe_q2->Fill(q2);
Dstpe_elep->Fill(elep);
}
if (meson==DSTM&&lepton==EP){
Dstme_q2->Fill(q2);
Dstme_elep->Fill(elep);
}
if (meson==DST0&&lepton==EM){
Dst0e_q2->Fill(q2);
Dst0e_elep->Fill(elep);
}
if (meson==DSTB&&lepton==EP){
Dst0Be_q2->Fill(q2);
Dst0Be_elep->Fill(elep);
}
if (meson==D1P1P&&lepton==EM){
D1P1pe_q2->Fill(q2);
D1P1pe_elep->Fill(elep);
}
if (meson==D1P1N&&lepton==EP){
D1P1me_q2->Fill(q2);
D1P1me_elep->Fill(elep);
}
if (meson==D1P10&&lepton==EM){
D1P10e_q2->Fill(q2);
D1P10e_elep->Fill(elep);
}
if (meson==D1P1B&&lepton==EP){
D1P10Be_q2->Fill(q2);
D1P10Be_elep->Fill(elep);
}
if (meson==D3P0P&&lepton==EM){
D3P0pe_q2->Fill(q2);
D3P0pe_elep->Fill(elep);
}
if (meson==D3P0N&&lepton==EP){
D3P0me_q2->Fill(q2);
D3P0me_elep->Fill(elep);
}
if (meson==D3P00&&lepton==EM){
D3P00e_q2->Fill(q2);
D3P00e_elep->Fill(elep);
}
if (meson==D3P0B&&lepton==EP){
D3P00Be_q2->Fill(q2);
D3P00Be_elep->Fill(elep);
}
if (meson==D3P1P&&lepton==EM){
D3P1pe_q2->Fill(q2);
D3P1pe_elep->Fill(elep);
}
if (meson==D3P1N&&lepton==EP){
D3P1me_q2->Fill(q2);
D3P1me_elep->Fill(elep);
}
if (meson==D3P10&&lepton==EM){
D3P10e_q2->Fill(q2);
D3P10e_elep->Fill(elep);
}
if (meson==D3P1B&&lepton==EP){
D3P10Be_q2->Fill(q2);
D3P10Be_elep->Fill(elep);
}
if (meson==D3P2P&&lepton==EM){
D3P2pe_q2->Fill(q2);
D3P2pe_elep->Fill(elep);
}
if (meson==D3P2N&&lepton==EP){
D3P2me_q2->Fill(q2);
D3P2me_elep->Fill(elep);
}
if (meson==D3P20&&lepton==EM){
D3P20e_q2->Fill(q2);
D3P20e_elep->Fill(elep);
}
if (meson==D3P2B&&lepton==EP){
D3P20Be_q2->Fill(q2);
D3P20Be_elep->Fill(elep);
}
}
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runKstarll(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("kstkmm.root", "RECREATE");
TH2F* _dalitz = new TH2F("h1","q^2! vs Elep",
70,0.0,3.5,60,0.0,30.0);
TH1F* _ctl = new TH1F("h2","ctl",50,-1.0,1.0);
TH1F* _q2 = new TH1F("h3","q2",50,0.0,25.0);
TH1F* _q2low = new TH1F("h4","q2 (low)",
50,0.0,1.0);
TH1F* _q2lowlow = new TH1F("h5","q2 (lowlow)",
50,0.0,0.00001);
TH1F* _phi = new TH1F("h6","phi",50,-EvtConst::pi,EvtConst::pi);
TH1F* _chi = new TH1F("h7","chi",50,0.0,EvtConst::twoPi);
TH1F* _chictl = new TH1F("h8","chictl",50,0.0,EvtConst::twoPi);
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
int count=1;
EvtVector4R kstar,l1,l2;
EvtVector4R k,pi,b;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/KSTARLL.DEC");
//EvtGen myGenerator(decay_name,pdttable_name,myRandomEngine);
myGenerator.readUDecay(udecay_name);
std::vector<double> q2low(5);
std::vector<double> q2high(5);
std::vector<int> counts(5);
//kee
//int n=4;
//q2low[0]=0.0; q2high[0]=4.5;
//q2low[1]=4.5; q2high[1]=8.41;
//q2low[2]=10.24; q2high[2]=12.96;
//q2low[3]=14.06; q2high[3]=30.0;
//kmm
//int n=4;
//q2low[0]=0.0; q2high[0]=4.5;
//q2low[1]=4.5; q2high[1]=9.0;
//q2low[2]=10.24; q2high[2]=12.96;
//q2low[3]=14.06; q2high[3]=30.0;
//K*ee
int n=5;
q2low[0]=0.0; q2high[0]=0.1;
q2low[1]=0.1; q2high[1]=4.5;
q2low[2]=4.5; q2high[2]=8.41;
q2low[3]=10.24; q2high[3]=12.96;
q2low[4]=14.06; q2high[4]=30.0;
//K*mm
//int n=5;
//q2low[0]=0.0; q2high[0]=0.1;
//q2low[1]=0.1; q2high[1]=4.5;
//q2low[2]=4.5; q2high[2]=9.0;
//q2low[3]=10.24; q2high[3]=12.96;
//q2low[4]=14.06; q2high[4]=30.0;
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(B0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
// root_part->printTree();
//root_part->getDaug(0)->printTree();
//root_part->getDaug(1)->printTree();
//root_part->getDaug(2)->printTree();
kstar=root_part->getDaug(0)->getP4Lab();
l1=root_part->getDaug(1)->getP4Lab();
l2=root_part->getDaug(2)->getP4Lab();
b=root_part->getP4();
k=root_part->getDaug(0)->getDaug(0)->getP4Lab();
pi=root_part->getDaug(0)->getDaug(1)->getP4Lab();
double qsq=(l1+l2).mass2();
for (int j=0;j<n;j++){
if (qsq>q2low[j]&&qsq<q2high[j]) counts[j]++;
}
_q2->Fill( (l1+l2).mass2() );
_q2low->Fill( (l1+l2).mass2() );
_q2lowlow->Fill( (l1+l2).mass2() );
_ctl->Fill(EvtDecayAngle((l1+l2+kstar),(l1+l2),l1));
_dalitz->Fill(l1.get(0),(l1+l2).mass2(),1.0);
root_part->deleteTree();
_phi->Fill(atan2(l1.get(1),l1.get(2)));
_chi->Fill(EvtDecayAngleChi(b,k,pi,l1,l2));
if (EvtDecayAngle((l1+l2+kstar),(l1+l2),l1)>0){
_chictl->Fill(EvtDecayAngleChi(b,k,pi,l1,l2));
}
EvtGenReport(EVTGEN_INFO,"EvtGen") << "count:"<<count<<" "<<(l1+l2).mass2()<<std::endl;
}while (count++<nevent);
for (int j=0;j<n;j++){
std::cout << "["<<q2low[j]<<".."<<q2high[j]<<"]="<<counts[j]<<std::endl;
}
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runKll(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("ksem.root", "RECREATE");
TH2F* _dalitz = new TH2F("h1","q^2! vs Elep",
70,0.0,3.5,60,0.0,30.0);
TH1F* _ctl = new TH1F("h2","ctl",50,-1.0,1.0);
TH1F* _q2 = new TH1F("h3","q2",50,0.0,25.0);
TH1F* _q2low = new TH1F("h4","q2 (low)",
50,0.0,1.0);
TH1F* _q2lowlow = new TH1F("h5","q2 (lowlow)",
50,0.0,0.00001);
TH1F* _phi = new TH1F("h6","phi",50,-EvtConst::pi,EvtConst::pi);
// TH1F* _chi = new TH1F("h7","chi",50,0.0,EvtConst::twoPi);
// TH1F* _chictl = new TH1F("h8","chictl",50,0.0,EvtConst::twoPi);
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
int count=1;
EvtVector4R k,l1,l2;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/KLL.DEC");
//EvtGen myGenerator(decay_name,pdttable_name,myRandomEngine);
myGenerator.readUDecay(udecay_name);
std::vector<double> q2low(5);
std::vector<double> q2high(5);
std::vector<int> counts(5);
//kee
// int n=4;
//q2low[0]=0.0; q2high[0]=4.5;
//q2low[1]=4.5; q2high[1]=8.41;
//q2low[2]=10.24; q2high[2]=12.96;
//q2low[3]=14.06; q2high[3]=30.0;
//kmm
int n=4;
q2low[0]=0.0; q2high[0]=4.5;
q2low[1]=4.5; q2high[1]=9.0;
q2low[2]=10.24; q2high[2]=12.96;
q2low[3]=14.06; q2high[3]=30.0;
//K*ee
//int n=5;
//q2low[0]=0.0; q2high[0]=0.1;
//q2low[1]=0.1; q2high[1]=4.5;
//q2low[2]=4.5; q2high[2]=8.41;
//q2low[3]=10.24; q2high[3]=12.96;
//q2low[4]=14.06; q2high[4]=30.0;
//K*mm
//int n=5;
//q2low[0]=0.0; q2high[0]=0.1;
//q2low[1]=0.1; q2high[1]=4.5;
//q2low[2]=4.5; q2high[2]=9.0;
//q2low[3]=10.24; q2high[3]=12.96;
//q2low[4]=14.06; q2high[4]=30.0;
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(B0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
// root_part->printTree();
//root_part->getDaug(0)->printTree();
//root_part->getDaug(1)->printTree();
//root_part->getDaug(2)->printTree();
k=root_part->getDaug(0)->getP4Lab();
l1=root_part->getDaug(1)->getP4Lab();
l2=root_part->getDaug(2)->getP4Lab();
//b=root_part->getP4();
// k=root_part->getDaug(0)->getDaug(0)->getP4Lab();
// pi=root_part->getDaug(0)->getDaug(1)->getP4Lab();
double qsq=(l1+l2).mass2();
for (int j=0;j<n;j++){
if (qsq>q2low[j]&&qsq<q2high[j]) counts[j]++;
}
_q2->Fill( (l1+l2).mass2() );
_q2low->Fill( (l1+l2).mass2() );
_q2lowlow->Fill( (l1+l2).mass2() );
_ctl->Fill(EvtDecayAngle((l1+l2+k),(l1+l2),l1));
_dalitz->Fill(l1.get(0),(l1+l2).mass2(),1.0);
root_part->deleteTree();
_phi->Fill(atan2(l1.get(1),l1.get(2)));
//_chi->Fill(EvtDecayAngleChi(b,k,pi,l1,l2));
//if (EvtDecayAngle((l1+l2+kstar),(l1+l2),l1)>0){
// _chictl->Fill(EvtDecayAngleChi(b,k,pi,l1,l2));
// }
EvtGenReport(EVTGEN_INFO,"EvtGen") << "count:"<<count<<" "<<(l1+l2).mass2()<<std::endl;
}while (count++<nevent);
for (int j=0;j<n;j++){
std::cout << "["<<q2low[j]<<".."<<q2high[j]<<"]="<<counts[j]<<std::endl;
}
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runHll(int nevent, EvtGen &myGenerator, char* mode) {
TString modename = mode;
TString filename;
filename = modename;
filename = filename + "_nnlo.root";
TFile *file=new TFile(filename, "RECREATE");
TString decname;
decname += modename;
decname.ToUpper();
decname = "exampleFiles/" + decname + ".DEC";
char udecay_name[100];
strcpy(udecay_name,decname);
TH2F* _dalitz = new TH2F("h1","q^2! vs Elep",
70,0.0,3.5,60,0.0,30.0);
TH1F* _ctl = new TH1F("h2","ctl",50,-1.0,1.0);
TH1F* _q2 = new TH1F("h3","q2",50,0.0,25.0);
TH1F* _q2low = new TH1F("h4","q2 (low)",
50,0.0,1.0);
TH1F* _q2lowlow = new TH1F("h5","q2 (lowlow)",
50,0.0,0.00001);
TH1F* _phi = new TH1F("h6","phi",50,-EvtConst::pi,EvtConst::pi);
TH1F* _chi = new TH1F("h7","chi",50,0.0,EvtConst::twoPi);
TH1F* _chictl = new TH1F("h8","chictl",50,0.0,EvtConst::twoPi);
EvtId B;
if (modename == "kee" || modename == "kmm" ||
modename == "kstksee" || modename == "kstksmm" ||
modename == "piee" || modename == "pimm" ||
modename == "rhoee" || modename == "rhomm")
{
B=EvtPDL::getId(std::string("B+"));
}
else
{
B=EvtPDL::getId(std::string("B0"));
}
int count = 1;
EvtVector4R b,h,l1,l2;
EvtVector4R hdaug1,hdaug2;
myGenerator.readUDecay(udecay_name);
std::vector<double> q2low(7);
std::vector<double> q2high(7);
std::vector<int> counts(7);
int n;
if (modename == "kee" || modename == "ksee" ||
modename == "piee" || modename == "pi0ee" ||
modename == "etaee" || modename == "etapee")
{
//kee
n = 6;
q2low[0] = 0.0; q2high[0] = 4.5;
q2low[1] = 4.5; q2high[1] = 8.41;
q2low[2] = 8.41; q2high[2] = 10.24;
q2low[3] = 10.24; q2high[3] = 12.96;
q2low[4] = 12.96; q2high[4] = 14.06;
q2low[5] = 14.06; q2high[5] = 30.0;
}
else if (modename == "kmm" || modename == "ksmm" ||
modename == "pimm" || modename == "pi0mm" ||
modename == "etamm" || modename == "etapmm")
{
//kmm
n = 6;
q2low[0] = 0.0; q2high[0] = 4.5;
q2low[1] = 4.5; q2high[1] = 9.0;
q2low[2] = 9.0; q2high[2] = 10.24;
q2low[3] = 10.24; q2high[3] = 12.96;
q2low[4] = 12.96; q2high[4] = 14.06;
q2low[5] = 14.06; q2high[5] = 30.0;
}
else if (modename == "kstkee" || modename == "kstksee" ||
modename == "rhoee" || modename == "rho0ee" || modename == "omegaee")
{
//K*ee
n = 7;
q2low[0] = 0.0; q2high[0] = 0.1;
q2low[1] = 0.1; q2high[1] = 4.5;
q2low[2] = 4.5; q2high[2] = 8.41;
q2low[3] = 8.41; q2high[3] = 10.24;
q2low[4] = 10.24; q2high[4] = 12.96;
q2low[5] = 12.96; q2high[5] = 14.06;
q2low[6] = 14.06; q2high[6] = 30.0;
}
else if (modename == "kstkmm" || modename == "kstksmm" ||
modename == "rhomm" || modename == "rho0mm" || modename == "omegamm")
{
//K*mm
n = 7;
q2low[0] = 0.0; q2high[0] = 0.1;
q2low[1] = 0.1; q2high[1] = 4.5;
q2low[2] = 4.5; q2high[2] = 9.0;
q2low[3] = 9.0; q2high[3] = 10.24;
q2low[4] = 10.24; q2high[4] = 12.96;
q2low[5] = 12.96; q2high[5] = 14.06;
q2low[6] = 14.06; q2high[6] = 30.0;
}
float q2binlow[n+1];
for (int i = 0; i < n; i++)
{
q2binlow[i] = q2low[i];
}
q2binlow[n] = 30.0;
TH1F* _q2var = new TH1F("h9", "q2var", n, q2binlow);
do
{
EvtVector4R p_init(EvtPDL::getMass(B),0.0,0.0,0.0);
EvtParticle* root_part = EvtParticleFactory::particleFactory(B,p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
// root_part->printTree();
//root_part->getDaug(0)->printTree();
//root_part->getDaug(1)->printTree();
//root_part->getDaug(2)->printTree();
h = root_part->getDaug(0)->getP4Lab();
l1 = root_part->getDaug(1)->getP4Lab();
l2 = root_part->getDaug(2)->getP4Lab();
double qsq=(l1+l2).mass2();
for (int j = 0 ; j < n; j++)
{
if (qsq > q2low[j] && qsq < q2high[j]) counts[j]++;
}
_q2->Fill( (l1+l2).mass2() );
_q2var->Fill( (l1+l2).mass2() );
_q2low->Fill( (l1+l2).mass2() );
_q2lowlow->Fill( (l1+l2).mass2() );
_ctl->Fill(EvtDecayAngle((l1+l2+h),(l1+l2),l1));
_dalitz->Fill(l1.get(0),(l1+l2).mass2(),1.0);
_phi->Fill(atan2(l1.get(1),l1.get(2)));
if (modename == "kstkee" || modename == "kstkmm" ||
modename == "kstksee" || modename == "kstksmm" ||
modename == "rhoee" || modename == "rhomm" ||
modename == "rho0ee" || modename == "rho0mm")
{
b = root_part->getP4();
hdaug1 = root_part->getDaug(0)->getDaug(0)->getP4Lab();
hdaug2 = root_part->getDaug(0)->getDaug(1)->getP4Lab();
_chi->Fill(EvtDecayAngleChi(b,hdaug1,hdaug2,l1,l2));
if (EvtDecayAngle((l1+l2+h),(l1+l2),l1) > 0)
{
_chictl->Fill(EvtDecayAngleChi(b,hdaug1,hdaug2,l1,l2));
}
}
if (count % 1000 == 0)
{
EvtGenReport(EVTGEN_INFO,"EvtGen") << "count:"<<count<<" "<<(l1+l2).mass2()<<std::endl;
}
root_part->deleteTree();
}
while (count++ < nevent);
for (int j = 0;j < n; j++)
{
std::cout << "[" << q2low[j] << ".." << q2high[j] << "] = " << counts[j] << std::endl;
}
file->Write();
file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runVectorIsr(int nevent, EvtGen &myGenerator) {
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
TFile *file=new TFile("vectorisr.root", "RECREATE");
TH1F* cosv = new TH1F("h1","Cos vector in e+e- frame",
50,-1.0,1.0);
TH1F* cosd1 = new TH1F("h2","Cos helang 1st dau of vector",
50,-1.0,1.0);
TH1F* cosd1d1 = new TH1F("h3","Cos helang 1st dau of 1st dau",
50,-1.0,1.0);
TH1F* cosd2d1 = new TH1F("h4","Cos helang 1st dau of 2nd dau",
50,-1.0,1.0);
TH2F* d1vsd1d1 = new TH2F("h5","Cos helangs d1 vs d1d1",
20,-1.0,1.0,20,-1.0,1.0);
TH2F* d2vsd2d1 = new TH2F("h6","Cos helangs d2 vs d2d1",
20,-1.0,1.0,20,-1.0,1.0);
TH2F* d1d1vsd2d1 = new TH2F("h7","Cos helangs d1d1 vs d2d1",
20,-1.0,1.0,20,-1.0,1.0);
TH1F* chidd = new TH1F("h8","Chi - angle between decay planes",
60,0.,360.0);
TH2F* chi12vsd1d1 = new TH2F("h9","Chi 1-2 vs d1d1",
30,0.,360.0,20,-1.0,1.0);
TH2F* chi12vsd2d1 = new TH2F("h10","Chi 1-2 vs d2d1",
30,0.,360.0,20,-1.0,1.0);
TH2F* chi21vsd1d1 = new TH2F("h11","Chi 2-1 vs d1d1",
30,0.,360.0,20,-1.0,1.0);
TH2F* chi21vsd2d1 = new TH2F("h12","Chi 2-1 vs d2d1",
30,0.,360.0,20,-1.0,1.0);
int count=1;
char udecay_name[100];
EvtVector4R cm, v, d1, d2, d1d1, d1d2, d2d1, d2d2;
strcpy(udecay_name,"exampleFiles/VECTORISR.DEC");
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(VPHO,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
cm=root_part->getP4Lab();
v=root_part->getDaug(0)->getP4Lab();
d1=root_part->getDaug(0)->getDaug(0)->getP4Lab();
d2=root_part->getDaug(0)->getDaug(1)->getP4Lab();
cosv->Fill(v.get(3)/v.d3mag());
double cosdecayd1=EvtDecayAngle(cm,v,d1);
double cosdecayd2=EvtDecayAngle(cm,v,d2);
cosd1->Fill(cosdecayd1);
// now get daughters of the daughters
//
// first daughter of first daughter
if (root_part->getDaug(0)->getDaug(0)->getNDaug()>=2){
d1d1=root_part->getDaug(0)->getDaug(0)->getDaug(0)->getP4Lab();
double cosdecayd1d1=EvtDecayAngle(v,d1,d1d1);
cosd1d1->Fill(cosdecayd1d1);
d1vsd1d1->Fill(cosdecayd1,cosdecayd1d1,1.0);
}
// first daughter of second daughter
if (root_part->getDaug(0)->getDaug(1)->getNDaug()>=2){
d2d1=root_part->getDaug(0)->getDaug(1)->getDaug(0)->getP4Lab();
double cosdecayd2d1=EvtDecayAngle(v,d2,d2d1);
cosd2d1->Fill(cosdecayd2d1);
d2vsd2d1->Fill(cosdecayd2,cosdecayd2d1,1.0);
if (root_part->getDaug(0)->getDaug(0)->getNDaug()>=2){
d1d1=root_part->getDaug(0)->getDaug(0)->getDaug(0)->getP4Lab();
double cosdecayd1d1=EvtDecayAngle(v,d1,d1d1);
d1d1vsd2d1->Fill(cosdecayd1d1,cosdecayd2d1,1.0);
//second daughters of daughters 1 and 2
d1d2=root_part->getDaug(0)->getDaug(0)->getDaug(1)->getP4Lab();
d2d2=root_part->getDaug(0)->getDaug(1)->getDaug(1)->getP4Lab();
double chi21=57.29578*EvtDecayAngleChi(v,d2d1,d2d2,d1d1,d1d2);
double chi12=57.29578*EvtDecayAngleChi(v,d1d1,d1d2,d2d1,d2d2);
chidd->Fill(chi12);
chi12vsd1d1->Fill(chi12,cosdecayd1d1,1.0);
chi12vsd2d1->Fill(chi12,cosdecayd2d1,1.0);
chi21vsd1d1->Fill(chi21,cosdecayd1d1,1.0);
chi21vsd2d1->Fill(chi21,cosdecayd2d1,1.0);
}
}
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runBsquark(int nevent, EvtGen &myGenerator) {
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
TFile *file=new TFile("bsquark.root", "RECREATE");
TH1F* elep = new TH1F("h1","Elep",50,0.0,1.5);
TH1F* q2 = new TH1F("h2","q2",50,0.0,3.0);
TH2F* dalitz=new TH2F("h3","q2 vs. Elep",50,0.0,1.5,50,0.0,3.0);
TH1F* elepbar = new TH1F("h11","Elep bar",
50,0.0,1.5);
TH1F* q2bar = new TH1F("h12","q2 bar",
50,0.0,3.0);
TH2F* dalitzbar=new TH2F("h13","q2 vs. Elep bar",50,0.0,1.5,50,0.0,3.0);
int count=1;
char udecay_name[100];
EvtVector4R cm, v, d1, d2, d1d1, d1d2, d2d1, d2d2;
strcpy(udecay_name,"exampleFiles/BSQUARK.DEC");
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(10.55,0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(VPHO,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle* p=root_part->nextIter();
while (p) {
if (p->getId()==B0) {
EvtParticle *dstar=p->getDaug(0);
EvtParticle *lepton=p->getDaug(1);
EvtParticle *sneutrino=p->getDaug(2);
EvtVector4R p4dstar=dstar->getP4();
EvtVector4R p4lepton=lepton->getP4();
EvtVector4R p4sneutrino=sneutrino->getP4();
elep->Fill(p4lepton.get(0));
q2->Fill((p4lepton+p4sneutrino).mass2());
dalitz->Fill(p4lepton.get(0),(p4lepton+p4sneutrino).mass2(),1.0);
}
if (p->getId()==B0B) {
EvtParticle *dstar=p->getDaug(0);
EvtParticle *lepton=p->getDaug(1);
EvtParticle *sneutrino=p->getDaug(2);
EvtVector4R p4dstar=dstar->getP4();
EvtVector4R p4lepton=lepton->getP4();
EvtVector4R p4sneutrino=sneutrino->getP4();
elepbar->Fill(p4lepton.get(0));
q2bar->Fill((p4lepton+p4sneutrino).mass2());
dalitzbar->Fill(p4lepton.get(0),(p4lepton+p4sneutrino).mass2(),1.0);
}
p=p->nextIter();
}
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runK3gamma(int nevent, EvtGen &myGenerator) {
static EvtId B0=EvtPDL::getId(std::string("B0"));
TFile *file=new TFile("k3gamma.root", "RECREATE");
TH1F* costheta = new TH1F("h1","cosTheta", 100,-1.0,1.0);
int count=1;
char udecay_name[100];
EvtVector4R cm, v, d1, d2, d1d1, d1d2, d2d1, d2d2;
strcpy(udecay_name,"exampleFiles/K3GAMMA.DEC");
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(B0,p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *k3=root_part->getDaug(0);
EvtParticle *k=k3->getDaug(0);
EvtVector4R p4b=root_part->getP4Lab();
EvtVector4R p4k3=k3->getP4Lab();
EvtVector4R p4k=k->getP4Lab();
costheta->Fill(EvtDecayAngle(p4b,p4k3,p4k));
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runLambda(int nevent, EvtGen &myGenerator) {
static EvtId LAMBDA=EvtPDL::getId(std::string("Lambda0"));
TFile *file=new TFile("lambda.root", "RECREATE");
TH1F* costheta = new TH1F("h1","cosTheta",100,-1.0,1.0);
int count=1;
char udecay_name[100];
EvtVector4R cm, v, d1, d2, d1d1, d1d2, d2d1, d2d2;
strcpy(udecay_name,"exampleFiles/LAMBDA.DEC");
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMass(LAMBDA),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(LAMBDA,p_init);
EvtSpinDensity rho;
rho.setDim(2);
rho.set(0,0,1.0);
rho.set(0,1,0.0);
rho.set(1,0,0.0);
rho.set(1,1,0.0);
root_part->setSpinDensityForwardHelicityBasis(rho);
myGenerator.generateDecay(root_part);
EvtParticle *p=root_part->getDaug(0);
EvtVector4R p4lambda=root_part->getP4Lab();
EvtVector4R p4p=p->getP4Lab();
costheta->Fill(p4p.get(3)/p4p.d3mag());
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runTauTauPiPi(int nevent, EvtGen &myGenerator) {
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
TFile *file=new TFile("tautaupipi.root", "RECREATE");
TH1F* cospi1 = new TH1F("h1","cos theta pi1", 50,-1.0,1.0);
TH1F* cospi2 = new TH1F("h2","cos theta pi2", 50,-1.0,1.0);
TH1F* costheta = new TH1F("h3","cos theta", 50,-1.0,1.0);
std::ofstream outmix;
outmix.open("tautaupipi.dat");
int count=1;
EvtVector4R tau1,tau2,pi1,pi2;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/TAUTAUPIPI.DEC");
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(VPHO,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
tau1=root_part->getDaug(0)->getP4Lab();
tau2=root_part->getDaug(1)->getP4Lab();
pi1=root_part->getDaug(0)->getDaug(0)->getP4Lab();
pi2=root_part->getDaug(1)->getDaug(0)->getP4Lab();
cospi1->Fill( EvtDecayAngle(tau1+tau2,tau1,pi1) );
cospi2->Fill( EvtDecayAngle(tau1+tau2,tau2,pi2) );
costheta->Fill(tau1.get(3)/tau1.d3mag());
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runTauTauEE(int nevent, EvtGen &myGenerator) {
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
TFile *file=new TFile("tautauee.root", "RECREATE");
TH1F* e1 = new TH1F("h1","e1",55,0.0,5.5);
TH1F* e2 = new TH1F("h2","e2",55,0.0,5.5);
TH2F* e1vse2 = new TH2F("h3","e1 vs e2",55,0.0,5.5,
55,0.0,5.5);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/TAUTAUEE.DEC");
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(VPHO, p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
e1->Fill(root_part->getDaug(0)->getDaug(0)->getP4Lab().get(0));
e2->Fill(root_part->getDaug(1)->getDaug(0)->getP4Lab().get(0));
e1vse2->Fill(root_part->getDaug(0)->getDaug(0)->getP4Lab().get(0),
root_part->getDaug(1)->getDaug(0)->getP4Lab().get(0),1.0);
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runTauTau2Pi2Pi(int nevent, EvtGen &myGenerator) {
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
TFile *file=new TFile("tautau2pi2pi.root", "RECREATE");
TH1F* e1 = new TH1F("h1","mrho",200,0.0,2.0);
TH1F* e2 = new TH1F("h2","coshel",200,-1.0,1.0);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/TAUTAU2PI2PI.DEC");
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(VPHO, p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtVector4R p4tau =root_part->getDaug(0)->getP4();
EvtVector4R p4rho =root_part->getDaug(0)->getDaug(0)->getP4()
+ root_part->getDaug(0)->getDaug(1)->getP4();
EvtVector4R p4pi = root_part->getDaug(0)->getDaug(0)->getP4();
e1->Fill(p4rho.mass());
double dcostheta=EvtDecayAngle(p4tau,p4rho,p4pi);
e2->Fill(dcostheta);
p4tau =root_part->getDaug(1)->getP4();
p4rho =root_part->getDaug(1)->getDaug(0)->getP4()
+ root_part->getDaug(1)->getDaug(1)->getP4();
p4pi = root_part->getDaug(1)->getDaug(0)->getP4();
e1->Fill(p4rho.mass());
dcostheta=EvtDecayAngle(p4tau,p4rho,p4pi);
e2->Fill(dcostheta);
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runTauTau3Pi3Pi(int nevent, EvtGen &myGenerator) {
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
TFile *file=new TFile("tautau3pi3pi.root", "RECREATE");
TH1F* e1 = new TH1F("h1","a1",200,0.0,2.0);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/TAUTAU3PI3PI.DEC");
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(VPHO, p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtVector4R p4tau =root_part->getDaug(0)->getP4();
EvtVector4R p4a1 =root_part->getDaug(0)->getDaug(0)->getP4()
+ root_part->getDaug(0)->getDaug(1)->getP4()
+ root_part->getDaug(0)->getDaug(2)->getP4();
e1->Fill(p4a1.mass());
p4tau =root_part->getDaug(1)->getP4();
p4a1 =root_part->getDaug(1)->getDaug(0)->getP4()
+ root_part->getDaug(1)->getDaug(1)->getP4()
+ root_part->getDaug(1)->getDaug(2)->getP4();
e1->Fill(p4a1.mass());
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runJPsiKstar(int nevent, EvtGen &myGenerator, int modeInt) {
std::ofstream outmix;
outmix.open("jpsikstar.dat");
int count=1;
char udecay_name[100];
if (modeInt==0) strcpy(udecay_name,"exampleFiles/JPSIKSTAR.DEC");
if (modeInt==1) strcpy(udecay_name,"exampleFiles/JPSIKSTAR1.DEC");
if (modeInt==2) strcpy(udecay_name,"exampleFiles/JPSIKSTAR2.DEC");
if (modeInt==3) strcpy(udecay_name,"exampleFiles/JPSIKSTAR3.DEC");
if (modeInt==4) strcpy(udecay_name,"exampleFiles/JPSIKSTAR4.DEC");
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
static EvtId GAMM=EvtPDL::getId(std::string("gamma"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
myGenerator.readUDecay(udecay_name);
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *btag,*bcp;
if (root_part->getDaug(0)->getNDaug()==3){
btag=root_part->getDaug(0);
bcp=root_part->getDaug(1);
}
else{
bcp=root_part->getDaug(0);
btag=root_part->getDaug(1);
}
EvtId tag;
if (btag->getId()==B0B){
tag=B0;
}
else{
tag=B0B;
}
EvtParticle *p_b,*p_psi,*p_kstar,*p_pi0,*p_kz,*p_ep,*p_em;
EvtVector4R p4_b,p4_psi,p4_kstar,p4_pi0,p4_kz,p4_ep,p4_em;
p_b=bcp;
p_psi=p_b->getDaug(0);
p_kstar=p_b->getDaug(1);
p_pi0=p_kstar->getDaug(0);
p_kz=p_kstar->getDaug(1);
p_ep=p_psi->getDaug(0);
p_em=p_psi->getDaug(1);
p4_b=p_b->getP4Lab();
p4_psi=p_psi->getP4Lab();
p4_kstar=p_kstar->getP4Lab();
p4_pi0=p_pi0->getP4Lab();
p4_kz=p_kz->getP4Lab();
p4_ep=p_ep->getP4Lab();
p4_em=p_em->getP4Lab();
outmix << tag.getId() << " ";
outmix << root_part->getDaug(0)->getLifetime() << " ";
outmix << root_part->getDaug(1)->getLifetime() << " ";
outmix << EvtDecayAngle(p4_b,p4_ep+p4_em,p4_ep) << " " ;
outmix << EvtDecayAngle(p4_b,p4_pi0+p4_kz,p4_pi0) << " " ;
outmix << EvtDecayAngleChi(p4_b,p4_pi0,p4_kz,
p4_ep, p4_em ) << "\n" ;
root_part->deleteTree();
}while (count++<10000);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runSVVCPLH(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("svvcplh.root", "RECREATE");
TH1F* t = new TH1F("h1","t",50,0.0,5.0);
TH1F* cospsi = new TH1F("h2","cos theta e+",50,-1.0,1.0);
TH1F* cosphi = new TH1F("h3","cos theta k+",50,-1.0,1.0);
TH1F* chi = new TH1F("h4","chi",50,0.0,2.0*EvtConst::pi);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/SVVCPLH.DEC");
myGenerator.readUDecay(udecay_name);
static EvtId BS0=EvtPDL::getId(std::string("B_s0"));
static EvtId BSB=EvtPDL::getId(std::string("anti-B_s0"));
do{
EvtVector4R p_init(EvtPDL::getMass(BS0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(BS0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p_b,*p_psi,*p_phi,*p_kp,*p_km,*p_ep,*p_em;
EvtVector4R p4_b,p4_psi,p4_phi,p4_kp,p4_km,p4_ep,p4_em;
p_b=root_part;
if (p_b->getNDaug()==1) p_b=p_b->getDaug(0);
p_psi=p_b->getDaug(0);
p_phi=p_b->getDaug(1);
p_kp=p_phi->getDaug(0);
p_km=p_phi->getDaug(1);
p_ep=p_psi->getDaug(0);
p_em=p_psi->getDaug(1);
p4_b=p_b->getP4Lab();
p4_psi=p_psi->getP4Lab();
p4_phi=p_phi->getP4Lab();
p4_kp=p_kp->getP4Lab();
p4_km=p_km->getP4Lab();
p4_ep=p_ep->getP4Lab();
p4_em=p_em->getP4Lab();
t->Fill(root_part->getLifetime());
cospsi->Fill(EvtDecayAngle(p4_b,p4_ep+p4_em,p4_ep));
cosphi->Fill(EvtDecayAngle(p4_b,p4_kp+p4_km,p4_kp));
chi->Fill(EvtDecayAngleChi(p4_b,p4_kp,p4_km,
p4_ep, p4_em ));
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runSVSCPLH(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("svscplh.root", "RECREATE");
TH1F* t = new TH1F("h1","t",200,-5.0,5.0);
TH1F* tB0tag = new TH1F("h2","dt B0 tag (ps)",200,-15.0,15.0);
TH1F* tB0Btag = new TH1F("h3","dt B0B tag (ps)",200,-15.0,15.0);
TH1F* ctheta = new TH1F("h4","costheta",50,-1.0,1.0);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/SVSCPLH.DEC");
myGenerator.readUDecay(udecay_name);
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
std::ofstream outmix;
outmix.open("svscplh.dat");
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p_tag,*p_cp,*p_jpsi,*p_ep;
EvtVector4R p4_tag,p4_cp,p4_jpsi,p4_ep;
p_tag=root_part->getDaug(0);
p_cp=root_part->getDaug(1);
p_jpsi=p_cp->getDaug(0);
p_ep=p_jpsi->getDaug(0);
p4_tag=p_tag->getP4Lab();
p4_cp=p_cp->getP4Lab();
p4_jpsi=p_jpsi->getP4Lab();
p4_ep=p_ep->getP4Lab();
double dt=p_cp->getLifetime()-p_tag->getLifetime();
dt=dt/(1e-12*3e11);
t->Fill(dt);
if (p_tag->getId()==B0){
tB0tag->Fill(dt);
outmix << dt << " 1"<<std::endl;
}
if (p_tag->getId()==B0B){
tB0Btag->Fill(dt);
outmix << dt << " -1"<<std::endl;
}
ctheta->Fill(EvtDecayAngle(p4_cp,p4_jpsi,p4_ep));
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runSSDCP(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("ssdcp.root", "RECREATE");
TH1F* t = new TH1F("h1","dt",100,-15.0,15.0);
TH1F* tB0tag = new TH1F("h2","dt B0 tag (ps)",100,-15.0,15.0);
TH1F* tB0Btag = new TH1F("h3","dt B0B tag (ps)",100,-15.0,15.0);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/SSDCP.DEC");
myGenerator.readUDecay(udecay_name);
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
std::ofstream outmix;
do{
EvtVector4R pinit(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,pinit);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p_tag,*p_cp,*p_jpsi;
EvtVector4R p4_tag,p4_cp,p4_jpsi,p4_ep;
p_tag=root_part->getDaug(0);
p_cp=root_part->getDaug(1);
p_jpsi=p_cp->getDaug(0);
//p_ep=p_jpsi->getDaug(0);
p4_tag=p_tag->getP4Lab();
p4_cp=p_cp->getP4Lab();
p4_jpsi=p_jpsi->getP4Lab();
//p4_ep=p_ep->getP4Lab();
double dt=p_cp->getLifetime()-p_tag->getLifetime();
dt=dt/(1e-12*EvtConst::c);
t->Fill(dt);
if (p_tag->getId()==B0){
tB0tag->Fill(dt);
}
if (p_tag->getId()==B0B){
tB0Btag->Fill(dt);
}
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runKstarstargamma(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("kstarstargamma.root", "RECREATE");
TH1F* m = new TH1F("h1","mkpi",100,0.5,2.5);
TH1F* ctheta = new TH1F("h2","ctheta",100,-1.0,1.0);
int count=1;
myGenerator.readUDecay("exampleFiles/KSTARSTARGAMMA.DEC");
static EvtId B0=EvtPDL::getId(std::string("B0"));
std::ofstream outmix;
do{
EvtVector4R pinit(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(B0,pinit);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p_kaon,*p_pion;
EvtVector4R p4_kaon,p4_pion;
p_kaon=root_part->getDaug(0);
p_pion=root_part->getDaug(1);
p4_kaon=p_kaon->getP4Lab();
p4_pion=p_pion->getP4Lab();
m->Fill((p4_kaon+p4_pion).mass());
ctheta->Fill(EvtDecayAngle(pinit,p4_kaon+p4_pion,p4_kaon));
//EvtGenReport(EVTGEN_INFO,"EvtGen") << "ctheta:"<<EvtDecayAngle(pinit,p4_kaon+p4_pion,p4_kaon)<<std::endl;
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runDSTARPI(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("dstarpi.root", "RECREATE");
TH1F* t = new TH1F("h1","dt",100,-15.0,15.0);
TH1F* tB0tagpip = new TH1F("h2","dt B0 tag pi+ (ps)",100,-15.0,15.0);
TH1F* tB0Btagpip = new TH1F("h3","dt B0B tag pi+(ps)",100,-15.0,15.0);
TH1F* tB0tagpim = new TH1F("h4","dt B0 tag pi- (ps)",100,-15.0,15.0);
TH1F* tB0Btagpim = new TH1F("h5","dt B0B tag pi- (ps)",100,-15.0,15.0);
int count=1;
myGenerator.readUDecay("exampleFiles/DSTARPI.DEC");
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
static EvtId PIP=EvtPDL::getId(std::string("pi+"));
static EvtId PIM=EvtPDL::getId(std::string("pi-"));
std::ofstream outmix;
do{
EvtVector4R pinit(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,pinit);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p_tag,*p_cp,*p_dstar,*p_pi;
p_tag=root_part->getDaug(0);
p_cp=root_part->getDaug(1);
p_dstar=p_cp->getDaug(1);
p_pi=p_cp->getDaug(0);
double dt=p_cp->getLifetime()-p_tag->getLifetime();
dt=dt/(1e-12*EvtConst::c);
t->Fill(dt);
if (p_tag->getId()==B0){
if (p_pi->getId()==PIP) tB0tagpip->Fill(dt);
if (p_pi->getId()==PIM) tB0tagpim->Fill(dt);
}
if (p_tag->getId()==B0B){
if (p_pi->getId()==PIP) tB0Btagpip->Fill(dt);
if (p_pi->getId()==PIM) tB0Btagpim->Fill(dt);
}
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runETACPHIPHI(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("etacphiphi.root", "RECREATE");
TH2F* cosphi12 = new TH2F("h1","cos phi1 vs phi2",
50,-1.0,1.0,50,-1.0,1.0);
TH1F* cosphi1 = new TH1F("h2","cos phi1",50,-1.0,1.0);
TH1F* cosphi2 = new TH1F("h3","cos phi2",50,-1.0,1.0);
TH1F* chi = new TH1F("h4","chi",50,0.0,2.0*EvtConst::pi);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/ETACPHIPHI.DEC");
myGenerator.readUDecay(udecay_name);
static EvtId ETAC=EvtPDL::getId(std::string("eta_c"));
do{
EvtVector4R p_init(EvtPDL::getMass(ETAC),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(ETAC,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p_etac,*p_phi1,*p_phi2,*p_kp1,*p_km1,*p_kp2,*p_km2;
EvtVector4R p4_etac,p4_phi1,p4_phi2,p4_kp1,p4_km1,p4_kp2,p4_km2;
p_etac=root_part;
p_phi1=p_etac->getDaug(0);
p_phi2=p_etac->getDaug(1);
p_kp1=p_phi1->getDaug(0);
p_km1=p_phi1->getDaug(1);
p_kp2=p_phi2->getDaug(0);
p_km2=p_phi2->getDaug(1);
p4_etac=p_etac->getP4Lab();
p4_phi1=p_phi1->getP4Lab();
p4_phi2=p_phi2->getP4Lab();
p4_kp1=p_kp1->getP4Lab();
p4_km1=p_km1->getP4Lab();
p4_kp2=p_kp2->getP4Lab();
p4_km2=p_km2->getP4Lab();
cosphi12->Fill(EvtDecayAngle(p4_etac,p4_phi1,p4_kp1),
EvtDecayAngle(p4_etac,p4_phi2,p4_kp2),1.0);
cosphi1->Fill(EvtDecayAngle(p4_etac,p4_phi1,p4_kp1));
cosphi2->Fill(EvtDecayAngle(p4_etac,p4_phi2,p4_kp2));
chi->Fill(EvtDecayAngleChi(p4_etac,p4_kp1,p4_km1,
p4_kp2, p4_km2));
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runVVPiPi(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("vvpipi.root", "RECREATE");
TH1F* cospsi = new TH1F("h1","cos theta J/psi ",50,-1.0,1.0);
TH1F* cose = new TH1F("h2","cos theta e+ ",50,-1.0,1.0);
TH1F* mpipi = new TH1F("h3","m pipi ",50,0.0,1.0);
TH2F* cosevspsi = new TH2F("h4","cos theta e+vs cos thete J/psi ",
25,-1.0,1.0,25,-1.0,1.0);
TH1F* cose1 = new TH1F("h5","cos theta e+ 1 ",50,-1.0,1.0);
TH1F* cose2 = new TH1F("h6","cos theta e+ 2 ",50,-1.0,1.0);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/VVPIPI.DEC");
myGenerator.readUDecay(udecay_name);
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(B0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p_b,*p_psip,*p_psi,*p_ep,*p_pi1,*p_pi2;
EvtVector4R p4_b,p4_psip,p4_psi,p4_ep,p4_pi1,p4_pi2;
p_b=root_part;
p_psip=p_b->getDaug(0);
p_psi=p_psip->getDaug(0);
p_pi1=p_psip->getDaug(1);
p_pi2=p_psip->getDaug(2);
p_ep=p_psi->getDaug(0);
p4_b=p_b->getP4Lab();
p4_psip=p_psip->getP4Lab();
p4_psi=p_psi->getP4Lab();
p4_pi1=p_pi1->getP4Lab();
p4_pi2=p_pi2->getP4Lab();
p4_ep=p_ep->getP4Lab();
cospsi->Fill(EvtDecayAngle(p4_b,p4_psip,p4_psi));
cose->Fill(EvtDecayAngle(p4_psip,p4_psi,p4_ep));
mpipi->Fill((p4_pi1+p4_pi2).mass());
cosevspsi->Fill(EvtDecayAngle(p4_b,p4_psip,p4_psi),
EvtDecayAngle(p4_psip,p4_psi,p4_ep),1.0);
if (std::fabs(EvtDecayAngle(p4_b,p4_psip,p4_psi))>0.95){
cose1->Fill(EvtDecayAngle(p4_psip,p4_psi,p4_ep));
}
if (std::fabs(EvtDecayAngle(p4_b,p4_psip,p4_psi))<0.05){
cose2->Fill(EvtDecayAngle(p4_psip,p4_psi,p4_ep));
}
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runSVVHelAmp(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("svvhelamp.root", "RECREATE");
TH1F* cospip = new TH1F("h1","cos theta pi+",
50,-1.0,1.0);
TH1F* cospim = new TH1F("h2","cos theta pi-",
50,-1.0,1.0);
TH1F* chi = new TH1F("h3","chi pi+ to pi- in D+ direction",
50,0.0,EvtConst::twoPi);
TH1F* chicospipp = new TH1F("h4","chi pi+ to pi- in D+ direction (cospip>0)",
50,0.0,EvtConst::twoPi);
TH1F* chicospipn = new TH1F("h5","chi pi+ to pi- in D+ direction (cospip<0",
50,0.0,EvtConst::twoPi);
TH1F* chipp = new TH1F("h6","chi pi+ to pi- in D+ direction (cospip>0,cospim>0)",
50,0.0,EvtConst::twoPi);
TH1F* chipn = new TH1F("h7","chi pi+ to pi- in D+ direction (cospip>0,cospim<0)",
50,0.0,EvtConst::twoPi);
TH1F* chinp = new TH1F("h8","chi pi+ to pi- in D+ direction (cospip<0,cospim>0)",
50,0.0,EvtConst::twoPi);
TH1F* chinn = new TH1F("h9","chi pi+ to pi- in D+ direction (cospip<0,cospim<0)",
50,0.0,EvtConst::twoPi);
TH1F* chinnnn = new TH1F("h10","chi pi+ to pi- in D+ direction (cospip<-0.5,cospim<-0.5)",
50,0.0,EvtConst::twoPi);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/SVVHELAMP.DEC");
myGenerator.readUDecay(udecay_name);
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(B0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p_b,*p_dstp,*p_dstm,*p_pip,*p_pim,*p_d0,*p_d0b;
EvtVector4R p4_b,p4_dstp,p4_dstm,p4_pip,p4_pim,p4_d0,p4_d0b;
p_b=root_part;
p_dstp=p_b->getDaug(0);
p_dstm=p_b->getDaug(1);
p_pip=p_dstp->getDaug(1);
p_pim=p_dstm->getDaug(1);
p_d0=p_dstp->getDaug(0);
p_d0b=p_dstm->getDaug(0);
p4_b=p_b->getP4Lab();
p4_dstp=p_dstp->getP4Lab();
p4_dstm=p_dstm->getP4Lab();
p4_pip=p_pip->getP4Lab();
p4_pim=p_pim->getP4Lab();
p4_d0=p_d0->getP4Lab();
p4_d0b=p_d0b->getP4Lab();
double costhpip=EvtDecayAngle(p4_b,p4_pip+p4_d0,p4_pip);
double costhpim=EvtDecayAngle(p4_b,p4_pim+p4_d0b,p4_pim);
double chiang=EvtDecayAngleChi(p4_b,p4_pip,p4_d0,p4_pim,p4_d0b);
cospip->Fill( costhpip );
cospim->Fill( costhpim );
chi->Fill( chiang );
if (costhpip>0) chicospipp->Fill( chiang );
if (costhpip<0) chicospipn->Fill( chiang );
if (costhpip>0 && costhpim>0) chipp->Fill( chiang );
if (costhpip>0 && costhpim<0) chipn->Fill( chiang );
if (costhpip<0 && costhpim>0) chinp->Fill( chiang );
if (costhpip<0 && costhpim<0) chinn->Fill( chiang );
if (costhpip<-0.5 && costhpim<-0.5) chinnnn->Fill( chiang );
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runPartWave(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("partwave.root", "RECREATE");
TH1F* cospip = new TH1F("h1","cos theta pi+",
50,-1.0,1.0);
TH1F* cospim = new TH1F("h2","cos theta pi-",
50,-1.0,1.0);
TH1F* chi = new TH1F("h3","chi pi+ to pi- in D+ direction",
50,0.0,EvtConst::twoPi);
TH1F* chicospipp = new TH1F("h4","chi pi+ to pi- in D+ direction (cospip>0)",
50,0.0,EvtConst::twoPi);
TH1F* chicospipn = new TH1F("h5","chi pi+ to pi- in D+ direction (cospip<0",
50,0.0,EvtConst::twoPi);
TH1F* chipp = new TH1F("h6","chi pi+ to pi- in D+ direction (cospip>0,cospim>0)",
50,0.0,EvtConst::twoPi);
TH1F* chipn = new TH1F("h7","chi pi+ to pi- in D+ direction (cospip>0,cospim<0)",
50,0.0,EvtConst::twoPi);
TH1F* chinp = new TH1F("h8","chi pi+ to pi- in D+ direction (cospip<0,cospim>0)",
50,0.0,EvtConst::twoPi);
TH1F* chinn = new TH1F("h9","chi pi+ to pi- in D+ direction (cospip<0,cospim<0)",
50,0.0,EvtConst::twoPi);
TH1F* chinnnn = new TH1F("h10","chi pi+ to pi- in D+ direction (cospip<-0.5,cospim<-0.5)",
50,0.0,EvtConst::twoPi);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/PARTWAVE.DEC");
myGenerator.readUDecay(udecay_name);
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(B0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p_b,*p_dstp,*p_dstm,*p_pip,*p_pim,*p_d0,*p_d0b;
EvtVector4R p4_b,p4_dstp,p4_dstm,p4_pip,p4_pim,p4_d0,p4_d0b;
p_b=root_part;
p_dstp=p_b->getDaug(0);
p_dstm=p_b->getDaug(1);
p_pip=p_dstp->getDaug(1);
p_pim=p_dstm->getDaug(1);
p_d0=p_dstp->getDaug(0);
p_d0b=p_dstm->getDaug(0);
p4_b=p_b->getP4Lab();
p4_dstp=p_dstp->getP4Lab();
p4_dstm=p_dstm->getP4Lab();
p4_pip=p_pip->getP4Lab();
p4_pim=p_pim->getP4Lab();
p4_d0=p_d0->getP4Lab();
p4_d0b=p_d0b->getP4Lab();
double costhpip=EvtDecayAngle(p4_b,p4_pip+p4_d0,p4_pip);
double costhpim=EvtDecayAngle(p4_b,p4_pim+p4_d0b,p4_pim);
double chiang=EvtDecayAngleChi(p4_b,p4_pip,p4_d0,p4_pim,p4_d0b);
cospip->Fill( costhpip );
cospim->Fill( costhpim );
chi->Fill( chiang );
if (costhpip>0) chicospipp->Fill( chiang );
if (costhpip<0) chicospipn->Fill( chiang );
if (costhpip>0 && costhpim>0) chipp->Fill( chiang );
if (costhpip>0 && costhpim<0) chipn->Fill( chiang );
if (costhpip<0 && costhpim>0) chinp->Fill( chiang );
if (costhpip<0 && costhpim<0) chinn->Fill( chiang );
if (costhpip<-0.5 && costhpim<-0.5) chinnnn->Fill( chiang );
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runPartWave2(int nevent, EvtGen &myGenerator) {
TFile file("partwave2.root", "RECREATE");
TH1F* cthetapi = new TH1F("h1","cos theta pi",
50,-1.0,1.0);
TH1F* cthetapi2 = new TH1F("h2","cos theta pi (|cosrho|<0.1)",
50,-1.0,1.0);
TH1F* cthetan = new TH1F("h3","cos thetan",
50,-1.0,1.0);
//TH1F* cthetan2 = new TH1F("h4","cos thetan costhetapi>0 ",
// 50,-1.0,1.0);
TH1F* cthetarho = new TH1F("h4","cos thetarho ",
50,-1.0,1.0);
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/PARTWAVE2.DEC");
myGenerator.readUDecay(udecay_name);
static EvtId JPSI=EvtPDL::getId(std::string("J/psi"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(B0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p_b,*p_jpsi,*p_rho,*p_pi1,*p_pi2;
EvtVector4R p4_b,p4_jpsi,p4_rho,p4_pi1,p4_pi2;
p_b=root_part;
p_jpsi=root_part->getDaug(0);
p_rho=0;
if (p_jpsi->getDaug(0)->getNDaug()==2){
p_rho=p_jpsi->getDaug(0);
}
if (p_jpsi->getDaug(1)->getNDaug()==2){
p_rho=p_jpsi->getDaug(1);
}
assert(p_rho!=0);
p_pi1=p_rho->getDaug(0);
p_pi2=p_rho->getDaug(1);
p4_b=p_b->getP4Lab();
p4_jpsi=p_jpsi->getP4Lab();
p4_rho=p_rho->getP4Lab();
p4_pi1=p_pi1->getP4Lab();
p4_pi2=p_pi2->getP4Lab();
double costhetan=EvtDecayPlaneNormalAngle(p4_b,p4_jpsi,p4_pi1,p4_pi2);
//EvtGenReport(EVTGEN_INFO,"EvtGen") << "costhetan:"<<costhetan<<std::endl;
cthetan->Fill( costhetan );
double costhpi=EvtDecayAngle(p4_jpsi,p4_rho,p4_pi1);
double costhrho=EvtDecayAngle(p4_b,p4_jpsi,p4_rho);
//EvtGenReport(EVTGEN_INFO,"EvtGen") << "costhetarho:"<<costhrho<<std::endl;
cthetarho->Fill( costhrho );
//if (((p4_rho.get(3)/p4_rho.d3mag()))<-0.95) cthetan2->Fill( costhetan );
cthetapi->Fill( costhpi );
if ((p4_rho.get(3)/p4_rho.d3mag())>0.9) {
cthetapi2->Fill( costhpi );
}
root_part->deleteTree();
}while (count++<nevent);
file.Write(); file.Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runTwoBody(int nevent, EvtGen &myGenerator, std::string decFile,
std::string rootFile) {
TFile *file=new TFile(rootFile.c_str(), "RECREATE");
int count=0;
myGenerator.readUDecay(decFile.c_str());
static EvtId B0=EvtPDL::getId(std::string("B0"));
vector<TH1F*> histograms;
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(B0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
//root_part->printTree();
myGenerator.generateDecay(root_part);
int nhist=0;
EvtParticle* p=root_part;
do {
int nDaug=p->getNDaug();
if (!(nDaug==0||nDaug==2)) {
EvtGenReport(EVTGEN_INFO,"EvtGen") << "nDaug="<<nDaug<<" but can only handle 0 or 2!"<<std::endl;
abort();
}
if (nDaug==2){
if (p->getParent()==0){
EvtVector4R p4=p->getDaug(0)->getP4();
double ctheta=p4.get(3)/p4.d3mag();
double phi=atan2(p4.get(2),p4.get(1));
if (count==0){
histograms.push_back(new TH1F("h1","cos theta",50,-1.0,1.0));
histograms.push_back(new TH1F("h2","phi",50,-EvtConst::pi
,EvtConst::pi));
}
histograms[nhist++]->Fill(ctheta);
histograms[nhist++]->Fill(phi);
}
else{
double ctheta=EvtDecayAngle(p->getParent()->getP4Lab(),
p->getP4Lab(),
p->getDaug(0)->getP4Lab());
if (count==0){
// char* tmp=new char[10];
// std::ostrstream strm(tmp,9);
// strm << (nhist+1) << '\0'<< std::endl;
// histograms.push_back(new TH1F(TString("h")+tmp,TString("cos theta")+tmp,50,-1.0,1.0));
std::ostringstream strm;
strm << (nhist+1);
- histograms.push_back( new TH1F(TString("h") + strm.str(),
- TString("cos theta") + strm.str(),
+ histograms.push_back( new TH1F(TString("h") + strm.str().c_str(),
+ TString("cos theta") + strm.str().c_str(),
50,-1.0,1.0) );
}
histograms[nhist++]->Fill(ctheta);
if (p->getDaug(0)->getNDaug()==2) {
double costhetan=EvtDecayPlaneNormalAngle(p->getParent()->getP4Lab(),
p->getP4Lab(),
p->getDaug(0)->getDaug(0)->getP4Lab(),
p->getDaug(0)->getDaug(1)->getP4Lab());
if (count==0){
// char* tmp=new char[10];
// std::ostrstream strm(tmp,9);
// strm << (nhist+1) << '\0'<< std::endl;
// histograms.push_back(new TH1F(TString("h")+tmp,TString("cos thetan")+tmp,50,-1.0,1.0));
std::ostringstream strm;
strm << (nhist+1);
- histograms.push_back( new TH1F(TString("h") + strm.str(),
- TString("cos theta") + strm.str(),
+ histograms.push_back( new TH1F(TString("h") + strm.str().c_str(),
+ TString("cos theta") + strm.str().c_str(),
50,-1.0,1.0) );
}
histograms[nhist++]->Fill(costhetan);
}
if (p->getDaug(1)->getNDaug()==2) {
double costhetan=EvtDecayPlaneNormalAngle(p->getParent()->getP4Lab(),
p->getP4Lab(),
p->getDaug(1)->getDaug(0)->getP4Lab(),
p->getDaug(1)->getDaug(1)->getP4Lab());
if (count==0){
// char* tmp=new char[10];
// std::ostrstream strm(tmp,9);
// strm << (nhist+1) << '\0'<< std::endl;
// histograms.push_back(new TH1F(TString("h")+tmp,TString("cos thetan")+tmp,50,-1.0,1.0));
std::ostringstream strm;
strm << (nhist+1);
- histograms.push_back( new TH1F(TString("h") + strm.str(),
- TString("cos theta") + strm.str(),
+ histograms.push_back( new TH1F(TString("h") + strm.str().c_str(),
+ TString("cos theta") + strm.str().c_str(),
50,-1.0,1.0) );
}
histograms[nhist++]->Fill(costhetan);
}
}
}
p=p->nextIter(root_part);
}while(p!=0);
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runPiPi(int nevent, EvtGen &myGenerator) {
std::ofstream outmix;
outmix.open("pipi.dat");
TFile *file=new TFile("pipi.root", "RECREATE");
TH1F* tB0Hist = new TH1F("h1","dt in B->pipi with B0 tag",50,-5.0,5.0);
TH1F* tB0BHist = new TH1F("h2","dt in B->pipi with B0B tag",50,-5.0,5.0);
TH1F* tB0 = new TH1F("h3","t in B->pipi for B0 tag",25,0.0,5.0);
TH1F* tB0B = new TH1F("h4","t in B->pipi for B0B tag",25,0.0,5.0);
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/PIPI.DEC");
//EvtGen myGenerator(decay_name,pdttable_name,myRandomEngine);
myGenerator.readUDecay(udecay_name);
EvtParticle *bcp,*btag;
int count=1;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
static EvtId GAMM=EvtPDL::getId(std::string("gamma"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
if (root_part->getDaug(0)->getNDaug()==3){
btag=root_part->getDaug(0);
bcp=root_part->getDaug(1);
}
else{
bcp=root_part->getDaug(0);
btag=root_part->getDaug(1);
}
EvtId tag; //cp tag
if (btag->getId()==B0B){
tag=B0;
}
else{tag=B0B;
}
// int a1=bcp->getDaug(0)->getId();
if (tag==B0) tB0Hist->Fill( bcp->getLifetime() - btag->getLifetime() );
if (tag==B0) tB0->Fill( btag->getLifetime() );
if (tag==B0B) tB0BHist->Fill( bcp->getLifetime() - btag->getLifetime() );
if (tag==B0B) tB0B->Fill( btag->getLifetime() );
outmix << bcp->getLifetime() << " " <<
btag->getLifetime() << " " << tag.getId() << std::endl;
root_part->deleteTree();
}while (count++<nevent);
outmix.close();
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runA1Pi(int nevent, EvtGen &myGenerator) {
std::ofstream outmix;
outmix.open("a1pi.dat");
int count=1;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/A1PI.DEC");
myGenerator.readUDecay(udecay_name);
EvtParticle *bcp,*btag;
EvtParticle *a1,*rho0,*pi1,*pi2,*pi3,*pi4;
EvtVector4R p4bcp,p4a1,p4rho0,p4pi1,p4pi2,p4pi3,p4pi4;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
if (root_part->getDaug(0)->getNDaug()==3){
btag=root_part->getDaug(0);
bcp=root_part->getDaug(1);
}
else{
bcp=root_part->getDaug(0);
btag=root_part->getDaug(1);
}
a1=bcp->getDaug(0);
pi1=bcp->getDaug(1);
rho0=a1->getDaug(0);
pi2=a1->getDaug(1);
pi3=rho0->getDaug(0);
pi4=rho0->getDaug(1);
p4bcp=bcp->getP4Lab();
p4a1=a1->getP4Lab();
p4pi1=pi1->getP4Lab();
p4rho0=rho0->getP4Lab();
p4pi2=pi2->getP4Lab();
p4pi3=pi3->getP4Lab();
p4pi4=pi4->getP4Lab();
EvtId tag; //cp tag
if (btag->getId()==B0B){
tag=B0;
}
else{
tag=B0B;
}
outmix << bcp->getLifetime() << " " <<
btag->getLifetime() << " " <<
EvtDecayAngle(p4bcp,p4rho0+p4pi2,p4rho0) << " "<<
EvtDecayAngle(p4a1,p4pi3+p4pi4,p4pi3) << " "<<
EvtPDL::getStdHep(tag) << std::endl;
root_part->deleteTree();
}while (count++<1000);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runCPTest(int nevent, EvtGen &myGenerator) {
std::ofstream outmix;
outmix.open("cptest.dat");
int count=1;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/CPTEST.DEC");
//EvtGen myGenerator(decay_name,pdttable_name,myRandomEngine);
myGenerator.readUDecay(udecay_name);
EvtParticle *bcp,*btag;
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle* root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
if (root_part->getDaug(0)->getNDaug()==3){
btag=root_part->getDaug(0);
bcp=root_part->getDaug(1);
}
else{
bcp=root_part->getDaug(0);
btag=root_part->getDaug(1);
}
EvtId tag; //cp tag
if (btag->getId()==B0B){
tag=B0;
}
else{
tag=B0B;
}
outmix << bcp->getLifetime() << " " <<
btag->getLifetime() << " " << tag.getId() << std::endl;
root_part->deleteTree();
}while (count++<nevent);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runBtoXsgamma(int nevent, EvtGen &myGenerator) {
static EvtId UPS4 = EvtPDL::getId(std::string("Upsilon(4S)"));
TFile *file=new TFile("BtoXsgamma.root", "RECREATE");
int count=1;
EvtParticle* root_part;
EvtVectorParticle *vector_part;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/BTOXSGAMMA.DEC");
myGenerator.readUDecay(udecay_name);
// Plot kinematics for b->s,gamma
int strangeid,antistrangeid;
int Bmulti,bId1a,bId1b,bId2a,bId2b,b1Id,b2Id;
do{
vector_part = new EvtVectorParticle;
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
vector_part->init(UPS4,p_init);
root_part = (EvtParticle *)vector_part;
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *B1 = root_part->getDaug(0);
Bmulti= B1->getNDaug();
if (Bmulti==1) B1 = B1->getDaug(0);
EvtId BId1a = B1->getDaug(0)->getId();
bId1a = EvtPDL::getStdHep(BId1a);
EvtId BId1b = B1->getDaug(1)->getId();
bId1b = EvtPDL::getStdHep(BId1b);
if (Bmulti==1) EvtGenReport(EVTGEN_INFO,"EvtGen") << "B1" << " bId1a=" << bId1a << " bId1b=" << bId1b << " ndaug=" << B1->getNDaug() << " Bid=" << EvtPDL::getStdHep(B1->getId()) << std::endl;
EvtParticle *B2 = root_part->getDaug(1);
Bmulti= B2->getNDaug();
if (Bmulti==1) B2 = B2->getDaug(0); // B has a daughter which is a string
EvtId BId2a = B2->getDaug(0)->getId();
bId2a = EvtPDL::getStdHep(BId2a);
EvtId BId2b = B2->getDaug(1)->getId();
bId2b = EvtPDL::getStdHep(BId2b);
if (Bmulti==1) EvtGenReport(EVTGEN_INFO,"EvtGen") << "B2" << " bId2a=" << bId2a << " bId2b=" << bId2b << " ndaug=" << B2->getNDaug() << " Bid=" << EvtPDL::getStdHep(B2->getId()) << std::endl;
EvtId B1Id = B1->getId();
b1Id = EvtPDL::getStdHep(B1Id);
EvtId B2Id = B2->getId();
b2Id = EvtPDL::getStdHep(B2Id);
strangeid=0;
antistrangeid=0;
if ((b1Id == 511) || (b1Id == -511) || (b2Id == 511) || (b2Id == -511)) {
strangeid=30343; antistrangeid=-30343;
} else if ((b1Id == 521) || (b1Id == -521) || (b2Id == 521) || (b2Id == -521)) {
strangeid=30353; antistrangeid=-30353;
} else if ((b1Id == 531) || (b1Id == -531) || (b2Id == 531) || (b2Id == -531)) {
strangeid=30363; antistrangeid=-30363;
}
EvtGenReport(EVTGEN_INFO,"EvtGen") << "bId1a "<<bId1a<<" bId1b "<<bId1b<<" bId2a "<<bId2a<<" bId2b "<<bId2b<<" for event "<<count<<std::endl;
EvtParticle *Bpeng = 0;
int bnum=0;
int pengcount=0;
if (((bId1a == strangeid) && (bId1b == 22)) || ((bId1a == antistrangeid) && (bId1b == 22))|| ((bId1b == strangeid) && (bId1a == 22)) || ((bId1b == antistrangeid) && (bId1a == 22))) {
Bpeng = B1;
bnum=1;
pengcount++;
}
if (((bId2a == strangeid) && (bId2b == 22)) || ((bId2a == antistrangeid) && (bId2b == 22)) || ((bId2b == strangeid) && (bId2a == 22)) || ((bId2b == antistrangeid) && (bId2a == 22))) {
Bpeng = B2;
bnum=2;
pengcount++;
}
if (pengcount == 0) {
Bpeng=B1;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "No penguin decay for event "<<count<<std::endl;
bnum=0;
} else if (pengcount == 2) {
Bpeng=B1;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Two penguin decays in event "<<count<<std::endl;
bnum=0;
}
Bmulti = Bpeng->getNDaug();
EvtParticle *Xs = Bpeng->getDaug(0);
EvtParticle *gam = Bpeng->getDaug(1);
EvtVector4R p4Xs = Xs->getP4Lab();
EvtId BId = Bpeng->getId();
EvtId XsId = Xs->getId();
int Xsmulti = Xs->getNDaug();
EvtId gamId = gam->getId();
//int bId = EvtPDL::getStdHep(BId);
//int XId = EvtPDL::getStdHep(XsId);
//int gId = EvtPDL::getStdHep(gamId);
//float XsMass = p4Xs.mass();
//double gmass = p4gam.mass();
//double genergy = p4gam.get(0);
// debug stuff: EvtGenReport(EVTGEN_INFO,"EvtGen") << "bnum=" << bnum << " pengcount=" << pengcount << " bId=" << bId << " Bmulti=" << Bmulti << " XsId=" << XId << " gId=" << gId << std::endl;
//need to change this to root...I don't have the energy now
//tuple->column("bnum", bnum);
//tuple->column("pengcount", pengcount);
//tuple->column("bId", bId);
//tuple->column("Bmulti", Bmulti);
//tuple->column("XsId", XId);
//tuple->column("gId", gId);
//tuple->column("XsMass", XsMass);
//tuple->column("Xsmulti", Xsmulti, 0,"Xs", HTRange<int>(0,200));
//tuple->column("gmass", gmass);
//tuple->column("genergy", genergy);
//HTValOrderedVector<int> XDaugId, XDaugNephewId;
//HTValOrderedVector<float> XsDaugMass, XsDaugNephewMass;
int nTot(0);
for(int i=0;i<Xsmulti;i++){
EvtParticle *XsDaug = Xs->getDaug(i);
EvtVector4R p4XsDaug = XsDaug->getP4Lab();
EvtId XsDaugId = XsDaug->getId();
//XDaugId.push_back(EvtPDL::getStdHep(XsDaugId));
//XsDaugMass.push_back( p4XsDaug.mass());
int Daumulti = XsDaug->getNDaug();
if(abs(EvtPDL::getStdHep(XsDaugId))==321||EvtPDL::getStdHep(XsDaugId)==310||EvtPDL::getStdHep(XsDaugId)==111||abs(EvtPDL::getStdHep(XsDaugId))==211||Daumulti==0){
nTot++;
EvtVector4R p4XsDaugNephew = XsDaug->getP4Lab();
EvtId XsDaugNephewId =XsDaug->getId() ;
//XDaugNephewId.push_back(EvtPDL::getStdHep(XsDaugId));
//XsDaugNephewMass.push_back( p4XsDaug.mass());
}else if(Daumulti!=0){
for(int k=0;k<Daumulti;k++){
EvtParticle *XsDaugNephew = XsDaug->getDaug(k);
EvtId XsDaugNephewId = XsDaugNephew->getId();
int Nephmulti = XsDaugNephew->getNDaug();
if(Nephmulti==0||abs(EvtPDL::getStdHep(XsDaugNephewId))==321||EvtPDL::getStdHep(XsDaugNephewId)==310||EvtPDL::getStdHep(XsDaugNephewId)==111||abs(EvtPDL::getStdHep(XsDaugNephewId))==211) {
nTot++;
EvtVector4R p4XsDaugNephew = XsDaugNephew->getP4Lab();
//XDaugNephewId.push_back(EvtPDL::getStdHep(XsDaugNephewId));
//XsDaugNephewMass.push_back( p4XsDaugNephew.mass());
}else{
for(int g=0;g<Nephmulti;g++){
nTot++;
EvtParticle *XsDaugNephewNephew = XsDaugNephew->getDaug(g);
EvtVector4R p4XsDaugNephewNephew = XsDaugNephewNephew->getP4Lab();
EvtId XsDaugNephewNephewId = XsDaugNephewNephew->getId();
//XDaugNephewId.push_back(EvtPDL::getStdHep(XsDaugNephewNephewId));
//XsDaugNephewMass.push_back( p4XsDaugNephewNephew.mass());
}
}
}
}
}
//tuple->column("XsDaugId", XDaugId,"Xsmulti", 0, "Xs");
//tuple->column("XsDaugMass", XsDaugMass,"Xsmulti", 0, "Xs");
//tuple->column("nTot", nTot, 0,"nTot", HTRange<int>(0,200));
//tuple->column("XsDaugNephewId", XDaugNephewId,"nTot", 0, "nTot");
//tuple->column("XsDaugNephewMass", XsDaugNephewMass,"nTot", 0, "nTot");
//tuple->dumpData();
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen")<<"End EvtGen. Ran on "<<nevent<<" events."<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runBtoK1273gamma(int nevent, EvtGen &myGenerator) {
static EvtId UPS4 = EvtPDL::getId(std::string("Upsilon(4S)"));
TFile *file=new TFile("BtoK1273gamma.root", "RECREATE");
//HepTuple *tuple = hfile.ntuple("BtoK1273gamma", 1);
int count=1;
EvtParticle *root_part;
EvtVectorParticle *vector_part;
char udecay_name[100];
strcpy(udecay_name,"exampleFiles/BTOK1273GAMMA.DEC");
myGenerator.readUDecay(udecay_name);
// Plot kinematics for b->s,gamma
int strangeid,antistrangeid;
int Bmulti,bId1a,bId1b,bId2a,bId2b,b1Id,b2Id;
do{
vector_part = new EvtVectorParticle;
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
vector_part->init(UPS4,p_init);
root_part = (EvtParticle *)vector_part;
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *B1 = root_part->getDaug(0);
Bmulti= B1->getNDaug();
if (Bmulti==1) B1 = B1->getDaug(0);
EvtId BId1a = B1->getDaug(0)->getId();
bId1a = EvtPDL::getStdHep(BId1a);
EvtId BId1b = B1->getDaug(1)->getId();
bId1b = EvtPDL::getStdHep(BId1b);
if (Bmulti==1) EvtGenReport(EVTGEN_INFO,"EvtGen") << "B1" << " bId1a=" << bId1a << " bId1b=" << bId1b << " ndaug=" << B1->getNDaug() << " Bid=" << EvtPDL::getStdHep(B1->getId()) << std::endl;
EvtParticle *B2 = root_part->getDaug(1);
Bmulti= B2->getNDaug();
if (Bmulti==1) B2 = B2->getDaug(0); // B has a daughter which is a string
EvtId BId2a = B2->getDaug(0)->getId();
bId2a = EvtPDL::getStdHep(BId2a);
EvtId BId2b = B2->getDaug(1)->getId();
bId2b = EvtPDL::getStdHep(BId2b);
if (Bmulti==1) EvtGenReport(EVTGEN_INFO,"EvtGen") << "B2" << " bId2a=" << bId2a << " bId2b=" << bId2b << " ndaug=" << B2->getNDaug() << " Bid=" << EvtPDL::getStdHep(B2->getId()) << std::endl;
EvtId B1Id = B1->getId();
b1Id = EvtPDL::getStdHep(B1Id);
EvtId B2Id = B2->getId();
b2Id = EvtPDL::getStdHep(B2Id);
strangeid=0;
antistrangeid=0;
if ((b1Id == 511) || (b1Id == -511) || (b2Id == 511) || (b2Id == -511)) {
strangeid=10313; antistrangeid=-10313;
} else if ((b1Id == 521) || (b1Id == -521) || (b2Id == 521) || (b2Id == -521)) {
strangeid=10323; antistrangeid=-10323;
}
EvtGenReport(EVTGEN_INFO,"EvtGen") << "bId1a "<<bId1a<<" bId1b "<<bId1b<<" bId2a "<<bId2a<<" bId2b "<<bId2b<<" for event "<<count<<std::endl;
EvtParticle *Bpeng = 0;
int bnum=0;
int pengcount=0;
if (((bId1a == strangeid) && (bId1b == 22)) || ((bId1a == antistrangeid) && (bId1b == 22))|| ((bId1b == strangeid) && (bId1a == 22)) || ((bId1b == antistrangeid) && (bId1a == 22))) {
Bpeng = B1;
bnum=1;
pengcount++;
}
if (((bId2a == strangeid) && (bId2b == 22)) || ((bId2a == antistrangeid) && (bId2b == 22)) || ((bId2b == strangeid) && (bId2a == 22)) || ((bId2b == antistrangeid) && (bId2a == 22))) {
Bpeng = B2;
bnum=2;
pengcount++;
}
if (pengcount == 0) {
Bpeng=B1;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "No penguin decay for event "<<count<<std::endl;
bnum=0;
} else if (pengcount == 2) {
Bpeng=B1;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Two penguin decays in event "<<count<<std::endl;
bnum=0;
}
Bmulti = Bpeng->getNDaug();
EvtParticle *Ks = Bpeng->getDaug(0);
//EvtParticle *gam = Bpeng->getDaug(1);
EvtVector4R p4Ks = Ks->getP4Lab();
//const EvtVector4R& p4gam = gam->getP4(); // gamma 4-mom in parent's rest frame
//EvtId BId = Bpeng->getId();
//EvtId KsId = Ks->getId();
int Ksmulti = Ks->getNDaug();
//EvtId gamId = gam->getId();
//int bId = EvtPDL::getStdHep(BId);
//int XId = EvtPDL::getStdHep(KsId);
//int gId = EvtPDL::getStdHep(gamId);
//double KsMass = p4Ks.mass();
//double gmass = p4gam.mass();
//double genergy = p4gam.get(0);
// debug stuff: EvtGenReport(EVTGEN_INFO,"EvtGen") << "bnum=" << bnum << " pengcount=" << pengcount << " bId=" << bId << " Bmulti=" << Bmulti << " KsId=" << XId << " gId=" << gId << std::endl;
//tuple->column("bnum", bnum);
//tuple->column("pengcount", pengcount);
//tuple->column("bId", bId);
//tuple->column("Bmulti", Bmulti);
//tuple->column("KsId", XId);
//tuple->column("gId", gId);
//tuple->column("KsMass", KsMass);
//tuple->column("Ksmulti", Ksmulti);
//tuple->column("gmass", gmass);
//tuple->column("genergy", genergy);
for(int i=0;i<Ksmulti;i++){
EvtParticle *KsDaug = Ks->getDaug(i);
EvtVector4R p4KsDaug = KsDaug->getP4Lab();
EvtId KsDaugId = KsDaug->getId();
//int XDaugId = EvtPDL::getStdHep(KsDaugId);
//double KsDaugMass = p4KsDaug.mass();
//tuple->column("KsDaugId", XDaugId);
//tuple->column("KsDaugMass", KsDaugMass);
}
//tuple->dumpData();
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen")<<"End EvtGen. Ran on "<<nevent<<" events."<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runCheckRotBoost(){
EvtDiracSpinor sp1,sp2;
//Generate ryd/lange random spinors.
sp1.set(EvtComplex(1.0,-2.0),
EvtComplex(3.0,1.0),
EvtComplex(-4.5,0.5),
EvtComplex(0.2,-0.5));
sp2.set(EvtComplex(0.1,-1.0),
EvtComplex(1.2,-0.5),
EvtComplex(3.6,1.8),
EvtComplex(-0.2,-0.6));
EvtComplex s=EvtLeptonSCurrent(sp1,sp2);
EvtComplex p=EvtLeptonPCurrent(sp1,sp2);
EvtVector4C a=EvtLeptonACurrent(sp1,sp2);
EvtVector4C v=EvtLeptonVCurrent(sp1,sp2);
EvtVector4C va=EvtLeptonVACurrent(sp1,sp2);
EvtTensor4C t=EvtLeptonTCurrent(sp1,sp2);
//start with boosts...
EvtVector4R ranBoost(2.0,0.4,-0.8,0.3);
EvtDiracSpinor sp1Boost=boostTo(sp1,ranBoost);
EvtDiracSpinor sp2Boost=boostTo(sp2,ranBoost);
EvtComplex sBoost=EvtLeptonSCurrent(sp1Boost,sp2Boost);
EvtComplex pBoost=EvtLeptonPCurrent(sp1Boost,sp2Boost);
EvtVector4C aBoost=EvtLeptonACurrent(sp1Boost,sp2Boost);
EvtVector4C vBoost=EvtLeptonVCurrent(sp1Boost,sp2Boost);
EvtVector4C vaBoost=EvtLeptonVACurrent(sp1Boost,sp2Boost);
EvtTensor4C tBoost=EvtLeptonTCurrent(sp1Boost,sp2Boost);
EvtVector4C aDirBoost=boostTo(a,ranBoost);
EvtVector4C vDirBoost=boostTo(v,ranBoost);
EvtVector4C vaDirBoost=boostTo(va,ranBoost);
EvtTensor4C tDirBoost(t);
tDirBoost.applyBoostTo(ranBoost);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Comparing after doing a random boost"<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Scalar "<<s<<" "<<sBoost<<s-sBoost<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "PseudoScalar "<<p<<" "<<pBoost<<p-pBoost<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "AxialVector "<<aDirBoost<<" "<<aBoost<<aDirBoost-aBoost<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Vector "<<vDirBoost<<" "<<vBoost<<vDirBoost-vBoost<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "V-A "<<vaDirBoost<<" "<<vaBoost<<vaDirBoost-vaBoost<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Tensor "<<tDirBoost<<" "<<tBoost<<tDirBoost-tBoost<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Done comparing after doing a random boost"<<std::endl;
//Now do rotations...
//start with boosts...
double alpha=0.4;
double beta=-0.61;
double gamma=3.0;
EvtDiracSpinor sp1Rot=rotateEuler(sp1,alpha,beta,gamma);
EvtDiracSpinor sp2Rot=rotateEuler(sp2,alpha,beta,gamma);
EvtComplex sRot=EvtLeptonSCurrent(sp1Rot,sp2Rot);
EvtComplex pRot=EvtLeptonPCurrent(sp1Rot,sp2Rot);
EvtVector4C aRot=EvtLeptonACurrent(sp1Rot,sp2Rot);
EvtVector4C vRot=EvtLeptonVCurrent(sp1Rot,sp2Rot);
EvtVector4C vaRot=EvtLeptonVACurrent(sp1Rot,sp2Rot);
EvtTensor4C tRot=EvtLeptonTCurrent(sp1Rot,sp2Rot);
EvtVector4C aDirRot(a);
EvtVector4C vDirRot(v);
EvtVector4C vaDirRot(va);
EvtTensor4C tDirRot(t);
aDirRot.applyRotateEuler(alpha,beta,gamma);
vDirRot.applyRotateEuler(alpha,beta,gamma);
vaDirRot.applyRotateEuler(alpha,beta,gamma);
tDirRot.applyRotateEuler(alpha,beta,gamma);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Comparing after doing a random rotation"<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Scalar "<<s<<" "<<sRot<<s-sRot<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "PseudoScalar "<<p<<" "<<pRot<<p-pRot<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "AxialVector "<<aDirRot<<" "<<aRot<<aDirRot-aRot<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Vector "<<vDirRot<<" "<<vRot<<vDirRot-vRot<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "V-A "<<vaDirRot<<" "<<vaRot<<vaDirRot-vaRot<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Tensor "<<tDirRot<<" "<<tRot<<tDirRot-tRot<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "Done comparing after doing a random rotation"<<std::endl;
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
int countInclusive(std::string name, EvtParticle *root_part,TH1F* mom, TH1F* mass){
EvtParticle *p = root_part;
int temp=0;
EvtId searchFor=EvtPDL::getId(name);
do{
EvtId type=p->getId();
if (type==searchFor) {
temp+=1;
if ( mom) mom->Fill(p->getP4Lab().d3mag());
if ( mass) mass->Fill(p->mass());
//if ( theBs.contains(p->getParent()->getId()) ) {
//dirPsimom->Fill(p->getP4Lab().d3mag());
//}
//EvtGenReport(EVTGEN_INFO,"EvtGen") << "LANGE " << p->getP4Lab().d3mag() << " " << p->getP4Lab().get(3)/p->getP4Lab().d3mag() << std::endl;
}
p=p->nextIter(root_part);
}while(p!=0);
return temp;
}
int countInclusiveSubTree(std::string name, EvtParticle *root_part,
EvtIdSet setIds,
TH1F* mom) {
int temp=0;
EvtParticle *p=root_part;
do{
if ( setIds.contains(p->getId()) ) {
temp+=countInclusive(name,p);
}
//p->printTree();
p=p->nextIter(root_part);
}while(p!=0);
//EvtGenReport(EVTGEN_INFO,"EvtGen") << "done"<<std::endl;
return temp;
}
int countInclusiveParent(std::string name, EvtParticle *root_part,
EvtIdSet setIds,
TH1F* mom) {
EvtParticle *p = root_part;
int temp=0;
EvtId searchFor=EvtPDL::getId(name);
do{
EvtId type=p->getId();
if (type==searchFor) {
if ( p->getParent() ) {
if ( setIds.contains(p->getParent()->getId()) ) {
temp+=1;
if ( mom) mom->Fill(p->getP4Lab().d3mag());
}
}
}
p=p->nextIter(root_part);
}while(p!=0);
return temp;
}
void runBMix(int nevent,EvtGen &myGenerator,std::string userFile,std::string rootFile) {
TFile *file=new TFile(rootFile.c_str(), "RECREATE");
TH1F* b0_b0 = new TH1F("h1","dt B0-B0",100,-15.0,15.0);
TH1F* b0b_b0b = new TH1F("h2","dt B0B-B0B",100,-15.0,15.0);
TH1F* b0b_b0 = new TH1F("h3","dt B0B-B0",100,-15.0,15.0);
TH1F* b0_b0b = new TH1F("h4","dt B0-B0B",100,-15.0,15.0);
int count=1;
myGenerator.readUDecay(userFile.c_str());
EvtId b0=EvtPDL::getId("B0");
EvtId b0b=EvtPDL::getId("anti-B0");
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
std::ofstream outmix;
TString outFileName(rootFile);
outFileName.ReplaceAll(".root", ".dat");
outmix.open(outFileName.Data());
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtId id1=root_part->getDaug(0)->getId();
EvtId id2=root_part->getDaug(1)->getId();
double t1=root_part->getDaug(0)->getLifetime();
double t2=root_part->getDaug(1)->getLifetime();
double dt=(t1-t2)/(1e-12*3e11);
if (id1==b0&&id2==b0) { b0_b0->Fill(dt); outmix << dt << " 1"<<std::endl;}
if (id1==b0b&&id2==b0b) {b0b_b0b->Fill(dt); outmix << dt << " 2"<<std::endl;}
if (id1==b0b&&id2==b0) {b0b_b0->Fill(dt); outmix << dt << " 0"<<std::endl;}
if (id1==b0&&id2==b0b) {b0_b0b->Fill(dt); outmix << dt << " 0"<<std::endl;}
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runDDalitz(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("ddalitz.root", "RECREATE");
TH2F* dalitz = new TH2F("h1","m^2!?[K-[p]+! vs m^2!?[K-[p]0!",
70,0.0,3.5,70,0.0,3.5);
TH2F* dalitz2 = new TH2F("h5","m^2!([p]^-![p]^0!) vs m^2!([K-[p]+!",
100,0.0,3.5,100,0.0,2.0);
TH1F* m12=new TH1F("h2","m?[K-[p]+!",100,0.0,3.0);
TH1F* m13=new TH1F("h3","m?[K-[p]0!",100,0.0,3.0);
TH1F* m23=new TH1F("h4","m?[[p]+[p]0!",100,0.0,2.0);
int count;
myGenerator.readUDecay("exampleFiles/DDALITZ.DEC");
count=1;
static EvtId D0=EvtPDL::getId(std::string("D0"));
std::ofstream outmix;
outmix.open("ddalitz.dat");
do{
EvtVector4R p_init(EvtPDL::getMass(D0),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(D0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
EvtVector4R p1=root_part->getDaug(0)->getP4Lab();
EvtVector4R p2=root_part->getDaug(1)->getP4Lab();
EvtVector4R p3=root_part->getDaug(2)->getP4Lab();
dalitz->Fill((p1+p2).mass2(),(p1+p3).mass2(),1.0);
dalitz2->Fill((p1+p2).mass2(),(p2+p3).mass2(),1.0);
m12->Fill((p1+p2).mass2());
m13->Fill((p1+p3).mass2());
m23->Fill((p2+p3).mass2());
outmix << (p1+p2).mass2() << " "<<(p2+p3).mass2()<<" "<<(p1+p3).mass2()<<std::endl;
root_part->deleteTree();
if(count == nevent-1) {
std::ofstream testi("testi.dat");
double val = m12->GetMean();
double errval = m12->GetMeanError();
testi << "evtgenlhc_test1 1 " << val << " " << errval << std::endl;
val = m23->GetMean();
errval = m23->GetMeanError();
testi << "evtgenlhc_test1 2 " << val << " " << errval << std::endl;
}
}while (count++<nevent);
file->Write(); file->Close();
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runPiPiPi(int nevent, EvtGen &myGenerator) {
std::ofstream outmix;
outmix.open("pipipi.dat");
int count;
EvtVector4R p4pip,p4pim,p4pi0;
myGenerator.readUDecay("exampleFiles/PIPIPI.DEC");
count=1;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
p4pip=root_part->getDaug(0)->getDaug(0)->getP4Lab();
p4pim=root_part->getDaug(0)->getDaug(1)->getP4Lab();
p4pi0=root_part->getDaug(0)->getDaug(2)->getP4Lab();
outmix << root_part->getDaug(0)->getLifetime() << " " <<
root_part->getDaug(1)->getLifetime()<< " ";
outmix << (p4pip+p4pim).mass2()<<" "<<
(p4pip+p4pi0).mass2()<<std::endl;
root_part->deleteTree();
}while (count++<10000);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runBHadronic(int nevent, EvtGen &myGenerator) {
std::ofstream outmix;
outmix.open("bhadronic.dat");
int count;
myGenerator.readUDecay("exampleFiles/BHADRONIC.DEC");
static EvtId B0=EvtPDL::getId(std::string("B0"));
count=1;
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(B0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p;
// root_part->printTree();
p=root_part;
do{
outmix << p->getId().getId()<<" "<< p->getP4Lab().d3mag() << std::endl;
p=p->nextIter();
}while(p!=0);
root_part->deleteTree();
}while (count++<nevent);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runSingleB(int nevent, EvtGen &myGenerator) {
int count;
static EvtId B0=EvtPDL::getId(std::string("B0"));
count=1;
do{
EvtVector4R p_init(EvtPDL::getMass(B0),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(B0,
p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
root_part->deleteTree();
}while (count++<nevent);
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runPiPiPiPi(int nevent, EvtGen &myGenerator) {
std::ofstream outmix;
outmix.open("pipipipi.dat");
int count;
EvtVector4R p4pi1,p4pi2,p4pi3,p4pi4;
myGenerator.readUDecay("exampleFiles/PIPIPIPI.DEC");
count=1;
EvtParticle *bcp,*btag;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
if (root_part->getDaug(0)->getNDaug()==3){
btag=root_part->getDaug(0);
bcp=root_part->getDaug(1);
}
else{
bcp=root_part->getDaug(0);
btag=root_part->getDaug(1);
}
EvtId tag; //cp tag
if (btag->getId()==B0B){
tag=B0;
}
else{
tag=B0B;
}
p4pi1=bcp->getDaug(0)->getP4Lab();
p4pi2=bcp->getDaug(1)->getP4Lab();
p4pi3=bcp->getDaug(2)->getP4Lab();
p4pi4=bcp->getDaug(3)->getP4Lab();
EvtVector4R p4bcp,p4rho0,p4a2;
p4rho0=p4pi1+p4pi2;
p4a2=p4rho0+p4pi3;
p4bcp=p4a2+p4pi4;
outmix << bcp->getLifetime() << " " <<
btag->getLifetime() << " " << tag.getId() <<" "
<< (p4pi1+p4pi2+p4pi3).mass() << " " <<
(p4pi1+p4pi2).mass() << " " <<
EvtDecayAngle(p4bcp,p4rho0+p4pi3,p4rho0) << " "<<
EvtDecayAngle(p4a2,p4pi1+p4pi2,p4pi1) << std::endl;
root_part->deleteTree();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "count:"<<count<<std::endl;
}while (count++<1000);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runA2Pi(int nevent, EvtGen &myGenerator) {
std::ofstream outmix;
outmix.open("a2pi.dat");
int count;
myGenerator.readUDecay("exampleFiles/A2PI.DEC");
EvtParticle *bcp,*btag;
EvtParticle *a2,*rho0,*pi1,*pi2,*pi3,*pi4;
EvtVector4R p4bcp,p4a2,p4rho0,p4pi1,p4pi2,p4pi3,p4pi4;
count=1;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
if (root_part->getDaug(0)->getNDaug()==3){
btag=root_part->getDaug(0);
bcp=root_part->getDaug(1);
}
else{
bcp=root_part->getDaug(0);
btag=root_part->getDaug(1);
}
a2=bcp->getDaug(0);
pi1=bcp->getDaug(1);
rho0=a2->getDaug(0);
pi2=a2->getDaug(1);
pi3=rho0->getDaug(0);
pi4=rho0->getDaug(1);
p4bcp=bcp->getP4Lab();
p4a2=a2->getP4Lab();
p4pi1=pi1->getP4Lab();
p4rho0=rho0->getP4Lab();
p4pi2=pi2->getP4Lab();
p4pi3=pi3->getP4Lab();
p4pi4=pi4->getP4Lab();
EvtId tag; //cp tag
if (btag->getId()==B0B){
tag=B0;
}
else{
tag=B0B;
}
outmix << bcp->getLifetime() << " " <<
btag->getLifetime() << " " <<
EvtDecayAngle(p4bcp,p4rho0+p4pi2,p4rho0) << " "<<
EvtDecayAngle(p4a2,p4pi3+p4pi4,p4pi3) << " "<<
EvtPDL::getStdHep(tag) << std::endl;
root_part->deleteTree();
}while (count++<1000);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runHelAmp(int nevent, EvtGen &myGenerator, std::string userFile,
std::string rootFile) {
TFile *file=new TFile(rootFile.c_str(), "RECREATE");
TH1F* costheta = new TH1F("h1","costheta",100,-1.0,1.0);
TH1F* costheta2 = new TH1F("h2","costheta2",100,-1.0,1.0);
int count;
myGenerator.readUDecay(userFile.c_str());
count=1;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtVector4R d14=root_part->getDaug(0)->getP4Lab();
double c=d14.get(3)/d14.d3mag();
costheta->Fill(c);
EvtVector4R p=root_part->getP4Lab();
EvtVector4R q=root_part->getDaug(0)->getP4Lab();
EvtVector4R d=root_part->getDaug(0)->getDaug(0)->getP4Lab();
costheta2->Fill(EvtDecayAngle(p,q,d));
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runHelAmp2(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("helamp2.root", "RECREATE");
TH1F* costheta = new TH1F("h1","costheta",100,-1.0,1.0);
int count;
myGenerator.readUDecay("exampleFiles/HELAMP2.DEC");
count=1;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtVector4R p=root_part->getDaug(0)->getDaug(0)->getP4Lab();
EvtVector4R q=root_part->getDaug(0)->getDaug(0)->getDaug(0)->getP4Lab();
EvtVector4R d=root_part->getDaug(0)->getDaug(0)->getDaug(0)->getDaug(0)->getP4Lab();
costheta->Fill(EvtDecayAngle(p,q,d));
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runD2Pi(int nevent, EvtGen &myGenerator) {
TFile *file=new TFile("d2pi.root", "RECREATE");
TH1F* cospi = new TH1F("h1","cos[Q]?[p]!",50,-1.0,1.0);
TH1F* ptpi = new TH1F("h2","Pt of pion",50,0.0,1.5);
TH1F* ppi = new TH1F("h3","P?[p]! in [Y](4S) rest frame",50,0.0,1.5);
int count;
myGenerator.readUDecay("exampleFiles/D2PI.DEC");
EvtParticle *b;
EvtParticle *d2,*pi;
EvtVector4R p4b,p4d2,p4pi;
count=1;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
b=root_part->getDaug(0);
d2=b->getDaug(0);
pi=d2->getDaug(1);
p4b=b->getP4Lab();
p4d2=d2->getP4Lab();
p4pi=pi->getP4Lab();
cospi->Fill(EvtDecayAngle(p4b,p4d2,p4pi));
ptpi->Fill(sqrt(p4pi.get(2)*p4pi.get(2)+p4pi.get(3)*p4pi.get(3)));
ppi->Fill(p4pi.d3mag());
root_part->deleteTree();
}while (count++<nevent);
file->Write(); file->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runPiPiCPT(int nevent, EvtGen &myGenerator) {
std::ofstream outmix;
outmix.open("pipicpt.dat");
int count;
myGenerator.readUDecay("exampleFiles/PIPICPT.DEC");
EvtParticle *bcp,*btag;
count=1;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
if (root_part->getDaug(0)->getNDaug()==3){
btag=root_part->getDaug(0);
bcp=root_part->getDaug(1);
}
else{
bcp=root_part->getDaug(0);
btag=root_part->getDaug(1);
}
EvtId tag; //cp tag
if (btag->getId()==B0B){
tag=B0;
}
else{
tag=B0B;
}
outmix << bcp->getLifetime() << " " <<
btag->getLifetime() << " " << tag.getId() << std::endl;
root_part->deleteTree();
}while (count++<10000);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runJpsiKs(int nevent, EvtGen &myGenerator) {
std::ofstream outmix;
outmix.open("jpsiks.dat");
int count;
myGenerator.readUDecay("exampleFiles/JPSIKS.DEC");
EvtParticle *bcp,*btag;
count=1;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId B0=EvtPDL::getId(std::string("B0"));
static EvtId B0B=EvtPDL::getId(std::string("anti-B0"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
if (root_part->getDaug(0)->getNDaug()==3){
btag=root_part->getDaug(0);
bcp=root_part->getDaug(1);
}
else{
bcp=root_part->getDaug(0);
btag=root_part->getDaug(1);
}
EvtId tag; //cp tag
if (btag->getId()==B0B){
tag=B0;
}
else{
tag=B0B;
}
outmix << bcp->getLifetime() << " " <<
btag->getLifetime() << " " << tag.getId() << std::endl;
root_part->deleteTree();
}while (count++<10000);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runDump(int nevent, EvtGen &myGenerator) {
int count;
int stable_list[1];
stable_list[0]=-1;
std::ofstream outmix;
outmix.open("dump.dat");
EvtParticle *p;
myGenerator.readUDecay("exampleFiles/DUMP.DEC");
count=1;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId PI0=EvtPDL::getId(std::string("pi0"));
static EvtId PIP=EvtPDL::getId(std::string("pi+"));
static EvtId PIM=EvtPDL::getId(std::string("pi-"));
static EvtId EP=EvtPDL::getId(std::string("e+"));
static EvtId KP=EvtPDL::getId(std::string("K+"));
static EvtId MUP=EvtPDL::getId(std::string("mu+"));
static EvtId EM=EvtPDL::getId(std::string("e-"));
static EvtId KM=EvtPDL::getId(std::string("K-"));
static EvtId MUM=EvtPDL::getId(std::string("mu-"));
static EvtId GAMM=EvtPDL::getId(std::string("gamma"));
static EvtId K0L=EvtPDL::getId(std::string("K_L0"));
do{
if (count==100*(count/100)) {
EvtGenReport(EVTGEN_INFO,"EvtGen") << count << std::endl;
}
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
p=root_part;
EvtParticle *otherb=root_part->getDaug(1);
EvtGenReport(EVTGEN_INFO,"EvtGen") <<"Event:"<<count<<std::endl;
root_part->printTree();
outmix << "B"
<<" "<<otherb->getP4Lab().get(0)
<<" "<<otherb->getP4Lab().get(1)
<<" "<<otherb->getP4Lab().get(2)
<<" "<<otherb->getP4Lab().get(3)<<std::endl;
do{
if (p->getId()==PIP||
p->getId()==EP||
p->getId()==KP||
p->getId()==MUP){
outmix << "chg +1"
<<" "<<p->getP4Lab().get(1)
<<" "<<p->getP4Lab().get(2)
<<" "<<p->getP4Lab().get(3)<<std::endl;
}
if (p->getId()==PIM||
p->getId()==EM||
p->getId()==KM||
p->getId()==MUM){
outmix << "chg -1"
<<" "<<p->getP4Lab().get(1)
<<" "<<p->getP4Lab().get(2)
<<" "<<p->getP4Lab().get(3)<<std::endl;
}
if (p->getId()==GAMM||
p->getId()==K0L){
outmix << "neu"
<<" "<<p->getP4Lab().get(1)
<<" "<<p->getP4Lab().get(2)
<<" "<<p->getP4Lab().get(3)<<std::endl;
}
p=p->nextIter();
}while(p!=0);
outmix << "event"<<std::endl;
root_part->deleteTree();
}while (count++<nevent);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runGenericCont(int nevent, EvtGen &myGenerator) {
int count;
std::ofstream outmix;
outmix.open("genericcont.dat");
EvtParticle *p;
myGenerator.readUDecay("exampleFiles/GENERIC.DEC");
count=1;
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
static EvtId VPHO=EvtPDL::getId(std::string("vpho"));
do{
if (count==1000*(count/1000)) {
EvtGenReport(EVTGEN_DEBUG,"EvtGen") << count << std::endl;
}
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(VPHO,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
p=root_part;
do{
outmix << p->getId().getId()<<" "<< p->getP4Lab().d3mag() << std::endl;
p=p->nextIter();
}while(p!=0);
//root_part->printTree();
root_part->deleteTree();
}while (count++<nevent);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runD1(int nevent, EvtGen &myGenerator) {
std::ofstream outmix;
outmix.open("d1.dat");
int count=1;
myGenerator.readUDecay("exampleFiles/D1.DEC");
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
EvtParticle *p_b,*p_d1,*p_e,*p_nu,*p_pi1,*p_dstar,*p_pi2,*p_d;
EvtVector4R p4_b,p4_d1,p4_e,p4_nu,p4_pi1,p4_dstar,p4_pi2,p4_d;
// code for analyzing B->D1 e nu ; D1 -> D* pi ; D* -> D pi
p_b=root_part->getDaug(1);
p_d1=p_b->getDaug(0);
p_e=p_b->getDaug(1);
p_nu=p_b->getDaug(2);
p_dstar=p_d1->getDaug(0);
p_pi1=p_d1->getDaug(1);
p_d=p_dstar->getDaug(0);
p_pi2=p_dstar->getDaug(1);
p4_b=p_b->getP4Lab();
p4_d1=p_d1->getP4Lab();
p4_e=p_e->getP4Lab();
p4_nu=p_nu->getP4Lab();
p4_dstar=p_dstar->getP4Lab();
p4_pi1=p_pi1->getP4Lab();
p4_pi2=p_pi2->getP4Lab();
p4_d=p_d->getP4Lab();
outmix << p4_e.get(0)<<" ";
outmix << (p4_e+p4_nu)*(p4_e+p4_nu)<<" ";
outmix << EvtDecayAngle(p4_b,p4_e+p4_nu,p4_e) << " " ;
outmix << EvtDecayAngle(p4_b,p4_dstar+p4_pi1,p4_dstar) << " " ;
outmix << EvtDecayAngle(p4_d1,p4_d+p4_pi2,p4_d) << " " ;
outmix << EvtDecayAngleChi(p4_b,p4_e,p4_nu,
p4_dstar, p4_pi1 ) << "\n" ;
root_part->deleteTree();
EvtGenReport(EVTGEN_DEBUG,"EvtGen") << "count:"<<count<<std::endl;
}while (count++<nevent);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runMix(int nevent, EvtGen &myGenerator) {
std::ofstream outmix;
outmix.open("mix.dat");
int count=1;
myGenerator.readUDecay("exampleFiles/MIX.DEC");
static EvtId UPS4=EvtPDL::getId(std::string("Upsilon(4S)"));
do{
EvtVector4R p_init(EvtPDL::getMass(UPS4),0.0,0.0,0.0);
EvtParticle *root_part=EvtParticleFactory::particleFactory(UPS4,
p_init);
root_part->setVectorSpinDensity();
myGenerator.generateDecay(root_part);
outmix << root_part->getDaug(0)->getId().getId() << " ";
outmix << root_part->getDaug(0)->getLifetime() << " ";
outmix << root_part->getDaug(1)->getId().getId() << " ";
outmix << root_part->getDaug(1)->getLifetime() << "\n";
root_part->deleteTree();
}while (count++<10000);
outmix.close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
void runBaryonic(int nEvent, EvtGen& myGenerator)
{
TFile* f = new TFile("baryonic.root", "RECREATE");
TH1D* q2Hist = new TH1D("q2Hist", "q square", 50, 0.0, 25.00);
EvtId BMINUS = EvtPDL::getId("B-");
EvtVector4R p_init(EvtPDL::getMass(BMINUS), 0.0, 0.0, 0.0);
EvtParticle* root_part;
myGenerator.readUDecay("exampleFiles/BARYONIC.DEC");
EvtVector4R l;
EvtVector4R p;
EvtVector4R g;
EvtVector4R sum;
for (int i=0; i<nEvent; ++i)
{
root_part = EvtParticleFactory::particleFactory(BMINUS, p_init);
root_part->setDiagonalSpinDensity();
myGenerator.generateDecay(root_part);
l = root_part->getDaug(0)->getP4Lab(); // lambda momentum
p = root_part->getDaug(1)->getP4Lab(); // pBar momentum
g = root_part->getDaug(2)->getP4Lab(); // gamma momentum
sum = p+g;
q2Hist->Fill( sum.mass2() );
root_part->deleteTree();
}
f->Write();
f->Close();
EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n";
}
diff --git a/validation/genExampleRootFiles.cc b/validation/genExampleRootFiles.cc
index 1bd2f2b..ec905c7 100644
--- a/validation/genExampleRootFiles.cc
+++ b/validation/genExampleRootFiles.cc
@@ -1,261 +1,260 @@
// Program to create ROOT files for EvtGen validation plots.
// This looks at the 1st generation daughters and stores 4-momenta
// info into a ROOT file for further analysis.
// Useful for Pythia, Photos and Tauola decay tests.
#include "EvtGen/EvtGen.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtParticleFactory.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtSimpleRandomEngine.hh"
#include "EvtGenBase/EvtMTRandomEngine.hh"
#include "EvtGenBase/EvtDecayBase.hh"
#ifdef EVTGEN_EXTERNAL
#include "EvtGenExternal/EvtExternalGenList.hh"
#endif
//#include "EvtGenBase/EvtHepMCEvent.hh"
#include "HepMC/GenEvent.h"
#include "TFile.h"
#include "TTree.h"
#include <iostream>
#include <string>
#include <vector>
using std::cout;
using std::endl;
-using std::string;
int main(int argc, char** argv) {
- string decayFileName("../DECAY_2010.DEC");
+ std::string decayFileName("../DECAY_2010.DEC");
if (argc > 1) {decayFileName = argv[1];}
cout<<"Decay file name is "<<decayFileName<<endl;
- string rootFileName("evtgenTest.root");
+ std::string rootFileName("evtgenTest.root");
if (argc > 2) {rootFileName = argv[2];}
cout<<"Root file name is "<<rootFileName<<endl;
- string parentName("Upsilon(4S)");
+ std::string parentName("Upsilon(4S)");
if (argc > 3) {parentName = argv[3];}
cout<<"Parent name is "<<parentName<<endl;
int nEvents(10);
if (argc > 4) {nEvents = atoi(argv[4]);}
bool useXml = false;
if(argc > 5) {useXml = (atoi(argv[5])==1);}
bool useEvtGenRandom = false;
if (argc > 6) {useEvtGenRandom = (atoi(argv[6])==1);}
cout<<"Number of events is "<<nEvents<<endl;
TFile* theFile = new TFile(rootFileName.c_str(), "recreate");
TTree* theTree = new TTree("Data", "Data");
TTree* nDaugTree = new TTree("nDaugTree", "nDaugTree");
TTree* dalitzTree = new TTree("dalitzTree", "dalitzTree");
theTree->SetDirectory(theFile);
nDaugTree->SetDirectory(theFile);
dalitzTree->SetDirectory(theFile);
int event(0), nDaug(0), daugId(0);
double E(0.0), p(0.0), px(0.0), py(0.0), pz(0.0);
double t(0.0), x(0.0), y(0.0), z(0.0);
double mass(0.0), lifetime(0.0);
double inv12(0.0), inv13(0.0), inv23(0.0);
double inv12Sq(0.0), inv13Sq(0.0), inv23Sq(0.0);
theTree->Branch("event", &event, "event/I");
theTree->Branch("nDaug", &nDaug, "nDaug/I");
theTree->Branch("id", &daugId, "id/I");
theTree->Branch("E", &E, "E/D");
theTree->Branch("p", &p, "p/D");
theTree->Branch("px", &px, "px/D");
theTree->Branch("py", &py, "py/D");
theTree->Branch("pz", &pz, "pz/D");
theTree->Branch("t", &t, "t/D");
theTree->Branch("x", &x, "x/D");
theTree->Branch("y", &y, "x/D");
theTree->Branch("z", &z, "x/D");
theTree->Branch("mass", &mass, "mass/D");
theTree->Branch("lifetime", &lifetime, "lifetime/D");
nDaugTree->Branch("event", &event, "event/I");
nDaugTree->Branch("nDaug", &nDaug, "nDaug/I");
dalitzTree->Branch("invMass12", &inv12, "invMass12/D");
dalitzTree->Branch("invMass13", &inv13, "invMass13/D");
dalitzTree->Branch("invMass23", &inv23, "invMass23/D");
dalitzTree->Branch("invMass12Sq", &inv12Sq, "invMass12Sq/D");
dalitzTree->Branch("invMass13Sq", &inv13Sq, "invMass13Sq/D");
dalitzTree->Branch("invMass23Sq", &inv23Sq, "invMass23Sq/D");
EvtParticle* baseParticle(0);
EvtParticle* theParent(0);
// Define the random number generator
EvtRandomEngine* myRandomEngine = 0;
#ifdef EVTGEN_CPP11
// Use the Mersenne-Twister generator (C++11 only)
myRandomEngine = new EvtMTRandomEngine();
#else
myRandomEngine = new EvtSimpleRandomEngine();
#endif
// Initialize the generator - read in the decay table and particle properties.
// For our validation purposes, we just want to read in one decay file and create
// plots from that.
EvtAbsRadCorr* radCorrEngine = 0;
std::list<EvtDecayBase*> extraModels;
#ifdef EVTGEN_EXTERNAL
EvtExternalGenList genList(true, "", "gamma", useEvtGenRandom);
radCorrEngine = genList.getPhotosModel();
extraModels = genList.getListOfModels();
#endif
int mixingType(1);
EvtGen myGenerator(decayFileName.c_str(), "../evt.pdl", myRandomEngine,
radCorrEngine, &extraModels, mixingType, useXml);
// If I wanted a user decay file, I would read it in now, e.g:
// myGenerator.readUDecay(otherFileName);
EvtId theId = EvtPDL::getId(parentName);
if (theId.getId() == -1 && theId.getAlias() == -1) {
cout<<"Error. Could not find valid EvtId for "<<parentName<<endl;
return -1;
}
static EvtId stringId = EvtPDL::getId(std::string("string"));
// Loop to create nEvents
int i;
for (i = 0; i < nEvents; i++) {
if (i%1000 == 0) {cout<<"Event number = "<<i+1<<" out of "<<nEvents<<std::endl;}
// Set up the parent particle
EvtVector4R pInit(EvtPDL::getMass(theId), 0.0, 0.0, 0.0);
baseParticle = EvtParticleFactory::particleFactory(theId, pInit);
if (baseParticle->getSpinStates() == 3) {baseParticle->setVectorSpinDensity();}
// Generate the event
myGenerator.generateDecay(baseParticle);
// Alternative way to generate decays and print out information:
// int PDGId = EvtPDL::getStdHep(theId);
// EvtVector4R origin(0.0, 0.0, 0.0, 0.0);
// EvtHepMCEvent* theEvent = myGenerator.generateDecay(PDGId, pInit, origin);
// HepMC::GenEvent* hepMCEvent = theEvent->getEvent();
// hepMCEvent->print();
// Extract other info from the HepMC event. Then delete it:
// delete theEvent;
// Now get the particle decay information, looping through daughter tracks (1st generation only)
// Find out if the first daughter is a string.
// If so, set this as the new parent particle.
EvtId daugEvtId(-1, -1);
EvtParticle* baseDaughter = baseParticle->getDaug(0);
if (baseDaughter != 0) {daugEvtId = baseDaughter->getId();}
if (daugEvtId == stringId) {
theParent = baseDaughter;
} else {
theParent = baseParticle;
}
nDaug = theParent->getNDaug();
int iDaug(0);
nDaugTree->Fill();
//theParent->printTree();
// Loop over the daughter tracks
for (iDaug = 0; iDaug < nDaug; iDaug++) {
EvtParticle* daug = theParent->getDaug(iDaug);
if (daug != 0) {
EvtVector4R p4Lab = daug->getP4Lab();
EvtVector4R pos4 = daug->get4Pos();
// PDG id
daugId = EvtPDL::getStdHep(daug->getId());
// 4-momenta
E = p4Lab.get(0);
px = p4Lab.get(1);
py = p4Lab.get(2);
pz = p4Lab.get(3);
p = sqrt(px*px + py*py + pz*pz);
// 4-position
t = pos4.get(0);
x = pos4.get(1);
y = pos4.get(2);
z = pos4.get(3);
mass = daug->mass();
lifetime = daug->getLifetime();
theTree->Fill();
} // Daughter exists
} // Number of daughters
if (nDaug == 3) {
EvtVector4R p4_d1 = theParent->getDaug(0)->getP4Lab();
EvtVector4R p4_d2 = theParent->getDaug(1)->getP4Lab();
EvtVector4R p4_d3 = theParent->getDaug(2)->getP4Lab();
inv12Sq = (p4_d1+p4_d2)*(p4_d1+p4_d2);
inv13Sq = (p4_d1+p4_d3)*(p4_d1+p4_d3);
inv23Sq = (p4_d2+p4_d3)*(p4_d2+p4_d3);
inv12 = sqrt(inv12Sq);
inv13 = sqrt(inv13Sq);
inv23 = sqrt(inv23Sq);
dalitzTree->Fill();
}
// Cleanup
baseParticle->deleteTree();
}
// Write out the TTree information to the ROOT file
theFile->cd();
theTree->Write();
nDaugTree->Write();
dalitzTree->Write();
theFile->Close();
cout<<"Done."<<endl;
delete myRandomEngine;
return 0;
}