diff --git a/EvtGenModels/EvtItgAbsFunction.hh b/EvtGenModels/EvtItgAbsFunction.hh index 1d60026..0274b68 100644 --- a/EvtGenModels/EvtItgAbsFunction.hh +++ b/EvtGenModels/EvtItgAbsFunction.hh @@ -1,72 +1,69 @@ //-------------------------------------------------------------------------- // // // Copyright Information: See EvtGen/COPYRIGHT // // 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. // // Module: EvtItgAbsFunction.hh // // Description: // Abstraction of a generic function for use in integration methods elsewhere // in this package. (Stolen and modified from the BaBar IntegrationUtils package // - author: Phil Strother). // // Modification history: // // Jane Tinslay March 21, 2001 Module adapted for use in // EvtGen // //------------------------------------------------------------------------ #ifndef EVTITGABSFUNCTION_HH #define EVTITGABSFUNCTION_HH //------------- // C Headers -- //------------- extern "C" { } class EvtItgAbsFunction { public: // Constructors EvtItgAbsFunction(double lowerRange, double upperRange); // Destructor virtual ~EvtItgAbsFunction( ) = default; virtual double value( double x) const; virtual double operator()(double x) const; // Selectors (const) inline double upperRange() const {return _upperRange;} inline double lowerRange() const {return _lowerRange;} inline void getRange(double &lower,double &upper) const { lower = _lowerRange; upper = _upperRange; } virtual void setCoeff(int, int, double)=0; virtual double getCoeff(int, int)=0; protected: virtual double myFunction(double x) const=0; void setRange(double x1,double x2) { _lowerRange=x1; _upperRange=x2; }; private: double _upperRange; double _lowerRange; - EvtItgAbsFunction( const EvtItgAbsFunction& ); // Copy Constructor - EvtItgAbsFunction& operator= ( const EvtItgAbsFunction& ); // Assignment op - }; #endif // EVTITGABSFUNCTION_HH diff --git a/EvtGenModels/EvtItgFourCoeffFcn.hh b/EvtGenModels/EvtItgFourCoeffFcn.hh index b024e54..b262d21 100644 --- a/EvtGenModels/EvtItgFourCoeffFcn.hh +++ b/EvtGenModels/EvtItgFourCoeffFcn.hh @@ -1,58 +1,54 @@ //-------------------------------------------------------------------------- // // // Copyright Information: See EvtGen/COPYRIGHT // // 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. // // Module: EvtItgFourCoeffFcn.hh // // Description: // Class describing a function with Four vectors of coefficients. // // Modification history: // // Jane Tinslay March 21, 2001 Module created // //------------------------------------------------------------------------ #ifndef EVTITFOURCOEFFFCN_HH #define EVTITFOURCOEFFFCN_HH #include #include "EvtGenModels/EvtItgAbsFunction.hh" class EvtItgFourCoeffFcn: public EvtItgAbsFunction { public: EvtItgFourCoeffFcn( double (*theFunction)(double, const std::vector &, const std::vector &, const std::vector &, const std::vector &), double lowerRange, double upperRange, const std::vector &coeffs1, const std::vector &coeffs2, const std::vector &coeffs3, const std::vector &coeffs4); void setCoeff(int, int, double) override; double getCoeff(int, int) override; protected: double myFunction(double x) const override; private: // Data members double (*_myFunction)(double x, const std::vector & coeffs1, const std::vector & coeffs2, const std::vector & coeffs3, const std::vector & coeffs4); - // Note: if your class needs a copy constructor or an assignment operator, - // make one of the following public and implement it. - EvtItgFourCoeffFcn( const EvtItgFourCoeffFcn& ); //// Copy Constructor - EvtItgFourCoeffFcn& operator= ( const EvtItgFourCoeffFcn& ); // Assignment op std::vector _coeffs1; std::vector _coeffs2; std::vector _coeffs3; std::vector _coeffs4; }; #endif // EvtITGPTRFUNCTION_HH diff --git a/EvtGenModels/EvtItgFunction.hh b/EvtGenModels/EvtItgFunction.hh index fba807d..c3b4c39 100644 --- a/EvtGenModels/EvtItgFunction.hh +++ b/EvtGenModels/EvtItgFunction.hh @@ -1,62 +1,58 @@ //-------------------------------------------------------------------------- // // Environment: // This software was developed for the BaBar collaboration. If you // use all or part of it, please give an appropriate acknowledgement. // // Copyright Information: See EvtGen/COPYRIGHT // Copyright (C) 1998 LBNL // //------------------------------------------------------------------------ #ifndef EVTITGFUNCTION_HH #define EVTITGFUNCTION_HH #include "EvtGenModels/EvtItgAbsFunction.hh" /** * Copyright (C) 1998 LBNL * * Generic function where the pointer to the function is available. * * The function is taken as type pointer to function returning double and * taking a double (the abscissa) and a const RWTValVector reference * (the parameter values of the function) as arguments. * * @see EvtItgFunctionEvtItgFunction * * @version $Id: EvtItgFunction.hh,v 1.2 2009-03-16 16:34:00 robbep Exp $ * * @author Phil Strother Originator */ class EvtItgFunction: public EvtItgAbsFunction { public: // Constructors EvtItgFunction( double (*theFunction)(double), double lowerRange, double upperRange); void setCoeff(int, int, double) override {}; double getCoeff(int, int) override {return 0.0;}; protected: // Helper functions double myFunction(double x) const override; private: // Data members double (*_myFunction)(double x); - // Note: if your class needs a copy constructor or an assignment operator, - // make one of the following public and implement it. - EvtItgFunction( const EvtItgFunction& ); // Copy Constructor - EvtItgFunction& operator= ( const EvtItgFunction& ); // Assignment op }; #endif // EvtITGFUNCTION_HH diff --git a/EvtGenModels/EvtItgPtrFunction.hh b/EvtGenModels/EvtItgPtrFunction.hh index cea27d4..8dd3a69 100644 --- a/EvtGenModels/EvtItgPtrFunction.hh +++ b/EvtGenModels/EvtItgPtrFunction.hh @@ -1,57 +1,53 @@ //-------------------------------------------------------------------------- // // // Copyright Information: See EvtGen/COPYRIGHT // // 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. // // Module: EvtItgPtrFunction.hh // // Description: // Class describing a function with one vector of coefficients. (Stolen and // modified from the BaBar IntegrationUtils package - author: Phil Strother). // // Modification history: // // Jane Tinslay March 21, 2001 Module adapted for use in // EvtGen // //------------------------------------------------------------------------ #ifndef EVTITGPTRFUNCTION_HH #define EVTITGPTRFUNCTION_HH #include #include "EvtGenModels/EvtItgAbsFunction.hh" class EvtItgPtrFunction: public EvtItgAbsFunction { public: EvtItgPtrFunction( double (*theFunction)(double, const std::vector &), double lowerRange, double upperRange, const std::vector &coeffs1); void setCoeff(int, int, double) override; double getCoeff(int, int) override; protected: double myFunction(double x) const override; private: // Data members double (*_myFunction)(double x, const std::vector & coeffs1); - // Note: if your class needs a copy constructor or an assignment operator, - // make one of the following public and implement it. - EvtItgPtrFunction( const EvtItgPtrFunction& ); //// Copy Constructor - EvtItgPtrFunction& operator= ( const EvtItgPtrFunction& ); // Assignment op std::vector _coeffs1; }; #endif // EVTITGPTRFUNCTION_HH diff --git a/EvtGenModels/EvtItgThreeCoeffFcn.hh b/EvtGenModels/EvtItgThreeCoeffFcn.hh index 37291d0..65c3517 100644 --- a/EvtGenModels/EvtItgThreeCoeffFcn.hh +++ b/EvtGenModels/EvtItgThreeCoeffFcn.hh @@ -1,65 +1,61 @@ //-------------------------------------------------------------------------- // // // Copyright Information: See EvtGen/COPYRIGHT // // 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. // // Module: EvtItgThreeCoeffFcn.hh // // Description: // Class describing a function with three vectors of coefficients. // // Modification history: // // Jane Tinslay March 21, 2001 Module created // //------------------------------------------------------------------------ #ifndef EVTITTHREECOEFFFCN_HH #define EVTITTHREECOEFFFCN_HH #include //------------- // C Headers -- //------------- extern "C" { } #include "EvtGenModels/EvtItgAbsFunction.hh" class EvtItgThreeCoeffFcn: public EvtItgAbsFunction { public: EvtItgThreeCoeffFcn( double (*theFunction)(double, const std::vector &, const std::vector &, const std::vector &), double lowerRange, double upperRange, const std::vector &coeffs1, const std::vector &coeffs2, const std::vector &coeffs3); void setCoeff(int, int, double) override; double getCoeff(int, int) override; protected: double myFunction(double x) const override; private: // Data members double (*_myFunction)(double x, const std::vector & coeffs1, const std::vector & coeffs2, const std::vector & coeffs3); - // Note: if your class needs a copy constructor or an assignment operator, - // make one of the following public and implement it. - EvtItgThreeCoeffFcn( const EvtItgThreeCoeffFcn& ); //// Copy Constructor - EvtItgThreeCoeffFcn& operator= ( const EvtItgThreeCoeffFcn& ); // Assignment op std::vector _coeffs1; std::vector _coeffs2; std::vector _coeffs3; }; #endif // EvtITGPTRFUNCTION_HH diff --git a/EvtGenModels/EvtItgTwoCoeffFcn.hh b/EvtGenModels/EvtItgTwoCoeffFcn.hh index f827f4b..218e25f 100644 --- a/EvtGenModels/EvtItgTwoCoeffFcn.hh +++ b/EvtGenModels/EvtItgTwoCoeffFcn.hh @@ -1,64 +1,59 @@ //-------------------------------------------------------------------------- // // // Copyright Information: See EvtGen/COPYRIGHT // // 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. // // Module: EvtItgTwoCoeffFcn.hh // // Description: // Class describing a function with two vectors of coefficients. // // Modification history: // // Jane Tinslay March 21, 2001 Module created // //------------------------------------------------------------------------ #ifndef EVTITTWOCOEFFFCN_HH #define EVTITTWOCOEFFFCN_HH #include //------------- // C Headers -- //------------- extern "C" { } #include "EvtGenModels/EvtItgAbsFunction.hh" class EvtItgTwoCoeffFcn: public EvtItgAbsFunction { public: EvtItgTwoCoeffFcn( double (*theFunction)(double, const std::vector &, const std::vector &), double lowerRange, double upperRange, const std::vector &coeffs1, const std::vector &coeffs2); void setCoeff(int, int, double) override; double getCoeff(int, int) override; protected: double myFunction(double x) const override; private: // Data members double (*_myFunction)(double x, const std::vector & coeffs1, const std::vector & coeffs2); - // Note: if your class needs a copy constructor or an assignment operator, - // make one of the following public and implement it. - EvtItgTwoCoeffFcn( const EvtItgTwoCoeffFcn& ); //// Copy Constructor - EvtItgTwoCoeffFcn& operator= ( const EvtItgTwoCoeffFcn& ); // Assignment op - std::vector _coeffs1; std::vector _coeffs2; }; #endif // EvtITGTWOCOEFFFUNCTION_HH diff --git a/EvtGenModels/EvtNoRadCorr.hh b/EvtGenModels/EvtNoRadCorr.hh index 9add671..bb4a2d2 100644 --- a/EvtGenModels/EvtNoRadCorr.hh +++ b/EvtGenModels/EvtNoRadCorr.hh @@ -1,45 +1,45 @@ //-------------------------------------------------------------------------- // // 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) 2012 University of Warwick, UK // // Module: EvtNoRadCorr // // Description: Create an empty radiative correction engine which does nothing. // This is required since the EvtGen constructor still needs at least one // concrete implementation of EvtAbsRadCorr for the case when Photos is not used. // // Modification history: // // John Back Sept 2012 Module created // //------------------------------------------------------------------------------ // #ifndef EVTNORADCORR_HH #define EVTNORADCORR_HH #include "EvtGenBase/EvtAbsRadCorr.hh" #include class EvtParticle; class EvtNoRadCorr : public EvtAbsRadCorr { public: EvtNoRadCorr() {;} virtual ~EvtNoRadCorr() {;} - virtual void doRadCorr(EvtParticle *p) {;} + virtual void doRadCorr(EvtParticle *) {;} private: }; #endif diff --git a/EvtGenModels/EvtSTSCP.hh b/EvtGenModels/EvtSTSCP.hh index 14bbe21..7b63dae 100644 --- a/EvtGenModels/EvtSTSCP.hh +++ b/EvtGenModels/EvtSTSCP.hh @@ -1,44 +1,44 @@ //-------------------------------------------------------------------------- // // 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: EvtGen/EvtSTSCP.hh // // Description: // // Modification history: // // DJL/RYD August 11, 1998 Module created // //------------------------------------------------------------------------ #ifndef EVTSTSCP_HH #define EVTSTSCP_HH #include "EvtGenBase/EvtDecayAmp.hh" class EvtParticle; class EvtSTSCP:public EvtDecayAmp { public: std::string getName() override; EvtDecayBase* clone() override; void init() override; void initProbMax() override; void decay(EvtParticle *p) override; std::string getParamName(int i) override; - std::string getDefaultName(int i) override; + std::string getParamDefault(int i) override; }; #endif diff --git a/EvtGenModels/EvtTVSPwave.hh b/EvtGenModels/EvtTVSPwave.hh index 6726485..6528483 100644 --- a/EvtGenModels/EvtTVSPwave.hh +++ b/EvtGenModels/EvtTVSPwave.hh @@ -1,50 +1,49 @@ //-------------------------------------------------------------------------- // // 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: EvtGen/EvtTVSPwave.hh // // Description: // // Modification history: // // DJL/RYD August 11, 1998 Module created // //------------------------------------------------------------------------ #ifndef EVTTVSPWAVE_HH #define EVTTVSPWAVE_HH #include "EvtGenBase/EvtDecayAmp.hh" class EvtParticle; //Class to handles decays of the form TENSOR -> VECTOR SCALAR, //which proceed via S,P, or D partial waves. There are //six arguements, which are the magnetude and phase of //each partial wave (S, P then D). Calls EvtTensorToVectorScalar class EvtTVSPwave:public EvtDecayAmp { public: std::string getName() override; EvtDecayBase* clone() override; void decay(EvtParticle *p) override; void init() override; void initProbMax() override; std::string getParamName(int i) override; std::string getParamDefault(int i) override; -}; }; #endif diff --git a/History.txt b/History.txt index b015629..ccba384 100644 --- a/History.txt +++ b/History.txt @@ -1,561 +1,564 @@ //========================================================================== // // History file for EvtGen // //=========================================================================== +30 Jan 2019 John Back + Fix modernization compiler errors and warnings. + 29 Jan 2019 Michal Kreps - Allow reading decay files which are missing end-of-line before end-of-file. + Allow reading decay files which are missing end-of-line before end-of-file. 21 & 24 Jan 2019 John Back Correct B 4-momentum in EvtSLDiBaryonAmp for the BToDiBaryonlnupQCD model, as well as including a parity factor for the vector amplitude terms. 21st December 2018 John Back Imported C++ modernization changes from Gerhard Raven (LHCb). 7th December 2018 John Back Added the EvtBLLNuL (BLLNUL) model that generates rare B -> ell ell nu ell decays, where ell = e or mu, courtesy of Anna Danilina and Nikolai Nikitin (LHCb). Removed the EvtB2MuMuMuNu (BUTOMMMN) model, since its now replaced by the more general BLLNuL one. 5th November 2018 John Back Added the BToDiBaryonlnupQCD model for generating B to p N* l nu decays, where N can be any (exited) charged baryon (spin 1/2 or 3/2), courtesy of Mark Smith and Ryan Newcombe (LHCb), with added code optimisations. 17th October 2018 John Back Added various decay models from LHCb EvtGenExtras package: EvtBcVHad ("BC_VHAD"), Evtbs2llGammaMNT ("BSTOGLLMNT"), Evtbs2llGammaISRFSR ("BSTOGLLISRFSR"), EvtbTosllMS ("BTOSLLMS"), EvtbTosllMSExt ("BTOSLLMSEXT"), EvtLb2Baryonlnu ("Lb2Baryonlnu"), EvtLb2plnuLCSR ("Lb2plnuLCSR"), EvtLb2plnuLQCD ("Lb2plnuLQCD"), EvtFlatSqDalitz ("FLATSQDALITZ") and EvtPhspFlatLifetime ("PHSPFLATLIFETIME"). 5th October 2018 John Back Updated setupEvtGen.sh to work with the new HepForge Phabricator site. 13th March 2018 John Back Updated EvtPythiaEngine to correctly handle updates of various particle properties so that Pythia uses the same information as EvtGen (evt.pdl) for the generic and alias PYTHIA decay model. 12th March 2018 John Back Updated EvtBcXMuNu models (X = Scalar, Vector, Tensor) to generate Bc to D0(star) mu nu decays, with associated form factors in EvtBCXFF, courtesy of Aleksei Luchinskii (LHCb). Also generalised the calculation of their maximum probabilities by reusing the CalcMaxProb method in EvtSemiLeptonicAmp, which now allows for different Q^2 binning (default remains at 25 bins). R01-07-00------------------------------------------------------------------- 13th December 2017 John Back New tag incorporating all changes below. Recommended external packages are HepMC 2.06.09, pythia 8.230, Photos++ 3.61 and Tauola++ 1.1.6c, as used in the setupEvtGen.sh script. 12th December 2017 John Back Changed Pythia validation example decay files to use Pythia8 codes. 6th December 2017 John Back Modified the examples to use DECAY.DEC (see 25th April 2016) instead of DECAY_2010.DEC. Changed EvtExternalGenList to assume Pythia8 codes are used in decay files by default, which is the case for DECAY.DEC. Also updated the setupEvtGen.sh script to work with Pythia 8.2x versions. 29th November 2017 John Back Modified EvtSVP, EvtVVP and EvtTVP models to handle both radiative and two-lepton decays, courtesy of Aleksei Luchinskii (LHCb). 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. 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/src/EvtGenBase/EvtGammaMatrix.cpp b/src/EvtGenBase/EvtGammaMatrix.cpp index ff3814b..7e88ca9 100644 --- a/src/EvtGenBase/EvtGammaMatrix.cpp +++ b/src/EvtGenBase/EvtGammaMatrix.cpp @@ -1,698 +1,696 @@ //-------------------------------------------------------------------------- // // 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: EvtGammaMatrix.cc // // Description: Make gamma matrices availible for the calc. of amplitudes, etc. // // Modification history: // // DJL/RYD September 25, 1996 Module created // //------------------------------------------------------------------------ // #include "EvtGenBase/EvtPatches.hh" #include #include #include #include "EvtGenBase/EvtComplex.hh" #include "EvtGenBase/EvtGammaMatrix.hh" #include "EvtGenBase/EvtDiracSpinor.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtTensor4C.hh" #include "EvtGenBase/EvtVector4C.hh" #include using std::endl; using std::ostream; EvtGammaMatrix::EvtGammaMatrix(){ int i,j; static EvtComplex zero(0.0,0.0); for(i=0;i<4;i++){ for(j=0;j<4;j++){ _gamma[i][j]=zero; } } } EvtGammaMatrix operator*(const EvtGammaMatrix& g, const EvtComplex& c) { return c*g; } EvtGammaMatrix operator*(const EvtComplex& c,const EvtGammaMatrix& g){ int i,j; EvtGammaMatrix temp; for(i=0;i<4;i++){ for(j=0;j<4;j++){ temp._gamma[i][j]=g._gamma[i][j]*c; } } return temp; } ostream& operator<<(ostream& s, const EvtGammaMatrix& g){ s<<"["< 3 || nu > 3) { EvtGenReport(EVTGEN_ERROR, "EvtSigmaTensor") << "Expected index between 0 and 3, but found " << nu << "!" << endl; assert(0); } return sigma[mu][nu]; } const EvtGammaMatrix& EvtGammaMatrix::sigmaLower(unsigned int mu, unsigned int nu) { const EvtComplex I(0, 1); EvtGammaMatrix a, b; static EvtGammaMatrix sigma[4][4]; static bool hasBeenCalled = false; static const EvtTensor4C eta = EvtTensor4C::g(); if (!hasBeenCalled) // has to be initialized only at the first call { // lower index for (int i=0; i<4; ++i) { a = eta.get(i, 0)*g0() + eta.get(i, 1)*g1() + eta.get(i, 2)*g2() + eta.get(i, 3)*g3(); for (int j=0; j<4; ++j) { b = eta.get(j, 0)*g0() + eta.get(j, 1)*g1() + eta.get(j, 2)*g2() + eta.get(j, 3)*g3(); sigma[i][j] = I/2 * (a*b - b*a); } } } return sigma[mu][nu]; } EvtGammaMatrix EvtGenFunctions::slash(const EvtVector4C& p) { return EvtGammaMatrix::g0()*p.get(0) + EvtGammaMatrix::g1()*p.get(1) + EvtGammaMatrix::g2()*p.get(2) + EvtGammaMatrix::g3()*p.get(3); } EvtGammaMatrix EvtGenFunctions::slash(const EvtVector4R& p) { return EvtGammaMatrix::g0()*p.get(0) + EvtGammaMatrix::g1()*p.get(1) + EvtGammaMatrix::g2()*p.get(2) + EvtGammaMatrix::g3()*p.get(3); } diff --git a/src/EvtGenModels/EvtRareLbToLll.cpp b/src/EvtGenModels/EvtRareLbToLll.cpp index ec66fe9..fcfa319 100644 --- a/src/EvtGenModels/EvtRareLbToLll.cpp +++ b/src/EvtGenModels/EvtRareLbToLll.cpp @@ -1,494 +1,495 @@ //-------------------------------------------------------------------------- // // 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) 2003 Caltech, UCSB // // Module: EvtRareLbToLll // // Description: // Implements the rare Lb --> Lambda^(*) ell ell models described in // http://arxiv.org/pdf/1108.6129.pdf // // Modification history: // // T. Blake November 2013 Created module // //------------------------------------------------------------------------ // #include #include "EvtGenModels/EvtRareLbToLll.hh" #include "EvtGenModels/EvtRareLbToLllFF.hh" #include "EvtGenModels/EvtRareLbToLllFFGutsche.hh" #include "EvtGenModels/EvtRareLbToLllFFlQCD.hh" #include "EvtGenBase/EvtComplex.hh" #include "EvtGenBase/EvtVector4C.hh" #include "EvtGenBase/EvtVector4R.hh" #include "EvtGenBase/EvtTensor4C.hh" #include "EvtGenBase/EvtDiracSpinor.hh" #include "EvtGenBase/EvtDiracParticle.hh" #include "EvtGenBase/EvtRaritaSchwinger.hh" #include "EvtGenBase/EvtSpinDensity.hh" #include "EvtGenBase/EvtPDL.hh" // The module name specification std::string EvtRareLbToLll::getName( ) { return "RareLbToLll" ; } // The implementation of the clone() method EvtDecayBase* EvtRareLbToLll::clone(){ return new EvtRareLbToLll; } void EvtRareLbToLll::init(){ checkNArg(1); // check that there are 3 daughteres checkNDaug(3); // Parent should be spin 1/2 Lambda_b0 const EvtSpinType::spintype spin = EvtPDL::getSpinType(getDaug(0)); if ( !( spin == EvtSpinType::DIRAC || spin == EvtSpinType::RARITASCHWINGER ) ) { EvtGenReport(EVTGEN_ERROR,"EvtGen") << " EvtRareLbToLll expects DIRAC or RARITASWINGER daughter " << std::endl; } // We expect that the second and third daughters // are the ell+ and ell- checkSpinDaughter(1,EvtSpinType::DIRAC); checkSpinDaughter(2,EvtSpinType::DIRAC); std::string model = getArgStr(0); if ( model == "Gutsche" ) { ffmodel_ = std::make_unique(); } else if ( model == "LQCD" ) { ffmodel_ = std::make_unique(); } else if ( model == "MR" ) { ffmodel_ = std::make_unique(); } else { EvtGenReport(EVTGEN_INFO ,"EvtGen") << " Unknown form-factor model, valid options are MR, LQCD, Gutsche." << std::endl; ::abort(); } wcmodel_ = std::make_unique(); ffmodel_->init(); return; } void EvtRareLbToLll::initProbMax(){ EvtGenReport(EVTGEN_INFO,"EvtGen") << " EvtRareLbToLll is finding maximum probability ... " << std::endl; m_maxProbability=0; if(m_maxProbability==0){ - auto parent = EvtDiracParticle{}; - parent.noLifeTime(); - parent.init(getParentId(),EvtVector4R(EvtPDL::getMass(getParentId()),0,0,0)); - parent.setDiagonalSpinDensity(); + EvtDiracParticle* parent = new EvtDiracParticle(); + parent->noLifeTime(); + parent->init(getParentId(),EvtVector4R(EvtPDL::getMass(getParentId()),0,0,0)); + parent->setDiagonalSpinDensity(); EvtAmp amp; EvtId daughters[3] = {getDaug(0),getDaug(1),getDaug(2)}; amp.init(getParentId(),3,daughters); - parent.makeDaughters(3,daughters); - EvtParticle *lambda = parent.getDaug(0); - EvtParticle *lep1 = parent.getDaug(1); - EvtParticle *lep2 = parent.getDaug(2); + parent->makeDaughters(3,daughters); + EvtParticle *lambda = parent->getDaug(0); + EvtParticle *lep1 = parent->getDaug(1); + EvtParticle *lep2 = parent->getDaug(2); lambda -> noLifeTime(); lep1 -> noLifeTime(); lep2 -> noLifeTime(); EvtSpinDensity rho; - rho.setDiag(parent.getSpinStates()); + rho.setDiag(parent->getSpinStates()); double M0 = EvtPDL::getMass(getParentId()); double mL = EvtPDL::getMass(getDaug(0)); double m1 = EvtPDL::getMass(getDaug(1)); double m2 = EvtPDL::getMass(getDaug(2)); double q2,pstar,elambda,theta; double q2min = (m1+m2)*(m1+m2); double q2max = (M0-mL)*(M0-mL); EvtVector4R p4lambda,p4lep1,p4lep2,boost; EvtGenReport(EVTGEN_INFO,"EvtGen") << " EvtRareLbToLll is probing whole phase space ..." << std::endl; int i,j; double prob=0; for(i=0;i<=100;i++){ q2 = q2min+i*(q2max-q2min)/100.; elambda = (M0*M0+mL*mL-q2)/2/M0; if(i==0) pstar = 0; else pstar = sqrt(q2-(m1+m2)*(m1+m2))*sqrt(q2-(m1-m2)*(m1-m2))/2/sqrt(q2); boost.set(M0-elambda,0,0,+sqrt(elambda*elambda-mL*mL)); if ( i != 100 ) { p4lambda.set(elambda,0,0,-sqrt(elambda*elambda-mL*mL)); } else { p4lambda.set(mL,0,0,0); } for(j=0;j<=45;j++){ theta = j*EvtConst::pi/45; p4lep1.set(sqrt(pstar*pstar+m1*m1),0,+pstar*sin(theta),+pstar*cos(theta)); p4lep2.set(sqrt(pstar*pstar+m2*m2),0,-pstar*sin(theta),-pstar*cos(theta)); //std::cout << "p1: " << p4lep1 << " p2: " << p4lep2 << " pstar: " << pstar << std::endl; if ( i != 100 ) // At maximal q2 we are already in correct frame as Lambda and W/Zvirtual are at rest { p4lep1 = boostTo(p4lep1,boost); p4lep2 = boostTo(p4lep2,boost); } lambda -> init(getDaug(0),p4lambda); lep1 -> init(getDaug(1),p4lep1 ); lep2 -> init(getDaug(2),p4lep2 ); - calcAmp(amp,&parent); + calcAmp(amp,parent); prob = rho.normalizedProb(amp.getSpinDensity()); //std::cout << "q2: " << q2 << " \t theta: " << theta << " \t prob: " << prob << std::endl; //std::cout << "p1: " << p4lep1 << " p2: " << p4lep2 << " q2-q2min: " << q2-(m1+m2)*(m1+m2) << std::endl; if(prob>m_maxProbability){ EvtGenReport(EVTGEN_INFO,"EvtGen") << " - probability " << prob << " found at q2 = " << q2 << " (" << 100*(q2-q2min)/(q2max-q2min) << " %) and theta = " << theta*180/EvtConst::pi << std::endl; m_maxProbability=prob; } } //::abort(); } //m_poleSize = 0.04*q2min; m_maxProbability *= 1.2; + delete parent; } setProbMax(m_maxProbability); EvtGenReport(EVTGEN_INFO,"EvtGen") << " EvtRareLbToLll set up maximum probability to " << m_maxProbability << std::endl; } void EvtRareLbToLll::decay( EvtParticle *parent ){ parent->initializePhaseSpace( getNDaug(),getDaugs() ); calcAmp( _amp2,parent ); } bool EvtRareLbToLll::isParticle( EvtParticle *parent ) const { static EvtIdSet partlist("Lambda_b0"); return partlist.contains( parent->getId() ); } void EvtRareLbToLll::calcAmp( EvtAmp& amp, EvtParticle *parent ){ //parent->setDiagonalSpinDensity(); EvtParticle* lambda = parent->getDaug(0); static EvtIdSet leptons("e-","mu-","tau-"); const bool isparticle = isParticle( parent ); EvtParticle* lp = 0; EvtParticle* lm = 0; if ( leptons.contains(parent->getDaug(1)->getId()) ) { lp = parent->getDaug(1); lm = parent->getDaug(2); } else { lp = parent->getDaug(2); lm = parent->getDaug(1); } EvtVector4R P; P.set(parent->mass(), 0.0,0.0,0.0); EvtVector4R q = lp->getP4() + lm->getP4(); double qsq = q.mass2(); // Leptonic currents EvtVector4C LV[2][2]; // \bar{\ell} \gamma^{\mu} \ell EvtVector4C LA[2][2]; // \bar{\ell} \gamma^{\mu} \gamma^{5} \ell for ( int i =0 ; i < 2; ++i ){ for ( int j = 0; j < 2; ++j ){ if ( isparticle ) { LV[i][j] = EvtLeptonVCurrent( lp->spParent(i), lm->spParent(j) ); LA[i][j] = EvtLeptonACurrent( lp->spParent(i), lm->spParent(j) ); } else { LV[i][j] = EvtLeptonVCurrent( lp->spParent(1-i), lm->spParent(1-j) ); LA[i][j] = EvtLeptonACurrent( lp->spParent(1-i), lm->spParent(1-j) ); } } } EvtRareLbToLllFF::FormFactors FF; //F, G, FT and GT ffmodel_->getFF( parent, lambda, FF ); EvtComplex C7eff = wcmodel_-> GetC7Eff(qsq); EvtComplex C9eff = wcmodel_-> GetC9Eff(qsq); EvtComplex C10eff = wcmodel_->GetC10Eff(qsq); EvtComplex AC[4]; EvtComplex BC[4]; EvtComplex DC[4]; EvtComplex EC[4]; // check to see if particle is same or opposite parity to Lb const int parity = ffmodel_->isNatural( lambda ) ? 1 : -1 ; // Lambda spin type const EvtSpinType::spintype spin = EvtPDL::getSpinType(lambda->getId()); static const double mb = 5.209; // Eq. 48 + 49 for ( unsigned int i = 0; i < 4; ++i ) { if ( parity > 0 ) { AC[i] = -2.*mb*C7eff*FF.FT_[i]/qsq + C9eff*FF.F_[i]; BC[i] = -2.*mb*C7eff*FF.GT_[i]/qsq - C9eff*FF.G_[i]; DC[i] = C10eff*FF.F_[i]; EC[i] = -C10eff*FF.G_[i]; } else { AC[i] = -2.*mb*C7eff*FF.GT_[i]/qsq - C9eff*FF.G_[i]; BC[i] = -2.*mb*C7eff*FF.FT_[i]/qsq + C9eff*FF.F_[i]; DC[i] = -C10eff*FF.G_[i]; EC[i] = C10eff*FF.F_[i]; } } // handle particle -> antiparticle in Hadronic currents const double cv = ( isparticle > 0 ) ? 1.0 : -1.0*parity; const double ca = ( isparticle > 0 ) ? 1.0 : +1.0*parity; const double cs = ( isparticle > 0 ) ? 1.0 : +1.0*parity; const double cp = ( isparticle > 0 ) ? 1.0 : -1.0*parity; if (EvtSpinType::DIRAC == spin ){ EvtVector4C H1[2][2]; // vector current EvtVector4C H2[2][2]; // axial-vector EvtVector4C T[6]; // Hadronic currents for ( int i =0 ; i < 2; ++i ){ for ( int j = 0; j < 2; ++j ){ HadronicAmp( parent, lambda, T, i, j ); H1[i][j] = ( cv*AC[0]*T[0] + ca*BC[0]*T[1] + cs*AC[1]*T[2] + cp*BC[1]*T[3] + cs*AC[2]*T[4] + cp*BC[2]*T[5] ); H2[i][j] = ( cv*DC[0]*T[0] + ca*EC[0]*T[1] + cs*DC[1]*T[2] + cp*EC[1]*T[3] + cs*DC[2]*T[4] + cp*EC[2]*T[5] ); } } // Spin sums int spins[4]; for ( int i =0; i < 2; ++i ) { for ( int ip = 0; ip < 2; ++ip ) { for ( int j = 0; j < 2; ++j ) { for ( int jp = 0; jp < 2; ++jp ) { spins[0] = i; spins[1] = ip; spins[2] = j; spins[3] = jp; EvtComplex M = H1[i][ip]*LV[j][jp] + H2[i][ip]*LA[j][jp]; amp.vertex( spins, M ); } } } } } else if ( EvtSpinType::RARITASCHWINGER == spin ) { EvtVector4C T[8]; EvtVector4C H1[2][4]; // vector current // swaped EvtVector4C H2[2][4]; // axial-vector // Build hadronic amplitude for ( int i =0 ; i < 2; ++i ){ for ( int j = 0; j < 4; ++j ){ H1[i][j] = ( cv*AC[0]*T[0] + ca*BC[0]*T[1] + cs*AC[1]*T[2] + cp*BC[1]*T[3] + cs*AC[2]*T[4] + cp*BC[2]*T[5] + cs*AC[3]*T[6] + cp*BC[3]*T[7] ); H2[i][j] = ( cv*DC[0]*T[0] + ca*EC[0]*T[1] + cs*DC[1]*T[2] + cp*EC[1]*T[3] + cs*DC[2]*T[4] + cp*EC[2]*T[5] + cs*DC[3]*T[6] + cp*EC[3]*T[7] ); } } // Spin sums int spins[4]; for ( int i =0; i < 2; ++i ) { for ( int ip = 0; ip < 4; ++ip ) { for ( int j = 0; j < 2; ++j ) { for ( int jp = 0; jp < 2; ++jp ) { spins[0] = i; spins[1] = ip; spins[2] = j; spins[3] = jp; EvtComplex M = H1[i][ip]*LV[j][jp] + H2[i][ip]*LA[j][jp]; amp.vertex( spins, M ); } } } } } else { EvtGenReport(EVTGEN_ERROR,"EvtGen" ) << " EvtRareLbToLll expects DIRAC or RARITASWINGER daughter " << std::endl; } return; } // spin 1/2 daughters void EvtRareLbToLll::HadronicAmp( EvtParticle* parent, EvtParticle* lambda, EvtVector4C* T, const int i, const int j ) { const EvtDiracSpinor Sfinal = lambda->spParent(j); const EvtDiracSpinor Sinit = parent->sp(i); const EvtVector4R L = lambda->getP4(); EvtVector4R P; P.set(parent->mass(), 0.0,0.0,0.0); const double Pm = parent->mass(); const double Lm = lambda->mass(); // \bar{u} \gamma^{\mu} u T[0] = EvtLeptonVCurrent( Sfinal, Sinit ); // \bar{u} \gamma^{\mu}\gamma^{5} u T[1] = EvtLeptonACurrent( Sfinal, Sinit ); // \bar{u} v^{\mu} u T[2] = EvtLeptonSCurrent( Sfinal, Sinit )*( P/Pm ); // \bar{u} v^{\mu} \gamma^{5} u T[3] = EvtLeptonPCurrent( Sfinal, Sinit )*( P/Pm ); // \bar{u} v^{\prime\mu} u T[4] = EvtLeptonSCurrent( Sfinal, Sinit )*( L/Lm ); // \bar{u} v^{\prime\mu} \gamma^{5} T[5] = EvtLeptonPCurrent( Sfinal, Sinit )*( L/Lm ); // Where: // v = p_{\Lambda_b}/m_{\Lambda_b} // v^{\prime} = p_{\Lambda}/m_{\Lambda} return; } // spin 3/2 daughters void EvtRareLbToLll::HadronicAmpRS( EvtParticle* parent, EvtParticle* lambda, EvtVector4C* T, const int i, const int j ) { const EvtRaritaSchwinger Sfinal = lambda->spRSParent(j); const EvtDiracSpinor Sinit = parent->sp(i); EvtVector4R P; P.set(parent->mass(), 0.0,0.0,0.0); const EvtVector4R L = lambda->getP4(); EvtTensor4C ID; ID.setdiag(1.0,1.0,1.0,1.0); EvtDiracSpinor Sprime; for(int i=0 ; i<4 ; i++ ){ Sprime.set_spinor(i,Sfinal.getVector(i)*P); } const double Pmsq = P.mass2(); const double Pm = parent->mass(); const double PmLm = Pm*lambda->mass(); EvtVector4C V1,V2; for(int i=0;i<4;i++){ V1.set(i,EvtLeptonSCurrent(Sfinal.getSpinor(i),Sinit)); V2.set(i,EvtLeptonPCurrent(Sfinal.getSpinor(i),Sinit)); } // \bar{u}_{alpha} v^{\alpha} \gamma^{\mu} u T[0] = EvtLeptonVCurrent(Sprime, Sinit)*(1/Pm); // \bar{u}_{alpha} v^{\alpha} \gamma^{\mu} \gamma^{5} u T[1] = EvtLeptonACurrent(Sprime, Sinit)*(1/Pm); // \bar{u}_{\alpha} v^{\alpha} v^{\mu} u T[2] = EvtLeptonSCurrent(Sprime, Sinit)*(P/Pmsq); // \bar{u}_{\alpha} v^{\alpha} v^{\mu} \gamma^{5} u T[3] = EvtLeptonPCurrent(Sprime, Sinit)*(P/Pmsq); // \bar{u}_{\alpha} v^{\alpha} v^{\prime \mu} u T[4] = EvtLeptonSCurrent(Sprime, Sinit)*(L/PmLm); // \bar{u}_{\alpha} v^{\alpha} v^{\prime \mu} \gamma^{5} u T[5] = EvtLeptonPCurrent(Sprime, Sinit)*(L/PmLm); // \bar{u}_{\alpha} g^{\alpha\mu} u T[6] = ID.cont2(V1); // \bar{u}_{\alpha} g^{\alpha\mu} \gamma^{5} u T[7] = ID.cont2(V2); // Where: // v = p_{\Lambda_b}/m_{\Lambda_b} // v^{\prime} = p_{\Lambda}/m_{\Lambda} return; } diff --git a/src/EvtGenModels/EvtVubNLO.cpp b/src/EvtGenModels/EvtVubNLO.cpp index 46151dc..d1ef4b9 100644 --- a/src/EvtGenModels/EvtVubNLO.cpp +++ b/src/EvtGenModels/EvtVubNLO.cpp @@ -1,642 +1,643 @@ //-------------------------------------------------------------------------- // // 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: // Copyright (C) 1998 Caltech, UCSB // // Module: EvtVubNLO.cc // // Description: Routine to decay B->Xulnu according to Bosch, Lange, Neubert, and Paz hep-ph/0402094 // Equation numbers refer to this paper // // Modification history: // // Riccardo Faccini Feb. 11, 2004 // //------------------------------------------------------------------------ // #include "EvtGenBase/EvtPatches.hh" #include #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtGenKine.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenModels/EvtVubNLO.hh" #include #include "EvtGenBase/EvtVector4R.hh" #include "EvtGenModels/EvtItgSimpsonIntegrator.hh" #include "EvtGenModels/EvtBtoXsgammaFermiUtil.hh" #include "EvtGenModels/EvtItgPtrFunction.hh" #include "EvtGenModels/EvtPFermi.hh" #include "EvtGenBase/EvtRandom.hh" #include "EvtGenBase/EvtDiLog.hh" +#include using std::cout; using std::endl; EvtVubNLO::~EvtVubNLO() { cout <<" max pdf : "<<_gmax< sCoeffs(11); sCoeffs[3] = _b; sCoeffs[4] = _mb; sCoeffs[5] = _mB; sCoeffs[6] = _idSF; sCoeffs[7] = lambda_SF(); sCoeffs[8] = mu_h(); sCoeffs[9] = mu_i(); sCoeffs[10] = 1.; _SFNorm = SFNorm(sCoeffs) ; // SF normalization; cout << " pdf 0.66, 1.32 , 4.32 "<0 && _masses[i] <= _masses[i-1]) { EvtGenReport(EVTGEN_ERROR,"EvtGen") << "EvtVubNLO generator expected " << " mass bins in ascending order!" << "Will terminate execution!"<= 0, but found: " <<_weights[i] < maxw ) maxw = _weights[i]; } if (maxw == 0) { EvtGenReport(EVTGEN_ERROR,"EvtGen") << "EvtVubNLO generator expected at least one " << " weight > 0, but found none! " << "Will terminate execution!"< u-bar specflav l+ nu EvtParticle *xuhad, *lepton, *neutrino; EvtVector4R p4; double pp,pm,pl,ml,El(0.0),Eh(0.0),sh(0.0); p->initializePhaseSpace(getNDaug(),getDaugs()); xuhad=p->getDaug(0); lepton=p->getDaug(1); neutrino=p->getDaug(2); _mB = p->mass(); ml = lepton->mass(); bool tryit = true; while (tryit) { // pm=(E_H+P_H) pm= EvtRandom::Flat(0.,1); pm= pow(pm,1./3.)*_mB; // pl=mB-2*El pl = EvtRandom::Flat(0.,1); pl=sqrt(pl)*pm; // pp=(E_H-P_H) pp = EvtRandom::Flat(0.,pl); _ntot++; El = (_mB-pl)/2.; Eh = (pp+pm)/2; sh = pp*pm; double pdf(0.); if (ppml&& sh > _masses[0]*_masses[0]&& _mB*_mB + sh - 2*_mB*Eh > ml*ml) { double xran = EvtRandom::Flat(0,_dGMax); pdf = tripleDiff(pp,pl,pm); // triple differential distribution // cout <<" P+,P-,Pl,Pdf= "<_dGMax){ EvtGenReport(EVTGEN_ERROR,"EvtGen") << "EvtVubNLO pdf above maximum: " <= xran ) tryit = false; if(pdf>_gmax)_gmax=pdf; } else { // cout <<" EvtVubNLO incorrect kinematics sh= "< _masses[j] ) j++; double w = _weights[j-1]; if ( w < xran1 ) tryit = true;// through away this candidate } } // cout <<" max prob "<init( getDaug(0), p4); // calculate the W 4 vector in the B Meson restrframe double apWB = ptmp; double pWB[4] = {_mB-Eh,-pHB[1],-pHB[2],-pHB[3]}; // first go in the W restframe and calculate the lepton and // the neutrino in the W frame double mW2 = _mB*_mB + sh - 2*_mB*Eh; // if(mW2<0.1){ // cout <<" low Q2! "< 1 ) ctL = 1; sttmp = sqrt(1-ctL*ctL); // eX' = eZ x eW double xW[3] = {-pWB[2],pWB[1],0}; // eZ' = eW double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB}; double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]); for (int j=0;j<2;j++) xW[j] /= lx; // eY' = eZ' x eX' double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]}; double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]); for (int j=0;j<3;j++) yW[j] /= ly; // p_lep = |p_lep| * ( sin(Theta) * cos(Phi) * eX' // + sin(Theta) * sin(Phi) * eY' // + cos(Theta) * eZ') for (int j=0;j<3;j++) pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j] + sttmp*sin(phL)*ptmp*yW[j] + ctL *ptmp*zW[j]; double apLW = ptmp; // boost them back in the B Meson restframe double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW; ptmp = sqrt(El*El-ml*ml); double ctLL = appLB/ptmp; if ( ctLL > 1 ) ctLL = 1; if ( ctLL < -1 ) ctLL = -1; double pLB[4] = {El,0,0,0}; double pNB[8] = {pWB[0]-El,0,0,0}; for (int j=1;j<4;j++) { pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j]; pNB[j] = pWB[j] - pLB[j]; } p4.set(pLB[0],pLB[1],pLB[2],pLB[3]); lepton->init( getDaug(1), p4); p4.set(pNB[0],pNB[1],pNB[2],pNB[3]); neutrino->init( getDaug(2), p4); return ; } double EvtVubNLO::tripleDiff ( double pp, double pl, double pm){ std::vector sCoeffs(11); sCoeffs[0] = pp; sCoeffs[1] = pl; sCoeffs[2] = pm; sCoeffs[3] = _b; sCoeffs[4] = _mb; sCoeffs[5] = _mB; sCoeffs[6] = _idSF; sCoeffs[7] = lambda_SF(); sCoeffs[8] = mu_h(); sCoeffs[9] = mu_i(); sCoeffs[10] = _SFNorm; // SF normalization; double c1=(_mB+pl-pp-pm)*(pm-pl); double c2=2*(pl-pp)*(pm-pl); double c3=(_mB-pm)*(pm-pp); double aF1=F10(sCoeffs); double aF2=F20(sCoeffs); double aF3=F30(sCoeffs); double td0=c1*aF1+c2*aF2+c3*aF3; auto func = EvtItgPtrFunction{&integrand, 0., _mB, sCoeffs}; auto jetSF = EvtItgSimpsonIntegrator{func,0.01,25}; double smallfrac=0.000001;// stop a bit before the end to avoid problems with numerical integration double tdInt = jetSF.evaluate(0,pp*(1-smallfrac)); double SU=U1lo(mu_h(),mu_i())*pow((pm-pp)/(_mB-pp),alo(mu_h(),mu_i())); double TD=(_mB-pp)*SU*(td0+tdInt); return TD; } double EvtVubNLO::integrand(double omega, const std::vector &coeffs){ //double pp=coeffs[0]; double c1=(coeffs[5]+coeffs[1]-coeffs[0]-coeffs[2])*(coeffs[2]-coeffs[1]); double c2=2*(coeffs[1]-coeffs[0])*(coeffs[2]-coeffs[1]); double c3=(coeffs[5]-coeffs[2])*(coeffs[2]-coeffs[0]); return c1*F1Int(omega,coeffs)+c2*F2Int(omega,coeffs)+c3*F3Int(omega,coeffs); } double EvtVubNLO::F10(const std::vector &coeffs){ double pp=coeffs[0]; double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); double mui=coeffs[9]; double muh=coeffs[8]; double z=1-y; double result= U1nlo(muh,mui)/ U1lo(muh,mui); result += anlo(muh,mui)*log(y); result += C_F(muh)*(-4*pow(log(y*coeffs[4]/muh),2)+10*log(y*coeffs[4]/muh)-4*log(y)-2*log(y)/(1-y)-4.0*EvtDiLog::DiLog(z)-pow(EvtConst::pi,2)/6.-12 ); result += C_F(mui)*(2*pow(log(y*coeffs[4]*pp/pow(mui,2)),2)-3*log(y*coeffs[4]*pp/pow(mui,2))+7-pow(EvtConst::pi,2) ); result *=shapeFunction(pp,coeffs); // changes due to SSF result += (-subS(coeffs)+2*subT(coeffs)+(subU(coeffs)-subV(coeffs))*(1/y-1.))/(coeffs[5]-pp); result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0]),2)*(-5*(lambda1()+3*lambda2())/6+2*(2*lambda1()/3-lambda2())/pow(y,2)); // result += (subS(coeffs)+subT(coeffs)+(subU(coeffs)-subV(coeffs))/y)/(coeffs[5]-pp); // this part has been added after Feb '05 //result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0]),2)*((lambda1()+3*lambda2())/6+2*(2*lambda1()/3-lambda2())/pow(y,2)); return result; } double EvtVubNLO::F1Int(double omega,const std::vector &coeffs){ double pp=coeffs[0]; double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); // mubar == mui return C_F(coeffs[9])*( (shapeFunction(omega,coeffs)-shapeFunction(pp,coeffs))*(4*log(y*coeffs[4]*(pp-omega)/pow(coeffs[9],2))-3)/(pp-omega)+ (g1(y,(pp-omega)/(coeffs[5]-coeffs[0]))/(coeffs[5]-pp)*shapeFunction(omega,coeffs)) ); } double EvtVubNLO::F20(const std::vector &coeffs){ double pp=coeffs[0]; double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); double result= C_F(coeffs[8])*log(y)/(1-y)*shapeFunction(pp,coeffs)- 1/y*(subS(coeffs)+2*subT(coeffs)-(subT(coeffs)+subV(coeffs))/y)/(coeffs[5]-pp); // added after Feb '05 result += shapeFunction(pp,coeffs)/pow((coeffs[5]-coeffs[0])*y,2)*(2*lambda1()/3+4*lambda2()-y*(7/6*lambda1()+3*lambda2())); return result; } double EvtVubNLO::F2Int(double omega,const std::vector &coeffs){ double pp=coeffs[0]; double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); return C_F(coeffs[9])*g3(y,(pp-omega)/(coeffs[5]-coeffs[0]))*shapeFunction(omega,coeffs)/(coeffs[5]-pp); } double EvtVubNLO::F30(const std::vector &coeffs){ double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); return shapeFunction(coeffs[0],coeffs)/pow((coeffs[5]-coeffs[0])*y,2)*(-2*lambda1()/3+lambda2()); } double EvtVubNLO::F3Int(double omega,const std::vector &coeffs){ double pp=coeffs[0]; double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); return C_F(coeffs[9])*g3(y,(pp-omega)/(coeffs[5]-coeffs[0]))/2*shapeFunction(omega,coeffs)/(coeffs[2]-coeffs[0]); } double EvtVubNLO::g1(double y, double x){ double result=(y*(-9+10*y)+x*x*(-12+13*y)+2*x*(-8+6*y+3*y*y))/y/pow(1+x,2)/(x+y); result -= 4*log((1+1/x)*y)/x; result -=2*log(1+y/x)*(3*pow(x,4)*(-2+y)-2*pow(y,3)-4*pow(x,3)*(2+y)-2*x*y*y*(4+y)-x*x*y*(12+4*y+y*y))/x/pow((1+x)*y,2)/(x+y); return result; } double EvtVubNLO::g2(double y, double x){ double result=y*(10*pow(x,4)+y*y+3*x*x*y*(10+y)+pow(x,3)*(12+19*y)+x*y*(8+4*y+y*y)); result -= 2*x*log(1+y/x)*(5*pow(x,4)+2*y*y+6*pow(x,3)*(1+2*y)+4*y*x*(1+2*y)+x*x*y*(18+5*y)); result *= 2/(pow(y*(1+x),2)*y*(x+y)); return result; } double EvtVubNLO::g3(double y, double x){ double result=(2*pow(y,3)*(-11+2*y)-10*pow(x,4)*(6-6*y+y*y)+x*y*y*(-94+29*y+2*y*y)+2*x*x*y*(-72+18*y+13*y*y)-pow(x,3)*(72+42*y-70*y*y+3*pow(y,3)))/(pow(y*(1+x),2)*y*(x+y)); result += 2*log(1+y/x)*(-6*x*pow(y,3)*(-5+y)+4*pow(y,4)+5*pow(x,5)*(6-6*y+y*y)-4*pow(x*y,2)*(-20+6*y+y*y)+pow(x,3)*y*(90-10*y-28*y*y+pow(y,3))+pow(x,4)*(36+36*y-50*y*y+4*pow(y,3)))/(pow((1+x)*y*y,2)*(x+y)); return result; } /* old version (before Feb 05 notebook from NNeubert double EvtVubNLO::F1Int(double omega,const std::vector &coeffs){ double pp=coeffs[0]; double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); // mubar == mui return C_F(coeffs[9])*( (shapeFunction(omega,coeffs)-shapeFunction(pp,coeffs))*(4*log(y*coeffs[4]*(pp-omega)/pow(coeffs[9],2))-3)/(pp-omega)- (1./y/(coeffs[5]-pp)*shapeFunction(omega,coeffs)*(5-6*y+4*(3-y)*log((pp-omega)/y/coeffs[4]))) ); } double EvtVubNLO::F2Int(double omega,const std::vector &coeffs){ double pp=coeffs[0]; double y=(coeffs[2]-coeffs[0])/(coeffs[5]-coeffs[0]); return C_F(coeffs[9])*shapeFunction(omega,coeffs)*(2-11/y-4/y*log((pp-omega)/y/coeffs[4]))/(coeffs[5]-pp); } double EvtVubNLO::F3(const std::vector &coeffs){ return C_F(coeffs[9])*shapeFunction(omega,coeffs)/(coeffs[2]-coeffs[0]); } */ double EvtVubNLO::SFNorm( const std::vector &/*coeffs*/){ double omega0=1.68;//normalization scale (mB-2*1.8) if(_idSF==1){ // exponential SF double omega0=1.68;//normalization scale (mB-2*1.8) return M0(mu_i(),omega0)*pow(_b,_b)/lambda_SF()/ (Gamma(_b)-Gamma(_b,_b*omega0/lambda_SF())); } else if(_idSF==2){ // Gaussian SF double c=cGaus(_b); return M0(mu_i(),omega0)*2/lambda_SF()/pow(c,-(1+_b)/2.)/ (Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c)); } else { EvtGenReport(EVTGEN_ERROR,"EvtGen") << "unknown SF "<<_idSF< &sCoeffs){ if( sCoeffs[6]==1){ return sCoeffs[10]*expShapeFunction(omega, sCoeffs); } else if( sCoeffs[6]==2) { return sCoeffs[10]*gausShapeFunction(omega, sCoeffs); } else { EvtGenReport(EVTGEN_ERROR,"EvtGen") << "EvtVubNLO : unknown shape function # " < &c){ return (lambda_bar(1.68)-c[0])*shapeFunction(c[0],c);} double EvtVubNLO::subT(const std::vector &c){ return -3*lambda2()*subS(c)/mu_pi2(1.68);} double EvtVubNLO::subU(const std::vector &c){ return -2*subS(c);} double EvtVubNLO::subV(const std::vector &c){ return -subT(c);} double EvtVubNLO::lambda_bar(double omega0){ if(_lbar<0){ if(_idSF==1){ // exponential SF double rat=omega0*_b/lambda_SF(); _lbar=lambda_SF()/_b*(Gamma(1+_b)-Gamma(1+_b,rat))/(Gamma(_b)-Gamma(_b,rat)); } else if(_idSF==2){ // Gaussian SF double c=cGaus(_b); _lbar=lambda_SF()*(Gamma(1+_b/2)-Gamma(1+_b/2,pow(omega0/lambda_SF(),2)*c))/(Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c))/sqrt(c); } } return _lbar; } double EvtVubNLO::mu_pi2(double omega0){ if(_mupi2<0){ if(_idSF==1){ // exponential SF double rat=omega0*_b/lambda_SF(); _mupi2= 3*(pow(lambda_SF()/_b,2)*(Gamma(2+_b)-Gamma(2+_b,rat))/(Gamma(_b)-Gamma(_b,rat))-pow(lambda_bar(omega0),2)); } else if(_idSF==2){ // Gaussian SF double c=cGaus(_b); double m1=Gamma((3+_b)/2)-Gamma((3+_b)/2,pow(omega0/lambda_SF(),2)*c); double m2=Gamma(1+_b/2)-Gamma(1+_b/2,pow(omega0/lambda_SF(),2)*c); double m3=Gamma((1+_b)/2)-Gamma((1+_b)/2,pow(omega0/lambda_SF(),2)*c); _mupi2= 3*pow(lambda_SF(),2)*(m1/m3-pow(m2/m3,2))/c; } } return _mupi2; } double EvtVubNLO::M0(double mui,double omega0){ double mf=omega0-lambda_bar(omega0); return 1+4*C_F(mui)*(-pow(log(mf/mui),2)-log(mf/mui)-pow(EvtConst::pi/2,2)/6.+mu_pi2(omega0)/3/pow(mf,2)*(log(mf/mui)-0.5)); } double EvtVubNLO::alphas(double mu){ double Lambda4=0.302932; double lg=2*log(mu/Lambda4); return 4*EvtConst::pi/lg/beta0()*(1-beta1()*log(lg)/pow(beta0(),2)/lg+pow(beta1()/lg,2)/pow(beta0(),4)*(pow(log(lg)-0.5,2)-1.25+beta2()*beta0()/pow(beta1(),2))); } double EvtVubNLO::gausShapeFunction ( double omega, const std::vector &sCoeffs){ double b=sCoeffs[3]; double l=sCoeffs[7]; double wL=omega/l; return pow(wL,b)*exp(-cGaus(b)*wL*wL); } double EvtVubNLO::expShapeFunction ( double omega, const std::vector &sCoeffs){ double b=sCoeffs[3]; double l=sCoeffs[7]; double wL=omega/l; return pow(wL,b-1)*exp(-b*wL); } double EvtVubNLO::Gamma(double z) { - auto gammaCoeffs = std::array{ 76.18009172947146, -86.50532032941677, - 24.01409824083091, -1.231739572450155, - 0.1208650973866179e-2, -0.5395239384953e-5 }; + auto gammaCoeffs = std::array{ 76.18009172947146, -86.50532032941677, + 24.01409824083091, -1.231739572450155, + 0.1208650973866179e-2, -0.5395239384953e-5 }; //Lifted from Numerical Recipies in C double y = z; double x = z; double tmp = x + 5.5; tmp = tmp - (x+0.5)*log(tmp); double ser=1.000000000190015; for (const auto& gammaCoeff : gammaCoeffs) { y += 1.0; ser += gammaCoeff/y; } return exp(-tmp+log(2.5066282746310005*ser/x)); } double EvtVubNLO::Gamma(double z, double tmin) { std::vector c(1); c[0]=z; auto func = EvtItgPtrFunction{&dgamma, tmin, 100., c}; auto jetSF = EvtItgSimpsonIntegrator{func,0.001}; return jetSF.evaluate(tmin,100.); } diff --git a/src/EvtGenModels/EvtbTosllMSFF.cpp b/src/EvtGenModels/EvtbTosllMSFF.cpp index 8bd67f7..97b3eab 100644 --- a/src/EvtGenModels/EvtbTosllMSFF.cpp +++ b/src/EvtGenModels/EvtbTosllMSFF.cpp @@ -1,626 +1,621 @@ //-------------------------------------------------------------------------- // // 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) 2000 Caltech, UCSB // // Module: EvtbTosllMSFF.cpp // Description: Form factors for B -> (P, V) ell^+ ell^- transitions according // to the paper: D.Melikhov, B.Stech, PRD62, 014006 (2000). // // Modification history: // // N.Nikitin (nnikit@mail.cern.ch) March 13, 2008 Module created // N.Nikitin (nnikit@mail.cern.ch) March 27, 2008 add \bar B -> \bar (K,K^*) transition ff // N.Nikitin (nnikit@mail.cern.ch) April 26, 2008 add \bar Bs -> phi transition ff // N.Nikitin (nnikit@mail.cern.ch) April 26, 2008 add \bar Bs -> K* transition ff // N.Nikitin (nnikit@mail.cern.ch) April 27, 2008 add \bar B -> \bar rho transition ff // N.Nikitin (nnikit@mail.cern.ch) Nvmbr 04, 2011 add \bar B -> omega transition ff // N.Nikitin (nnikit@mail.cern.ch) Dec 16, 2011 add \bar B -> \bar K_1(1270) transition ff (from H.Hatanaka and Kwei-Chou Yang, PRD78, 074007 (2008)) // N.Nikitin (nnikit@mail.cern.ch) Dec 16, 2011 add \bar B -> \bar K_1(1400) transition ff (from H.Hatanaka and Kwei-Chou Yang, PRD78, 074007 (2008)) // N.Nikitin (nnikit@mail.cern.ch) May 11, 2012 add \bar Bs -> f_0(980) transition ff with NLO corrections from Table II (see P.Colangelo et al., PRD81, 074001 (2010)) // //------------------------------------------------------------------------ #include "EvtGenBase/EvtPatches.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtComplex.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenModels/EvtbTosllMSFF.hh" #include #include EvtbTosllMSFF::EvtbTosllMSFF(){} double EvtbTosllMSFF::equation9_10(double ff0, double M2, double q2, - double sigma1, double sigma2, int eq_num) - { - double ff; - - ff=1.0; + double sigma1, double sigma2, int eq_num) { + double ff=1.0; - switch(eq_num) - { - case 9: - ff=1./(1.-q2/M2); - [[fallthrough]]; - case 10: - ff=ff*ff0/(1.-sigma1*q2/M2+sigma2*pow(q2,2)/pow(M2,2)); - break; - default: - EvtGenReport(EVTGEN_ERROR,"EvtGen") << - "In the function EvtbTosllMSFF::equation9_10 \n" << - "the parameter eq_num non equal to the 9 or 10! \n" << - "eq_num =" << eq_num < K transition form factors" // << std::endl; } // B -> \pi transition form factors if((parent == EvtPDL::getId(std::string("B+"))&& daught == EvtPDL::getId(std::string("pi+")))|| (parent == EvtPDL::getId(std::string("B-"))&& daught == EvtPDL::getId(std::string("pi-")))|| (parent == EvtPDL::getId(std::string("B0"))&& daught == EvtPDL::getId(std::string("pi0")))|| (parent == EvtPDL::getId(std::string("anti-B0"))&& daught == EvtPDL::getId(std::string("pi0")))){ double ff0[] ={0.29, 0.29, 0.28}; double sigma1[]={0.48, 0.76, 0.48}; double sigma2[]={0.00, 0.28, 0.00}; int eq_num[]={ 9, 10, 9}; double M_P2=5.27*5.27; // GeV^2 for B^0 - meson double M_V2=5.32*5.32; // GeV^2 for B^* - meson fp = equation9_10(ff0[0], M_P2, t, sigma1[0], sigma2[0], eq_num[0]); f0 = equation9_10(ff0[1], M_V2, t, sigma1[1], sigma2[1], eq_num[1]); ft = equation9_10(ff0[2], M_P2, t, sigma1[2], sigma2[2], eq_num[2]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n B -> pi transition form factors" // << std::endl; } // B_d -> \eta transition form factors if((parent == EvtPDL::getId(std::string("B0"))&& daught == EvtPDL::getId(std::string("eta")))|| (parent == EvtPDL::getId(std::string("anti-B0"))&& daught == EvtPDL::getId(std::string("eta")))){ double ff0[] ={0.36, 0.36, 0.36}; double sigma1[]={0.60, 0.80, 0.58}; double sigma2[]={0.20, 0.40, 0.18}; int eq_num[]={ 9, 10, 9}; double M_P2=5.27*5.27; // GeV^2 for B_d^0 - meson double M_V2=5.32*5.32; // GeV^2 for B_d^* - meson fp = equation9_10(ff0[0], M_P2, t, sigma1[0], sigma2[0], eq_num[0]); fp = -0.5*fp; f0 = equation9_10(ff0[1], M_V2, t, sigma1[1], sigma2[1], eq_num[1]); f0 = -0.5*f0; ft = equation9_10(ff0[2], M_P2, t, sigma1[2], sigma2[2], eq_num[2]); ft = -0.5*ft; models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n Bd -> eta transition form factors" // << std::endl; } // B_d -> \eta' transition form factors if((parent == EvtPDL::getId(std::string("B0"))&& daught == EvtPDL::getId(std::string("eta'")))|| (parent == EvtPDL::getId(std::string("anti-B0"))&& daught == EvtPDL::getId(std::string("eta'")))){ double ff0[] ={0.36, 0.36, 0.39}; double sigma1[]={0.60, 0.80, 0.58}; double sigma2[]={0.20, 0.45, 0.18}; int eq_num[]={ 9, 10, 9}; double M_P2=5.27*5.27; // GeV^2 for B_d^0 - meson double M_V2=5.32*5.32; // GeV^2 for B_d^* - meson fp = equation9_10(ff0[0], M_P2, t, sigma1[0], sigma2[0], eq_num[0]); f0 = equation9_10(ff0[1], M_V2, t, sigma1[1], sigma2[1], eq_num[1]); ft = equation9_10(ff0[2], M_P2, t, sigma1[2], sigma2[2], eq_num[2]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n Bd -> eta transition form factors" // << std::endl; } // B_s -> \eta transition form factors if((parent == EvtPDL::getId(std::string("B_s0"))&& daught == EvtPDL::getId(std::string("eta")))|| (parent == EvtPDL::getId(std::string("anti-B_s0"))&& daught == EvtPDL::getId(std::string("eta")))){ double ff0[] ={0.36, 0.36, 0.36}; double sigma1[]={0.60, 0.80, 0.58}; double sigma2[]={0.20, 0.40, 0.18}; int eq_num[]={ 9, 10, 9}; double M_P2=5.37*5.37; // GeV^2 for B_s^0 - meson double M_V2=5.42*5.42; // GeV^2 for B_s^* - meson fp = equation9_10(ff0[0], M_P2, t, sigma1[0], sigma2[0], eq_num[0]); f0 = equation9_10(ff0[1], M_V2, t, sigma1[1], sigma2[1], eq_num[1]); ft = equation9_10(ff0[2], M_P2, t, sigma1[2], sigma2[2], eq_num[2]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n Bs -> eta transition form factors" // << std::endl; } // B_s -> \eta' transition form factors if((parent == EvtPDL::getId(std::string("B_s0"))&& daught == EvtPDL::getId(std::string("eta'")))|| (parent == EvtPDL::getId(std::string("anti-B_s0"))&& daught == EvtPDL::getId(std::string("eta'")))){ double ff0[] ={0.36, 0.36, 0.39}; double sigma1[]={0.60, 0.80, 0.58}; double sigma2[]={0.20, 0.45, 0.18}; int eq_num[]={ 9, 10, 9}; double M_P2=5.37*5.37; // GeV^2 for B_s^0 - meson double M_V2=5.42*5.42; // GeV^2 for B_s^* - meson fp = equation9_10(ff0[0], M_P2, t, sigma1[0], sigma2[0], eq_num[0]); f0 = equation9_10(ff0[1], M_V2, t, sigma1[1], sigma2[1], eq_num[1]); ft = equation9_10(ff0[2], M_P2, t, sigma1[2], sigma2[2], eq_num[2]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n Bs -> eta transition form factors" // << std::endl; } // B_s -> f_0(980) transition form factors if((parent == EvtPDL::getId(std::string("B_s0"))&& daught == EvtPDL::getId(std::string("f_0")))|| (parent == EvtPDL::getId(std::string("anti-B_s0"))&& daught == EvtPDL::getId(std::string("f_0")))){ double ff0[] ={0.238, 0.238, 0.308}; double sigma1[]={ 1.50, 0.53, 1.46}; double sigma2[]={ 0.58, -0.36, 0.58}; int eq_num[]={ 10, 10, 10}; double M_P2=5.366*5.366; // GeV^2 for B_s^0 - meson fp = 0.0 - equation9_10(ff0[0], M_P2, t, sigma1[0], sigma2[0], eq_num[0]); f0 = 0.0 - equation9_10(ff0[1], M_P2, t, sigma1[1], sigma2[1], eq_num[1]); ft = equation9_10(ff0[2], M_P2, t, sigma1[2], sigma2[2], eq_num[2]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n Bs -> eta transition form factors" // << std::endl; } // EvtGenReport(EVTGEN_NOTICE,"EvtGen") << "\n models_counter = " << models_counter // << "\n Scalar form-factors at q^2 = " << t // << " for B -> P transition:" // << "\n fp = " << fp // << "\n f0 = " << f0 // << "\n ft = " << ft << std::endl; if(models_counter!=1){ EvtGenReport(EVTGEN_ERROR,"EvtGen") << "\n In the function EvtbTosllMSFF::getScalarFF(...) \n" << "the parameter models_counter not equal 1! \n" << "models_counter = " << models_counter < \bar K* transition form factors if((parent == EvtPDL::getId(std::string("B+"))&& daught == EvtPDL::getId(std::string("K*+")))|| (parent == EvtPDL::getId(std::string("B-"))&& daught == EvtPDL::getId(std::string("K*-")))|| (parent == EvtPDL::getId(std::string("B0"))&& daught == EvtPDL::getId(std::string("K*0")))|| (parent == EvtPDL::getId(std::string("anti-B0"))&& daught == EvtPDL::getId(std::string("anti-K*0")))){ double ff0[] ={0.44, 0.45, 0.36, 0.32, 0.39, 0.39, 0.27}; double sigma1[]={0.45, 0.46, 0.64, 1.23, 0.45, 0.72, 1.31}; double sigma2[]={0.00, 0.00, 0.36, 0.38, 0.00, 0.62, 0.41}; int eq_num[]={ 9, 9, 10, 10, 9, 10, 10}; double M_P2=5.37*5.37; // GeV^2 for B^0_s - meson double M_V2=5.42*5.42; // GeV^2 for B^*_s - meson v =equation9_10(ff0[0], M_P2, t, sigma1[0], sigma2[0], eq_num[0]); a0=equation9_10(ff0[1], M_P2, t, sigma1[1], sigma2[1], eq_num[1]); a1=equation9_10(ff0[2], M_V2, t, sigma1[2], sigma2[2], eq_num[2]); a2=equation9_10(ff0[3], M_V2, t, sigma1[3], sigma2[3], eq_num[3]); t1=equation9_10(ff0[4], M_P2, t, sigma1[4], sigma2[4], eq_num[4]); t2=equation9_10(ff0[5], M_V2, t, sigma1[5], sigma2[5], eq_num[5]); t3=equation9_10(ff0[6], M_V2, t, sigma1[6], sigma2[6], eq_num[6]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n barB -> barK* transition form factors" // << std::endl; } // \bar B -> \bar\rho transition form factors if((parent == EvtPDL::getId(std::string("B+"))&& daught == EvtPDL::getId(std::string("rho+")))|| (parent == EvtPDL::getId(std::string("B-"))&& daught == EvtPDL::getId(std::string("rho-")))|| (parent == EvtPDL::getId(std::string("B0"))&& daught == EvtPDL::getId(std::string("rho0")))|| (parent == EvtPDL::getId(std::string("anti-B0"))&& daught == EvtPDL::getId(std::string("rho0")))){ double ff0[] ={0.31, 0.30, 0.26, 0.24, 0.27, 0.27, 0.19}; double sigma1[]={0.59, 0.54, 0.73, 1.40, 0.60, 0.74, 1.42}; double sigma2[]={0.00, 0.00, 0.10, 0.50, 0.00, 0.19, 0.51}; int eq_num[]={ 9, 9, 10, 10, 9, 10, 10}; double M_P2=5.27*5.27; // GeV^2 for B - meson double M_V2=5.32*5.32; // GeV^2 for B^* - meson v =equation9_10(ff0[0], M_P2, t, sigma1[0], sigma2[0], eq_num[0]); a0=equation9_10(ff0[1], M_P2, t, sigma1[1], sigma2[1], eq_num[1]); a1=equation9_10(ff0[2], M_V2, t, sigma1[2], sigma2[2], eq_num[2]); a2=equation9_10(ff0[3], M_V2, t, sigma1[3], sigma2[3], eq_num[3]); t1=equation9_10(ff0[4], M_P2, t, sigma1[4], sigma2[4], eq_num[4]); t2=equation9_10(ff0[5], M_V2, t, sigma1[5], sigma2[5], eq_num[5]); t3=equation9_10(ff0[6], M_V2, t, sigma1[6], sigma2[6], eq_num[6]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n barB -> bar rho transition form factors" // << std::endl; } // \bar B -> \omega transition form factors (exactly as for \bar B -> \rho^0 ff!) if((parent == EvtPDL::getId(std::string("B0"))&& daught == EvtPDL::getId(std::string("omega")))|| (parent == EvtPDL::getId(std::string("anti-B0"))&& daught == EvtPDL::getId(std::string("omega")))){ double ff0[] ={0.31, 0.30, 0.26, 0.24, 0.27, 0.27, 0.19}; double sigma1[]={0.59, 0.54, 0.73, 1.40, 0.60, 0.74, 1.42}; double sigma2[]={0.00, 0.00, 0.10, 0.50, 0.00, 0.19, 0.51}; int eq_num[]={ 9, 9, 10, 10, 9, 10, 10}; double M_P2=5.27*5.27; // GeV^2 for B - meson double M_V2=5.32*5.32; // GeV^2 for B^* - meson v =equation9_10(ff0[0], M_P2, t, sigma1[0], sigma2[0], eq_num[0]); a0=equation9_10(ff0[1], M_P2, t, sigma1[1], sigma2[1], eq_num[1]); a1=equation9_10(ff0[2], M_V2, t, sigma1[2], sigma2[2], eq_num[2]); a2=equation9_10(ff0[3], M_V2, t, sigma1[3], sigma2[3], eq_num[3]); t1=equation9_10(ff0[4], M_P2, t, sigma1[4], sigma2[4], eq_num[4]); t2=equation9_10(ff0[5], M_V2, t, sigma1[5], sigma2[5], eq_num[5]); t3=equation9_10(ff0[6], M_V2, t, sigma1[6], sigma2[6], eq_num[6]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n barB -> omega transition form factors" // << std::endl; } // \bar Bs -> phi transition form factors if((parent == EvtPDL::getId(std::string("B_s0"))&& daught == EvtPDL::getId(std::string("phi")))|| (parent == EvtPDL::getId(std::string("anti-B_s0"))&& daught == EvtPDL::getId(std::string("phi")))){ double ff0[] ={0.44, 0.42, 0.34, 0.31, 0.38, 0.38, 0.26}; double sigma1[]={0.62, 0.55, 0.73, 1.30, 0.62, 0.83, 1.41}; double sigma2[]={0.20, 0.12, 0.42, 0.52, 0.20, 0.71, 0.57}; int eq_num[]={ 9, 9, 10, 10, 9, 10, 10}; double M_P2=5.37*5.37; // GeV^2 for B^0_s - meson double M_V2=5.42*5.42; // GeV^2 for B^*_s - meson v =equation9_10(ff0[0], M_P2, t, sigma1[0], sigma2[0], eq_num[0]); a0=equation9_10(ff0[1], M_P2, t, sigma1[1], sigma2[1], eq_num[1]); a1=equation9_10(ff0[2], M_V2, t, sigma1[2], sigma2[2], eq_num[2]); a2=equation9_10(ff0[3], M_V2, t, sigma1[3], sigma2[3], eq_num[3]); t1=equation9_10(ff0[4], M_P2, t, sigma1[4], sigma2[4], eq_num[4]); t2=equation9_10(ff0[5], M_V2, t, sigma1[5], sigma2[5], eq_num[5]); t3=equation9_10(ff0[6], M_V2, t, sigma1[6], sigma2[6], eq_num[6]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n barBs -> phi transition form factors" // << std::endl; } // \bar Bs -> K* (without \bar !) transition form factors if((parent == EvtPDL::getId(std::string("B_s0"))&& daught == EvtPDL::getId(std::string("anti-K*0")))|| (parent == EvtPDL::getId(std::string("anti-B_s0"))&& daught == EvtPDL::getId(std::string("K*0")))){ double ff0[] ={0.38, 0.37, 0.29, 0.26, 0.32, 0.32, 0.23}; double sigma1[]={0.66, 0.60, 0.86, 1.32, 0.66, 0.98, 1.42}; double sigma2[]={0.30, 0.16, 0.60, 0.54, 0.31, 0.90, 0.62}; int eq_num[]={ 9, 9, 10, 10, 9, 10, 10}; double M_P2=5.27*5.27; // GeV^2 for B - meson double M_V2=5.32*5.32; // GeV^2 for B^* - meson v =equation9_10(ff0[0], M_P2, t, sigma1[0], sigma2[0], eq_num[0]); a0=equation9_10(ff0[1], M_P2, t, sigma1[1], sigma2[1], eq_num[1]); a1=equation9_10(ff0[2], M_V2, t, sigma1[2], sigma2[2], eq_num[2]); a2=equation9_10(ff0[3], M_V2, t, sigma1[3], sigma2[3], eq_num[3]); t1=equation9_10(ff0[4], M_P2, t, sigma1[4], sigma2[4], eq_num[4]); t2=equation9_10(ff0[5], M_V2, t, sigma1[5], sigma2[5], eq_num[5]); t3=equation9_10(ff0[6], M_V2, t, sigma1[6], sigma2[6], eq_num[6]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n barBs -> K* transition form factors" // << std::endl; } // \bar B -> \bar K_1(1270) transition form factors // See the paper: H.Hatanaka and Kwei-Chou Yang, PRD78, 074007 (2008) if((parent == EvtPDL::getId(std::string("B+"))&& daught == EvtPDL::getId(std::string("K_1+")))|| (parent == EvtPDL::getId(std::string("B-"))&& daught == EvtPDL::getId(std::string("K_1-")))|| (parent == EvtPDL::getId(std::string("B0"))&& daught == EvtPDL::getId(std::string("K_10")))|| (parent == EvtPDL::getId(std::string("anti-B0"))&& daught == EvtPDL::getId(std::string("anti-K_10")))){ double ff0A[] ={0.450, 0.340, 0.41, 0.22, 0.31, 0.310, 0.28}; double sigma1A[]={1.600, 0.635, 1.51, 2.40, 2.01, 0.629, 1.36}; double sigma2A[]={0.974, 0.211, 1.18, 1.78, 1.50, 0.387, 0.72}; double ff0B[] ={-0.37, -0.29, -0.17, -0.45, -0.25, -0.250, -0.11}; double sigma1B[]={ 1.72, 0.729, 0.919, 1.34, 1.59, 0.378, -1.61}; double sigma2B[]={0.912, 0.074, 0.855, 0.69, 0.79, -0.755, 10.2}; int eq_num[]={ 10, 10, 10, 10, 10, 10, 10}; double MM2 =5.279*5.279; // GeV^2 double MB =5.279; // GeV double MK1 =1.272; // GeV double MK1A=1.31; // GeV double MK1B=1.34; // GeV double sinK=sin(thetaK); // sin(-34^o) double cosK=cos(thetaK); // cos(-34^o) double a, v0, v1, v2; a = sinK*equation9_10(ff0A[0], MM2, t, sigma1A[0], sigma2A[0], eq_num[0])*(MB+MK1)/(MB+MK1A); a =a+cosK*equation9_10(ff0B[0], MM2, t, sigma1B[0], sigma2B[0], eq_num[0])*(MB+MK1)/(MB+MK1B); v0= sinK*equation9_10(ff0A[1], MM2, t, sigma1A[1], sigma2A[1], eq_num[1])*MK1A/MK1; v0=v0+cosK*equation9_10(ff0B[1], MM2, t, sigma1B[1], sigma2B[1], eq_num[1])*MK1B/MK1; v1= sinK*equation9_10(ff0A[2], MM2, t, sigma1A[2], sigma2A[2], eq_num[2])*(MB+MK1A)/(MB+MK1); v1=v1+cosK*equation9_10(ff0B[2], MM2, t, sigma1B[2], sigma2B[2], eq_num[2])*(MB+MK1B)/(MB+MK1); v2= sinK*equation9_10(ff0A[3], MM2, t, sigma1A[3], sigma2A[3], eq_num[3])*(MB+MK1)/(MB+MK1A); v2=v2+cosK*equation9_10(ff0B[3], MM2, t, sigma1B[3], sigma2B[3], eq_num[3])*(MB+MK1)/(MB+MK1B); v =a; a0=v0; a1=v1; a2=v2; t1= sinK*equation9_10(ff0A[4], MM2, t, sigma1A[4], sigma2A[4], eq_num[4]); t1=t1+cosK*equation9_10(ff0B[4], MM2, t, sigma1B[4], sigma2B[4], eq_num[4]); t2= sinK*equation9_10(ff0A[5], MM2, t, sigma1A[5], sigma2A[5], eq_num[5])*(MB*MB-MK1A*MK1A)/(MB*MB-MK1*MK1); t2=t2+cosK*equation9_10(ff0B[5], MM2, t, sigma1B[5], sigma2B[5], eq_num[5])*(MB*MB-MK1B*MK1B)/(MB*MB-MK1*MK1); t3= sinK*equation9_10(ff0A[6], MM2, t, sigma1A[6], sigma2A[6], eq_num[6]); t3=t3+cosK*equation9_10(ff0B[6], MM2, t, sigma1B[6], sigma2B[6], eq_num[6]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n barB -> bar K_1(1270) transition form factors" // << std::endl; } // \bar B -> \bar K_1(1400) transition form factors // See the paper: H.Hatanaka and Kwei-Chou Yang, PRD78, 074007 (2008) if((parent == EvtPDL::getId(std::string("B+"))&& daught == EvtPDL::getId(std::string("K'_1+")))|| (parent == EvtPDL::getId(std::string("B-"))&& daught == EvtPDL::getId(std::string("K'_1-")))|| (parent == EvtPDL::getId(std::string("B0"))&& daught == EvtPDL::getId(std::string("K'_10")))|| (parent == EvtPDL::getId(std::string("anti-B0"))&& daught == EvtPDL::getId(std::string("anti-K'_10")))){ double ff0A[] ={0.450, 0.340, 0.41, 0.22, 0.31, 0.310, 0.28}; double sigma1A[]={1.600, 0.635, 1.51, 2.40, 2.01, 0.629, 1.36}; double sigma2A[]={0.974, 0.211, 1.18, 1.78, 1.50, 0.387, 0.72}; double ff0B[] ={-0.37, -0.29, -0.17, -0.45, -0.25, -0.250, -0.11}; double sigma1B[]={ 1.72, 0.729, 0.919, 1.34, 1.59, 0.378, -1.61}; double sigma2B[]={0.912, 0.074, 0.855, 0.69, 0.79, -0.755, 10.2}; int eq_num[]={ 10, 10, 10, 10, 10, 10, 10}; double MM2 =5.279*5.279; // GeV^2 double MB =5.279; // GeV double MK1 =1.403; // GeV double MK1A=1.31; // GeV double MK1B=1.34; // GeV double sinK=sin(thetaK); // sin(-34^o) double cosK=cos(thetaK); // cos(-34^o) double a, v0, v1, v2; a = cosK*equation9_10(ff0A[0], MM2, t, sigma1A[0], sigma2A[0], eq_num[0])*(MB+MK1)/(MB+MK1A); a =a-sinK*equation9_10(ff0B[0], MM2, t, sigma1B[0], sigma2B[0], eq_num[0])*(MB+MK1)/(MB+MK1B); v0= cosK*equation9_10(ff0A[1], MM2, t, sigma1A[1], sigma2A[1], eq_num[1])*MK1A/MK1; v0=v0-sinK*equation9_10(ff0B[1], MM2, t, sigma1B[1], sigma2B[1], eq_num[1])*MK1B/MK1; v1= cosK*equation9_10(ff0A[2], MM2, t, sigma1A[2], sigma2A[2], eq_num[2])*(MB+MK1A)/(MB+MK1); v1=v1-sinK*equation9_10(ff0B[2], MM2, t, sigma1B[2], sigma2B[2], eq_num[2])*(MB+MK1B)/(MB+MK1); v2= cosK*equation9_10(ff0A[3], MM2, t, sigma1A[3], sigma2A[3], eq_num[3])*(MB+MK1)/(MB+MK1A); v2=v2-sinK*equation9_10(ff0B[3], MM2, t, sigma1B[3], sigma2B[3], eq_num[3])*(MB+MK1)/(MB+MK1B); v =a; a0=v0; a1=v1; a2=v2; t1= cosK*equation9_10(ff0A[4], MM2, t, sigma1A[4], sigma2A[4], eq_num[4]); t1=t1-sinK*equation9_10(ff0B[4], MM2, t, sigma1B[4], sigma2B[4], eq_num[4]); t2= cosK*equation9_10(ff0A[5], MM2, t, sigma1A[5], sigma2A[5], eq_num[5])*(MB*MB-MK1A*MK1A)/(MB*MB-MK1*MK1); t2=t2-sinK*equation9_10(ff0B[5], MM2, t, sigma1B[5], sigma2B[5], eq_num[5])*(MB*MB-MK1B*MK1B)/(MB*MB-MK1*MK1); t3= cosK*equation9_10(ff0A[6], MM2, t, sigma1A[6], sigma2A[6], eq_num[6]); t3=t3-sinK*equation9_10(ff0B[6], MM2, t, sigma1B[6], sigma2B[6], eq_num[6]); models_counter=models_counter+1; // EvtGenReport(EVTGEN_NOTICE,"EvtGen") <<"\n The function EvtbTosllMSFF::getVectorFF(...) passed." // << "\n barB -> bar K_1(1270) transition form factors" // << std::endl; } // EvtGenReport(EVTGEN_NOTICE,"EvtGen") << "\n models_counter = " << models_counter // << "\n Vector form-factors at q^2 = " << t // << " for B -> V transition:" // << "\n v = " << v // << "\n a0 = " << a0 // << "\n a1 = " << a1 // << "\n a2 = " << a2 // << "\n t1 = " << t1 // << "\n t2 = " << t2 // << "\n t3 = " << t3 << std::endl; if(models_counter!=1){ EvtGenReport(EVTGEN_ERROR,"EvtGen") << "\n In the function EvtbTosllMSFF::getVectorFF(...) \n" << "the parameter models_counter not equal 1! \n" << "models_counter = " << models_counter < return m_u; // i=2 => return m_d; // i=3 => return m_s; // i=4 => return m_c; // i=5 => return m_b; double EvtbTosllMSFF::getQuarkMass(int i){ double qm=0.0; switch(i) { case 1: qm=0.23; // m_u break; case 2: qm=0.23; // m_d = m_u break; case 3: qm=0.35; // m_s break; case 4: qm=1.45; // m_c break; case 5: qm=4.85; // m_b break; default: EvtGenReport(EVTGEN_ERROR,"EvtGen") << "In the function EvtbTosllMSFF::getQuarkMass \n" << "the parameter i not equal 1, 2, 3, 4 or 5! \n" << "i =" << i < 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 #include #include #include #include #include #include #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 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:"< extraModels; #ifdef EVTGEN_EXTERNAL bool convertPythiaCodes(false); bool useEvtGenRandom(true); EvtExternalGenList genList(convertPythiaCodes, "", "gamma", useEvtGenRandom); radCorrEngine = genList.getPhotosModel(); extraModels = genList.getListOfModels(); #endif EvtGen myGenerator("../DECAY.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++setVectorSpinDensity(); myGenerator.generateDecay(root_part); root_part->printTree(); root_part->deleteTree(); }while (count++setVectorSpinDensity(); myGenerator.generateDecay(root_part); root_part->deleteTree(); }while (count++setVectorSpinDensity(); myGenerator.generateDecay(root_part); EvtParticle *p = root_part; do{ 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++Write(); file->Close(); EvtGenReport(EVTGEN_INFO,"EvtGen") << "SUCCESS\n"; } -void runMassCheck(int nevent, EvtGen &myGenerator, +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++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 <deleteTree(); }while (count++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++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++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++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++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++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;ireadDecayFile(std::string("../DECAY.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++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;ik2) { if ( tline>1) { dList[tline-3]=*dListItem; dListItem=new std::vector; 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; } 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()<setVectorSpinDensity(); } else{ root_part->setDiagonalSpinDensity(); } myGenerator.generateDecay(root_part); EvtParticle* p = root_part; std::vector fs=findFinalState(p); int j; for(j=0; j<(tline-1); j++) { std::vector temp=dList[j]; if ( temp.size() == fs.size() ) { bool foundIt=true; unsigned int k, l; std::vector alreadyUsed(temp.size()); for(k=0; kFill(0.5); //EvtGenReport(EVTGEN_INFO,"EvtGen") << "found one " << j << std::endl; dListNum[j]++;} } } root_part->deleteTree(); count++; }while(count findFinalState( EvtParticle *tree){ EvtParticle *p=tree; std::vector 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()<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++semileptonic events="<semileptonic events="<semileptonic events="<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; 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."< histo1(parser.getNToken()); std::vector histo2(parser.getNToken()); std::vector 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 temp(parser.getNToken()/2,0); std::vector tempB(parser.getNToken()/2,0); std::vector tempB0B(parser.getNToken()/2,0); std::vector tempB0(parser.getNToken()/2,0); std::vector tempBP(parser.getNToken()/2,0); std::vector tempBM(parser.getNToken()/2,0); std::vector tempD(parser.getNToken()/2,0); std::vector tempD0B(parser.getNToken()/2,0); std::vector tempD0(parser.getNToken()/2,0); std::vector tempDP(parser.getNToken()/2,0); std::vector 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()< 0 ) { root_part=EvtParticleFactory::particleFactory(UPS4,p_init); } else{ root_part=EvtParticleFactory::particleFactory(VPHO, p_init); } root_part->setVectorSpinDensity(); myGenerator.generateDecay(root_part); //EvtStdHep stdhep; //stdhep.init(); //root_part->makeStdHep(stdhep); //EvtGenReport(EVTGEN_INFO,"EvtGen") <printTree(); root_part->deleteTree(); }while (count++ 0 ) { root_part=EvtParticleFactory::particleFactory(UPS4,p_init); } else{ root_part=EvtParticleFactory::particleFactory(VPHO, p_init); } root_part->setVectorSpinDensity(); myGenerator.generateDecay(root_part); root_part->deleteTree(); }while (count++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++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++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 EP=EvtPDL::getId(std::string("e+")); static EvtId EM=EvtPDL::getId(std::string("e-")); 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 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 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++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 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 q2low(5); std::vector q2high(5); std::vector 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;jq2low[j]&&qsqFill( (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:"< q2low(5); std::vector q2high(5); std::vector 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;jq2low[j]&&qsqFill( (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:"< q2low(7); std::vector q2high(7); std::vector counts(7); int n(0); 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:"<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++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++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++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++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++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++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++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++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 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); + }while (count++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++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"<getId()==B0B){ tB0Btag->Fill(dt); outmix << dt << " -1"<Fill(EvtDecayAngle(p4_cp,p4_jpsi,p4_ep)); root_part->deleteTree(); }while (count++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++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:"<deleteTree(); }while (count++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_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++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++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")); 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++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")); 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++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")); 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++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 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:"<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:"<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++ 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="<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().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().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().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++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 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++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); + }while (count++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++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 "<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(0,200)); //tuple->column("gmass", gmass); //tuple->column("genergy", genergy); //HTValOrderedVector XDaugId, XDaugNephewId; //HTValOrderedVector XsDaugMass, XsDaugNephewMass; int nTot(0); for(int i=0;igetDaug(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;kgetDaug(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;ggetDaug(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(0,200)); //tuple->column("XsDaugNephewId", XDaugNephewId,"nTot", 0, "nTot"); //tuple->column("XsDaugNephewMass", XsDaugNephewMass,"nTot", 0, "nTot"); //tuple->dumpData(); root_part->deleteTree(); }while (count++Write(); file->Close(); EvtGenReport(EVTGEN_INFO,"EvtGen")<<"End EvtGen. Ran on "<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 "<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;igetDaug(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++Write(); file->Close(); EvtGenReport(EVTGEN_INFO,"EvtGen")<<"End EvtGen. Ran on "<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) { + 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"<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"<Fill(dt); outmix << dt << " 2"<Fill(dt); outmix << dt << " 0"<Fill(dt); outmix << dt << " 0"<deleteTree(); }while (count++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()<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++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()<deleteTree(); - }while (count++<10000); + }while (count++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++setDiagonalSpinDensity(); myGenerator.generateDecay(root_part); root_part->deleteTree(); }while (count++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:"<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); + }while (count++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++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++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++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); + }while (count++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); + }while (count++setVectorSpinDensity(); myGenerator.generateDecay(root_part); p=root_part; EvtParticle *otherb=root_part->getDaug(1); EvtGenReport(EVTGEN_INFO,"EvtGen") <<"Event:"<printTree(); outmix << "B" <<" "<getP4Lab().get(0) <<" "<getP4Lab().get(1) <<" "<getP4Lab().get(2) <<" "<getP4Lab().get(3)<getId()==PIP|| p->getId()==EP|| p->getId()==KP|| p->getId()==MUP){ outmix << "chg +1" <<" "<getP4Lab().get(1) <<" "<getP4Lab().get(2) <<" "<getP4Lab().get(3)<getId()==PIM|| p->getId()==EM|| p->getId()==KM|| p->getId()==MUM){ outmix << "chg -1" <<" "<getP4Lab().get(1) <<" "<getP4Lab().get(2) <<" "<getP4Lab().get(3)<getId()==GAMM|| p->getId()==K0L){ outmix << "neu" <<" "<getP4Lab().get(1) <<" "<getP4Lab().get(2) <<" "<getP4Lab().get(3)<nextIter(); }while(p!=0); outmix << "event"<deleteTree(); }while (count++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++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:"<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); + }while (count++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/test/example1.cc b/test/example1.cc index a70bfcf..66a56f9 100644 --- a/test/example1.cc +++ b/test/example1.cc @@ -1,92 +1,92 @@ // // Sample test program for running EvtGen // #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/EvtReport.hh" #include "EvtGenBase/EvtHepMCEvent.hh" #include "EvtGenBase/EvtSimpleRandomEngine.hh" #include "EvtGenBase/EvtMTRandomEngine.hh" #include "EvtGenBase/EvtAbsRadCorr.hh" #include "EvtGenBase/EvtDecayBase.hh" #ifdef EVTGEN_EXTERNAL #include "EvtGenExternal/EvtExternalGenList.hh" #endif #include #include #include -int main(int argc, char** argv) { +int main(int /*argc*/, char** /*argv*/) { EvtParticle* parent(0); // Define the random number generator EvtRandomEngine* eng = 0; #ifdef EVTGEN_CPP11 // Use the Mersenne-Twister generator (C++11 only) eng = new EvtMTRandomEngine(); #else eng = new EvtSimpleRandomEngine(); #endif EvtRandom::setRandomEngine(eng); EvtAbsRadCorr* radCorrEngine = 0; std::list extraModels; #ifdef EVTGEN_EXTERNAL bool convertPythiaCodes(false); bool useEvtGenRandom(true); EvtExternalGenList genList(convertPythiaCodes, "", "gamma", useEvtGenRandom); radCorrEngine = genList.getPhotosModel(); extraModels = genList.getListOfModels(); #endif //Initialize the generator - read in the decay table and particle properties EvtGen myGenerator("../DECAY.DEC","../evt.pdl", eng, radCorrEngine, &extraModels); //If I wanted a user decay file, I would read it in now. //myGenerator.readUDecay("../user.dec"); static EvtId UPS4 = EvtPDL::getId(std::string("Upsilon(4S)")); int nEvents(100); // Loop to create nEvents, starting from an Upsilon(4S) int i; for (i = 0; i < nEvents; i++) { std::cout<<"Event number "<setVectorSpinDensity(); // Generate the event myGenerator.generateDecay(parent); // Write out the results EvtHepMCEvent theEvent; theEvent.constructEvent(parent); HepMC::GenEvent* genEvent = theEvent.getEvent(); genEvent->print(std::cout); parent->deleteTree(); } delete eng; return 0; } diff --git a/test/exampleWriteHepMC.cc b/test/exampleWriteHepMC.cc index 7689f93..d086d65 100644 --- a/test/exampleWriteHepMC.cc +++ b/test/exampleWriteHepMC.cc @@ -1,123 +1,123 @@ // // Sample test program for running EvtGen // #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/EvtReport.hh" #include "EvtGenBase/EvtHepMCEvent.hh" #include "EvtGenBase/EvtSimpleRandomEngine.hh" #include "EvtGenBase/EvtMTRandomEngine.hh" #include "EvtGenBase/EvtAbsRadCorr.hh" #include "EvtGenBase/EvtDecayBase.hh" #ifdef EVTGEN_EXTERNAL #include "EvtGenExternal/EvtExternalGenList.hh" #endif #include #include #include #include bool filter(HepMC::GenEvent* event) { bool hasLepton = false; bool hasCharm = false; for (HepMC::GenEvent::particle_iterator it = event->particles_begin(); it != event->particles_end(); ++it) { if ( abs((*it)->pdg_id()) == 11 || abs((*it)->pdg_id()) == 13 ) { hasLepton = true; } int id = abs((*it)->pdg_id()); if ( id > 400 && id < 500 ) { hasCharm = true; } } return (hasLepton && (!hasCharm)); } -int main(int argc, char** argv) { +int main(int /*argc*/, char** /*argv*/) { EvtParticle* parent(0); // Define the random number generator EvtRandomEngine* eng = 0; #ifdef EVTGEN_CPP11 // Use the Mersenne-Twister generator (C++11 only) eng = new EvtMTRandomEngine(); #else eng = new EvtSimpleRandomEngine(); #endif EvtRandom::setRandomEngine(eng); EvtAbsRadCorr* radCorrEngine = 0; std::list extraModels; #ifdef EVTGEN_EXTERNAL bool convertPythiaCodes(false); bool useEvtGenRandom(true); EvtExternalGenList genList(convertPythiaCodes, "", "gamma", useEvtGenRandom); radCorrEngine = genList.getPhotosModel(); extraModels = genList.getListOfModels(); #endif //Initialize the generator - read in the decay table and particle properties EvtGen myGenerator("../DECAY.DEC","../evt.pdl", eng, radCorrEngine, &extraModels); //If I wanted a user decay file, I would read it in now. myGenerator.readUDecay("exampleFiles/Btousemileptonic.dec"); static EvtId UPS4 = EvtPDL::getId(std::string("Upsilon(4S)")); int nEvents(40000); std::ofstream hepmcFile("hepMCtest"); // Loop to create nEvents, starting from an Upsilon(4S) int i; for (i = 0; i < nEvents; i++) { std::cout<<"Event number "<setVectorSpinDensity(); // Generate the event myGenerator.generateDecay(parent); // Write out the results EvtHepMCEvent theEvent; theEvent.constructEvent(parent); HepMC::GenEvent* genEvent = theEvent.getEvent(); if ( filter(genEvent) ) { // genEvent->print(hepmcFile); hepmcFile<<(*genEvent); } parent->deleteTree(); } hepmcFile.close(); delete eng; return 0; }