diff --git a/test/Makefile b/test/Makefile index 3843d25..8e150cf 100644 --- a/test/Makefile +++ b/test/Makefile @@ -1,110 +1,152 @@ + ############################################################################ + # Copyright 2017 University of Warwick # + # # + # Licensed under the Apache License, Version 2.0 (the "License"); # + # you may not use this file except in compliance with the License. # + # You may obtain a copy of the License at # + # # + # http://www.apache.org/licenses/LICENSE-2.0 # + # # + # Unless required by applicable law or agreed to in writing, software # + # distributed under the License is distributed on an "AS IS" BASIS, # + # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # + # See the License for the specific language governing permissions and # + # limitations under the License. # + # # + # Laura++ package authors: # + # John Back # + # Paul Harrison # + # Thomas Latham # + # # + ############################################################################ + # # + # ------------------------------------- # + # Standalone Makefile for Laura++ tests # + # ------------------------------------- # + # # + # Instructions # + # - Review 'external configuration' section below # + # to match systems compilers setup # + # # + # - Make sure the ROOTSYS environment variable is set and points # + # to your ROOT release or the root-config script is in your PATH # + # # + # - run 'make ' # + # # + # Build targets # + # bin - make all examples (default) # + # - make the specific example # + # clean - delete all intermediate and final build objects # + # # + ############################################################################ + # --- External configuration ---------------------------------- # first check that ROOTSYS is defined ifndef ROOTSYS ROOTSYS := $(shell root-config --prefix) ROOTBINDIR := $(shell root-config --bindir) ifeq ($(ROOTSYS), ) $(error running of root-config failed or reported null value) endif else ROOTBINDIR := $(ROOTSYS)/bin endif # If you don't want to build one or more tests, just add them to the SKIPLIST variable SKIPLIST = ROOTCONFIG := $(ROOTBINDIR)/root-config ARCH := $(shell $(ROOTCONFIG) --arch) PLATFORM := $(shell $(ROOTCONFIG) --platform) INCLUDES = WORKDIR = tmp ifeq ($(findstring linux, $(ARCH)),linux) # This set here should work for Linux. CXX = g++ LD = g++ CXXFLAGS = -g -O2 -Wall -Wextra -Wshadow -Woverloaded-virtual -Werror -fPIC MFLAGS = -MM LDFLAGS = -g SOFLAGS = -shared endif ifeq ($(ARCH),macosx64) # For Mac OS X you may need to put -m64 in CXXFLAGS and SOFLAGS. CXX = g++ LD = g++ CXXFLAGS = -g -O3 -Wall -Wextra -Wshadow -Woverloaded-virtual -Werror -fPIC -m64 MFLAGS = -MM LDFLAGS = -g SOFLAGS = -m64 -dynamiclib -single_module -undefined dynamic_lookup endif # --- Internal configuration ---------------------------------- INCDIR=../inc DEPDIR=$(WORKDIR)/dependencies OBJDIR=$(WORKDIR)/objects ROOTCFLAGS := $(shell $(ROOTCONFIG) --cflags) ROOTLIBS := $(shell $(ROOTCONFIG) --libs) ROOTLIBS += -lEG ROOTLIBS += -lMinuit ROOTLIBS += -lTreePlayer ifeq ($(strip $(SKIPLIST)),) ROOTLIBS += -lRooFitCore ROOTLIBS += -lRooFit endif LAURALIBDIR=$(shell pwd | xargs dirname)/lib LAURALIB = $(LAURALIBDIR)/libLaura++.so INCLUDES += -I$(INCDIR) CXXFLAGS += $(INCLUDES) CXXFLAGS += $(ROOTCFLAGS) default: bin # List of all source files CCLIST:=$(filter-out $(SKIPLIST), $(wildcard *.cc)) # List of all source files that contain main functions BINCCLIST:=$(shell egrep -l "^[[:space:]]*int[[:space:]]*main\>" $(CCLIST)) # List of all object files to build BINOLIST:=$(patsubst %.cc,%.o,$(addprefix $(OBJDIR)/,$(notdir $(BINCCLIST)))) # List of all dependency files to make DLIST:=$(patsubst %.cc,%.d,$(addprefix $(DEPDIR)/,$(notdir $(CCLIST)))) # List of all binary files to make BINLIST:=$(patsubst %.cc,%,$(notdir $(BINCCLIST))) # Implicit rule making all dependency Makefiles included at the end of this makefile $(DEPDIR)/%.d: %.cc @echo "Making $@" @mkdir -p $(DEPDIR) @set -e; $(CXX) $(MFLAGS) $(CXXFLAGS) $< \ | sed 's#\($(notdir $*)\)\.o[ :]*#$(OBJDIR)/\1.o $@ : #g' > $@; \ [ -s $@ ] || rm -f $@ # Implicit rule to compile all sources $(OBJDIR)/%.o : %.cc @echo "Compiling $<" @mkdir -p $(OBJDIR) @$(CXX) $(CXXFLAGS) -c $< -o $@ # Rule to compile all binaries % : $(OBJDIR)/%.o $(LAURALIB) @echo "Linking $@" @$(CXX) $(LDFLAGS) $< -o $@ $(LAURALIB) $(ROOTLIBS) bin: $(BINLIST) clean: rm -f $(BINLIST) rm -rf $(WORKDIR) .PHONY : bin default clean -include $(DLIST) diff --git a/test/TestCovariant.cc b/test/TestCovariant.cc index e90dd6d..24005cd 100644 --- a/test/TestCovariant.cc +++ b/test/TestCovariant.cc @@ -1,140 +1,163 @@ +/* +Copyright 2017 University of Warwick + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +Laura++ package authors: +John Back +Paul Harrison +Thomas Latham +*/ + #include #include #include #include "TFile.h" #include "TH2.h" #include "TString.h" #include "TTree.h" #include "LauSimpleFitModel.hh" #include "LauBkgndDPModel.hh" #include "LauDaughters.hh" #include "LauEffModel.hh" #include "LauIsobarDynamics.hh" #include "LauMagPhaseCoeffSet.hh" #include "LauBelleNR.hh" #include "LauResonanceMaker.hh" #include "LauVetoes.hh" void usage( std::ostream& out, const TString& progName ) { out<<"Usage:\n"; out<"< if ( argc < 2 ) { usage( std::cerr, argv[0] ); return EXIT_FAILURE; } const Int_t spin = atoi(argv[1]); if ( spin < 0 || spin > 3 ) { std::cerr << "Only valid for spin 0, 1, 2 or 3" << std::endl; return EXIT_FAILURE; } // If you want to use square DP histograms for efficiency, // backgrounds or you just want the square DP co-ordinates // stored in the toy MC ntuple then set this to kTRUE Bool_t squareDP = kFALSE; // This defines the DP => decay is B+ -> K+ pi+ pi- // Particle 1 = K+ // Particle 2 = pi+ // Particle 3 = pi- // The DP is defined in terms of m13Sq and m23Sq LauDaughters* daughters = new LauDaughters("B+", "K+", "pi+", "pi-", squareDP); // Optionally apply some vetoes to the DP LauVetoes* vetoes = new LauVetoes(); // Define the efficiency model (defaults to unity everywhere) // Can optionally provide a histogram to model variation over DP // (example syntax given in commented-out section) LauEffModel* effModel = new LauEffModel(daughters, vetoes); // Create the isobar model // Set the values of the Blatt-Weisskopf barrier radii and whether they are fixed or floating LauResonanceMaker& resMaker = LauResonanceMaker::get(); resMaker.setSpinFormalism( LauAbsResonance::Covariant ); resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Parent, 5.0 ); resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Light, 4.0 ); resMaker.fixBWRadius( LauBlattWeisskopfFactor::Parent, kTRUE ); resMaker.fixBWRadius( LauBlattWeisskopfFactor::Light, kFALSE ); LauIsobarDynamics* sigModel = new LauIsobarDynamics(daughters, effModel); LauAbsResonance* reson(0); std::map waves; waves[0] = "S"; waves[1] = "P"; waves[2] = "D"; waves[3] = "F"; const TString nrName = TString::Format( "BelleNR_%swave", waves[spin].Data() ); reson = sigModel->addResonance(nrName, 2, LauAbsResonance::BelleNR); reson->setResonanceParameter("alpha", 0.0); LauBelleNR* belleNR = dynamic_cast( reson ); belleNR->enforceLegendreSpinFactors(kFALSE); // Reset the maximum signal DP ASq value // This will be automatically adjusted to avoid bias or extreme // inefficiency if you get the value wrong but best to set this by // hand once you've found the right value through some trial and // error. sigModel->setASqMaxValue(0.32); // Create the fit model LauSimpleFitModel* fitModel = new LauSimpleFitModel(sigModel); // Create the complex coefficients for the isobar model // Here we're using the magnitude and phase form: // c_j = a_j exp(i*delta_j) std::vector coeffset; coeffset.push_back( new LauMagPhaseCoeffSet(nrName, 1.00, 0.00, kTRUE, kTRUE) ); for (std::vector::iterator iter=coeffset.begin(); iter!=coeffset.end(); ++iter) { fitModel->setAmpCoeffSet(*iter); } // Set the signal yield and define whether it is fixed or floated const Double_t nSig = 100000.0; LauParameter * nSigEvents = new LauParameter("nSigEvents",nSig,-2.0*nSig,2.0*nSig,kFALSE); fitModel->setNSigEvents(nSigEvents); // Set the number of experiments to generate or fit and which // experiment to start with Int_t nExpt(1); Int_t firstExpt(0); fitModel->setNExpts( nExpt, firstExpt ); // Switch on/off calculation of asymmetric errors. fitModel->useAsymmFitErrors(kFALSE); // Randomise initial fit values for the signal mode fitModel->useRandomInitFitPars(kTRUE); const Bool_t haveBkgd = ( fitModel->nBkgndClasses() > 0 ); // Switch on/off Poissonian smearing of total number of events fitModel->doPoissonSmearing(haveBkgd); // Switch on/off Extended ML Fit option fitModel->doEMLFit(haveBkgd); // Set the names of the files to read/write TString dataFile = TString::Format("gen-covtest-spin%d.root", spin); TString treeName("genResults"); TString rootFileName("dummy.root"); TString tableFileName("genResults"); // Execute the generation/fit fitModel->run( "gen", dataFile, treeName, rootFileName, tableFileName ); return EXIT_SUCCESS; } diff --git a/test/TestCovariant2.cc b/test/TestCovariant2.cc index 0998e44..ded8e33 100644 --- a/test/TestCovariant2.cc +++ b/test/TestCovariant2.cc @@ -1,172 +1,195 @@ +/* +Copyright 2017 University of Warwick + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +Laura++ package authors: +John Back +Paul Harrison +Thomas Latham +*/ + #include #include #include #include "TFile.h" #include "TH2.h" #include "TString.h" #include "TTree.h" #include "LauSimpleFitModel.hh" #include "LauBkgndDPModel.hh" #include "LauDaughters.hh" #include "LauEffModel.hh" #include "LauIsobarDynamics.hh" #include "LauMagPhaseCoeffSet.hh" #include "LauResonanceMaker.hh" #include "LauVetoes.hh" void usage( std::ostream& out, const TString& progName ) { out<<"Usage:\n"; out< \n"; } int main( int argc, char** argv ) { // Process command-line arguments // Usage: // ./TestCovariant2 if ( argc < 4 ) { usage( std::cerr, argv[0] ); return EXIT_FAILURE; } TString resName = argv[1]; TString formalism = argv[2]; TString barrier = argv[3]; // If you want to use square DP histograms for efficiency, // backgrounds or you just want the square DP co-ordinates // stored in the toy MC ntuple then set this to kTRUE Bool_t squareDP = kFALSE; // This defines the DP => decay is B+ -> K+ pi+ pi- // Particle 1 = K+ // Particle 2 = pi+ // Particle 3 = pi- // The DP is defined in terms of m13Sq and m23Sq LauDaughters* daughters = new LauDaughters("B+", "K+", "pi+", "pi-", squareDP); // Optionally apply some vetoes to the DP LauVetoes* vetoes = new LauVetoes(); // Define the efficiency model (defaults to unity everywhere) // Can optionally provide a histogram to model variation over DP // (example syntax given in commented-out section) LauEffModel* effModel = new LauEffModel(daughters, vetoes); // Create the isobar model // Set the values of the Blatt-Weisskopf barrier radii and whether they are fixed or floating LauResonanceMaker& resMaker = LauResonanceMaker::get(); if ( formalism == "Covariant" ) { resMaker.setSpinFormalism( LauAbsResonance::Covariant ); } else if ( formalism == "Zemach_P" ) { resMaker.setSpinFormalism( LauAbsResonance::Zemach_P ); } else if ( formalism == "Zemach_Pstar" ) { resMaker.setSpinFormalism( LauAbsResonance::Zemach_Pstar ); } else if ( formalism == "Legendre" ) { resMaker.setSpinFormalism( LauAbsResonance::Legendre ); } else { std::cerr << "Unknown formalism: " << formalism << std::endl; return EXIT_FAILURE; } resMaker.setBWBachelorRestFrame( LauBlattWeisskopfFactor::ResonanceFrame ); if ( barrier == "noBWfactors" ) { resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Parent, 0.0 ); resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Light, 0.0 ); } else if ( barrier == "noNumeratorBWfactors" ) { resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Parent, 0.0 ); resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Light, 4.0 ); } else if ( barrier == "noParentNumeratorBWfactor" ) { resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Parent, 0.0 ); resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Light, 4.0 ); } else if ( barrier == "pParentNumeratorBWfactor" ) { resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Parent, 4.0 ); resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Light, 4.0 ); } else if ( barrier == "pstarParentNumeratorBWfactor" ) { resMaker.setBWBachelorRestFrame( LauBlattWeisskopfFactor::ParentFrame ); resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Parent, 4.0 ); resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Light, 4.0 ); } else if ( barrier == "covParentNumeratorBWfactor" ) { resMaker.setBWBachelorRestFrame( LauBlattWeisskopfFactor::Covariant ); resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Parent, 4.0 ); resMaker.setDefaultBWRadius( LauBlattWeisskopfFactor::Light, 4.0 ); } else { std::cerr << "Unknown barrier: " << barrier << std::endl; return EXIT_FAILURE; } LauIsobarDynamics* sigModel = new LauIsobarDynamics(daughters, effModel); LauAbsResonance* reson(0); Int_t resPairAmpInt(1); if ( resName.Contains("K*") ) { resPairAmpInt = 2; } reson = sigModel->addResonance(resName, resPairAmpInt, LauAbsResonance::RelBW); //reson->ignoreMomenta(kTRUE);// this is for turning off the mass-dependent width if ( barrier == "noBWfactors" || barrier == "noNumeratorBWfactors" ) { reson->ignoreBarrierScaling(kTRUE); } const TString sanitisedName = reson->getSanitisedName(); // Reset the maximum signal DP ASq value // This will be automatically adjusted to avoid bias or extreme // inefficiency if you get the value wrong but best to set this by // hand once you've found the right value through some trial and // error. sigModel->setASqMaxValue(0.45); // Create the fit model LauSimpleFitModel* fitModel = new LauSimpleFitModel(sigModel); // Create the complex coefficients for the isobar model // Here we're using the magnitude and phase form: // c_j = a_j exp(i*delta_j) std::vector coeffset; coeffset.push_back( new LauMagPhaseCoeffSet(resName, 1.00, 0.00, kTRUE, kTRUE) ); for (std::vector::iterator iter=coeffset.begin(); iter!=coeffset.end(); ++iter) { fitModel->setAmpCoeffSet(*iter); } // Set the signal yield and define whether it is fixed or floated const Double_t nSig = 100000.0; LauParameter * nSigEvents = new LauParameter("nSigEvents",nSig,-2.0*nSig,2.0*nSig,kFALSE); fitModel->setNSigEvents(nSigEvents); // Set the number of experiments to generate or fit and which // experiment to start with Int_t nExpt(1); Int_t firstExpt(0); fitModel->setNExpts( nExpt, firstExpt ); // Switch on/off calculation of asymmetric errors. fitModel->useAsymmFitErrors(kFALSE); // Randomise initial fit values for the signal mode fitModel->useRandomInitFitPars(kTRUE); const Bool_t haveBkgd = ( fitModel->nBkgndClasses() > 0 ); // Switch on/off Poissonian smearing of total number of events fitModel->doPoissonSmearing(haveBkgd); // Switch on/off Extended ML Fit option fitModel->doEMLFit(haveBkgd); // Set the names of the files to read/write TString dataFile = TString::Format( "%s/gen-covtest-%s-withMassDepWidth-%s.root", sanitisedName.Data(), barrier.Data(), formalism.Data() ); TString treeName("genResults"); TString rootFileName("dummy.root"); TString tableFileName("genResults"); // Execute the generation/fit fitModel->run( "gen", dataFile, treeName, rootFileName, tableFileName ); return EXIT_SUCCESS; } diff --git a/test/TestNewKinematicsMethods.cc b/test/TestNewKinematicsMethods.cc index 19ebea1..d12526a 100644 --- a/test/TestNewKinematicsMethods.cc +++ b/test/TestNewKinematicsMethods.cc @@ -1,304 +1,327 @@ +/* +Copyright 2017 University of Warwick + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +/* +Laura++ package authors: +John Back +Paul Harrison +Thomas Latham +*/ + #include #include #include #include "TFile.h" #include "TRandom.h" #include "TTree.h" #include "LauDaughters.hh" #include "LauKinematics.hh" #include "LauRandom.hh" struct DPpoint { DPpoint(); DPpoint(const LauKinematics* kine); void reset(); void set(const LauKinematics* kine); bool operator==(const DPpoint& rhs) const; Double_t m12Sq; Double_t m13Sq; Double_t m23Sq; Double_t m12; Double_t m13; Double_t m23; Double_t c12; Double_t c13; Double_t c23; Double_t mPrime; Double_t thPrime; Double_t p1_12; Double_t p3_12; Double_t p2_23; Double_t p1_23; Double_t p1_13; Double_t p2_13; Double_t p1_Parent; Double_t p2_Parent; Double_t p3_Parent; }; DPpoint::DPpoint() : m12Sq(0.0), m13Sq(0.0), m23Sq(0.0), m12(0.0), m13(0.0), m23(0.0), c12(0.0), c13(0.0), c23(0.0), mPrime(0.0), thPrime(0.0), p1_12(0.0), p3_12(0.0), p2_23(0.0), p1_23(0.0), p1_13(0.0), p2_13(0.0), p1_Parent(0.0), p2_Parent(0.0), p3_Parent(0.0) { } DPpoint::DPpoint( const LauKinematics* kine ) : m12Sq(kine->getm12Sq()), m13Sq(kine->getm13Sq()), m23Sq(kine->getm23Sq()), m12(kine->getm12()), m13(kine->getm13()), m23(kine->getm23()), c12(kine->getc12()), c13(kine->getc13()), c23(kine->getc23()), mPrime(kine->getmPrime()), thPrime(kine->getThetaPrime()), p1_12(kine->getp1_12()), p3_12(kine->getp3_12()), p2_23(kine->getp2_23()), p1_23(kine->getp1_23()), p1_13(kine->getp1_13()), p2_13(kine->getp2_13()), p1_Parent(kine->getp1_Parent()), p2_Parent(kine->getp2_Parent()), p3_Parent(kine->getp3_Parent()) { } void DPpoint::reset() { m12Sq = 0.0; m13Sq = 0.0; m23Sq = 0.0; m12 = 0.0; m13 = 0.0; m23 = 0.0; c12 = 0.0; c13 = 0.0; c23 = 0.0; mPrime = 0.0; thPrime = 0.0; p1_12 = 0.0; p3_12 = 0.0; p2_23 = 0.0; p1_23 = 0.0; p1_13 = 0.0; p2_13 = 0.0; p1_Parent = 0.0; p2_Parent = 0.0; p3_Parent = 0.0; } void DPpoint::set( const LauKinematics* kine ) { m12Sq = kine->getm12Sq(); m13Sq = kine->getm13Sq(); m23Sq = kine->getm23Sq(); m12 = kine->getm12(); m13 = kine->getm13(); m23 = kine->getm23(); c12 = kine->getc12(); c13 = kine->getc13(); c23 = kine->getc23(); mPrime = kine->getmPrime(); thPrime = kine->getThetaPrime(); p1_12 = kine->getp1_12(); p3_12 = kine->getp3_12(); p2_23 = kine->getp2_23(); p1_23 = kine->getp1_23(); p1_13 = kine->getp1_13(); p2_13 = kine->getp2_13(); p1_Parent = kine->getp1_Parent(); p2_Parent = kine->getp2_Parent(); p3_Parent = kine->getp3_Parent(); } bool DPpoint::operator==( const DPpoint& rhs ) const { bool ok = true; if ( TMath::Abs( m12Sq - rhs.m12Sq ) > 1e-6 ) { std::cerr << "Problem with m12Sq: " << m12Sq << " != " << rhs.m12Sq << std::endl; ok = false; } if ( TMath::Abs( m13Sq - rhs.m13Sq ) > 1e-6 ) { std::cerr << "Problem with m13Sq: " << m13Sq << " != " << rhs.m13Sq << std::endl; ok = false; } if ( TMath::Abs( m23Sq - rhs.m23Sq ) > 1e-6 ) { std::cerr << "Problem with m23Sq: " << m23Sq << " != " << rhs.m23Sq << std::endl; ok = false; } if ( TMath::Abs( m12 - rhs.m12 ) > 1e-6 ) { std::cerr << "Problem with m12: " << m12 << " != " << rhs.m12 << std::endl; ok = false; } if ( TMath::Abs( m13 - rhs.m13 ) > 1e-6 ) { std::cerr << "Problem with m13: " << m13 << " != " << rhs.m13 << std::endl; ok = false; } if ( TMath::Abs( m23 - rhs.m23 ) > 1e-6 ) { std::cerr << "Problem with m23: " << m23 << " != " << rhs.m23 << std::endl; ok = false; } if ( TMath::Abs( c12 - rhs.c12 ) > 1e-6 ) { std::cerr << "Problem with c12: " << c12 << " != " << rhs.c12 << std::endl; ok = false; } if ( TMath::Abs( c13 - rhs.c13 ) > 1e-6 ) { std::cerr << "Problem with c13: " << c13 << " != " << rhs.c13 << std::endl; ok = false; } if ( TMath::Abs( c23 - rhs.c23 ) > 1e-6 ) { std::cerr << "Problem with c23: " << c23 << " != " << rhs.c23 << std::endl; ok = false; } if ( TMath::Abs( mPrime - rhs.mPrime ) > 1e-6 ) { std::cerr << "Problem with mPrime: " << mPrime << " != " << rhs.mPrime << std::endl; ok = false; } if ( TMath::Abs( thPrime - rhs.thPrime ) > 1e-6 ) { std::cerr << "Problem with thPrime: " << thPrime << " != " << rhs.thPrime << std::endl; ok = false; } if ( TMath::Abs( p1_12 - rhs.p1_12 ) > 1e-6 ) { std::cerr << "Problem with p1_12: " << p1_12 << " != " << rhs.p1_12 << std::endl; ok = false; } if ( TMath::Abs( p3_12 - rhs.p3_12 ) > 1e-6 ) { std::cerr << "Problem with p3_12: " << p3_12 << " != " << rhs.p3_12 << std::endl; ok = false; } if ( TMath::Abs( p2_23 - rhs.p2_23 ) > 1e-6 ) { std::cerr << "Problem with p2_23: " << p2_23 << " != " << rhs.p2_23 << std::endl; ok = false; } if ( TMath::Abs( p1_23 - rhs.p1_23 ) > 1e-6 ) { std::cerr << "Problem with p1_23: " << p1_23 << " != " << rhs.p1_23 << std::endl; ok = false; } if ( TMath::Abs( p1_13 - rhs.p1_13 ) > 1e-6 ) { std::cerr << "Problem with p1_13: " << p1_13 << " != " << rhs.p1_13 << std::endl; ok = false; } if ( TMath::Abs( p2_13 - rhs.p2_13 ) > 1e-6 ) { std::cerr << "Problem with p2_13: " << p2_13 << " != " << rhs.p2_13 << std::endl; ok = false; } if ( TMath::Abs( p1_Parent - rhs.p1_Parent ) > 1e-6 ) { std::cerr << "Problem with p1_Parent: " << p1_Parent << " != " << rhs.p1_Parent << std::endl; ok = false; } if ( TMath::Abs( p2_Parent - rhs.p2_Parent ) > 1e-6 ) { std::cerr << "Problem with p2_Parent: " << p2_Parent << " != " << rhs.p2_Parent << std::endl; ok = false; } if ( TMath::Abs( p3_Parent - rhs.p3_Parent ) > 1e-6 ) { std::cerr << "Problem with p3_Parent: " << p3_Parent << " != " << rhs.p3_Parent << std::endl; ok = false; } return ok; } int main( /*int argc, char** argv*/ ) { // If you want to use square DP histograms for efficiency, // backgrounds or you just want the square DP co-ordinates // stored in the toy MC ntuple then set this to kTRUE Bool_t squareDP = kTRUE; // This defines the DP => decay is B+ -> pi+ pi+ pi- // Particle 1 = pi+ // Particle 2 = pi+ // Particle 3 = pi- // The DP is defined in terms of m13Sq and m23Sq LauDaughters* daughters = new LauDaughters("B+", "pi+", "pi+", "pi-", squareDP); LauKinematics* kinematics = daughters->getKinematics(); Double_t m13Sq(0.0), m23Sq(0.0); DPpoint pos1; DPpoint pos2; bool ok(true); for ( UInt_t i(0); i < 10; ++i ) { pos1.reset(); pos2.reset(); kinematics->genFlatPhaseSpace( m13Sq, m23Sq ); kinematics->updateKinematics( m13Sq, m23Sq ); pos1.set(kinematics); const Double_t m12 = kinematics->getm12(); const Double_t c12 = kinematics->getc12(); const Double_t m13 = kinematics->getm13(); const Double_t c13 = kinematics->getc13(); const Double_t m23 = kinematics->getm23(); const Double_t c23 = kinematics->getc23(); // this is just to make sure that we reset all variables so we can notice if things don't get set properly kinematics->genFlatPhaseSpace( m13Sq, m23Sq ); kinematics->updateKinematics( m13Sq, m23Sq ); kinematics->updateKinematicsFrom12( m12, c12 ); pos2.set(kinematics); ok = pos1 == pos2; if ( ! ok ) { std::cerr << "Problem setting from m12, c12" << std::endl; } // this is just to make sure that we reset all variables so we can notice if things don't get set properly kinematics->genFlatPhaseSpace( m13Sq, m23Sq ); kinematics->updateKinematics( m13Sq, m23Sq ); kinematics->updateKinematicsFrom13( m13, c13 ); pos2.set(kinematics); ok = pos1 == pos2; if ( ! ok ) { std::cerr << "Problem setting from m13, c13" << std::endl; } // this is just to make sure that we reset all variables so we can notice if things don't get set properly kinematics->genFlatPhaseSpace( m13Sq, m23Sq ); kinematics->updateKinematics( m13Sq, m23Sq ); kinematics->updateKinematicsFrom23( m23, c23 ); pos2.set(kinematics); ok = pos1 == pos2; if ( ! ok ) { std::cerr << "Problem setting from m23, c23" << std::endl; } } if ( ! ok ) { return EXIT_FAILURE; } return EXIT_SUCCESS; }