Index: trunk/src/TwoBodyVectorMesonDecay.h
===================================================================
--- trunk/src/TwoBodyVectorMesonDecay.h (revision 0)
+++ trunk/src/TwoBodyVectorMesonDecay.h (revision 434)
@@ -0,0 +1,48 @@
+//==============================================================================
+// TwoBodyVectorMesonDecay.h
+//
+// Copyright (C) 2019 Tobias Toll and Thomas Ullrich
+//
+// This file is part of Sartre.
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation.
+// This program is distributed in the hope that it will be useful,
+// but without any warranty; without even the implied warranty of
+// merchantability or fitness for a particular purpose. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see .
+//
+// Author: Thomas Ullrich
+// $Date:$
+// $Author:$
+//==============================================================================
+#ifndef TwoBodyVectorMesonDecay_h
+#define TwoBodyVectorMesonDecay_h
+
+#include
+#include
+#include "TLorentzVector.h"
+#include "EventGeneratorSettings.h"
+#include "Event.h"
+
+class TwoBodyVectorMesonDecay {
+public:
+ TwoBodyVectorMesonDecay();
+
+ // Decay with polarization of the virtual photon taken into account (SCHC approximation).
+ pair decayVectorMeson(TLorentzVector& vm, Event& event, int daughterID);
+
+ // Simple 2-body decay flat in phase space
+ pair decayVectorMeson(TLorentzVector& vm, int daughterID);
+
+private:
+ double cosTheta(double, int);
+
+private:
+ TRandom3 *mRandom;
+ EventGeneratorSettings* mSettings;
+};
+#endif
Index: trunk/src/TwoBodyVectorMesonDecay.cpp
===================================================================
--- trunk/src/TwoBodyVectorMesonDecay.cpp (revision 0)
+++ trunk/src/TwoBodyVectorMesonDecay.cpp (revision 434)
@@ -0,0 +1,187 @@
+//==============================================================================
+// TwoBodyVectorMesonDecay.cpp
+//
+// Copyright (C) 2019 Tobias Toll and Thomas Ullrich
+//
+// This file is part of Sartre.
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation.
+// This program is distributed in the hope that it will be useful,
+// but without any warranty; without even the implied warranty of
+// merchantability or fitness for a particular purpose. See the
+// GNU General Public License for more details.
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see .
+//
+// Author: Thomas Ullrich
+// $Date:$
+// $Author:$
+//==============================================================================
+#include "TwoBodyVectorMesonDecay.h"
+#include "TGenPhaseSpace.h"
+#include "TVector3.h"
+#include
+#include
+
+
+TwoBodyVectorMesonDecay::TwoBodyVectorMesonDecay()
+{
+ mSettings = EventGeneratorSettings::instance();
+ mRandom = EventGeneratorSettings::randomGenerator();
+}
+
+pair
+TwoBodyVectorMesonDecay::decayVectorMeson(TLorentzVector& vm, int daughterID)
+{
+ //
+ // Simple 2-body decay flat in phase space
+ //
+ TGenPhaseSpace decay;
+
+ pair daughters;
+ double daughterMass = mSettings->lookupPDG(daughterID)->Mass();
+ double daughterMasses[2] = {daughterMass, daughterMass};
+
+ if (decay.SetDecay(vm, 2, daughterMasses)) {
+ double weight = decay.Generate(); // weight is always 1 here
+ if ((weight-1) > numeric_limits::epsilon()) {
+ cout << "TwoBodyVectorMesonDecay::decayVectorMeson(): Warning weight != 1, weight = " << weight << endl;
+ }
+ daughters.first = *decay.GetDecay(0);
+ daughters.second = *decay.GetDecay(1);
+ }
+ else {
+ cout << "TwoBodyVectorMesonDecay::decayVectorMeson(): Warning, kinematics of vector meson does not allow decay!" << endl;
+ }
+
+ return daughters;
+}
+
+pair
+TwoBodyVectorMesonDecay::decayVectorMeson(TLorentzVector& vm, Event& event, int daughterID)
+{
+ //
+ // Decay correctly treated with polarization of the virtual photon taken into account
+ // (SCHC approximation).
+ // This code was developed with the help of Athira Vijayakumar and Barak Schmookler
+ // from Stony Brook University.
+ //
+ // Essentially we create the daughters in the rest frame of the vector meson.
+ // Angels are generated randomly, flat in phi, and cos(theta) according to the
+ // respective distribution found in literature (see cosTheta()). At the end we
+ // boost back in the lab system.
+ //
+
+ pair daughters;
+
+ double daughterMass = mSettings->lookupPDG(daughterID)->Mass();
+
+ //
+ // Get matrix element from the cross-section ratio sig_L/sig_T.
+ // An alternative is to not use the ratio from Sartre directly but calculate
+ // it using Eq. 16 in PLB 449, 328. However, this is (i) model dependent and
+ // (ii) would introduce value that is not always identical to the one from
+ // Sartre. Tests by Athira showed that the two are very close so we use the
+ // actual value generated by Sartre, which is stored in the Event structure.
+ // Note also that the PLB version is t-integrated while Sartre's one is not.
+ //
+ double polarizationParameter = (1 - event.y)/(1 - event.y + (event.y*event.y)/2);
+ double a = event.crossSectionRatioLT * polarizationParameter;
+ double r0400 = a/(1+a);
+
+ //
+ // Generate random angles with the proper distributions
+ //
+ TRandom3 *random = EventGeneratorSettings::randomGenerator();
+ double phi = random->Uniform(2*M_PI);
+ double costh = cosTheta(r0400, daughterID);
+ if (fabs(costh) > 1) {
+ cout << "TwoBodyVectorMesonDecay::decayVectorMeson(): Error, falling back to simple decay scheme." << endl;
+ return decayVectorMeson(vm, daughterID);
+ }
+ double theta = acos(costh);
+
+ //
+ // Daughters in center-of-mass of vector meson,
+ // relying on both daughters having the same mass (and spin).
+ //
+ double p = sqrt(vm.M2()/4 - daughterMass*daughterMass);
+ double E = sqrt(p*p + daughterMass*daughterMass);
+ TVector3 p3(sin(theta)*cos(phi)*p, sin(theta)*sin(phi)*p, cos(theta)*p);
+ daughters.first = TLorentzVector(p3, E);
+ daughters.second = TLorentzVector(-p3, E);
+
+ //
+ // Boost to lab
+ //
+ TVector3 labVec = vm.Vect()*(1/vm.Energy());
+ daughters.first.Boost(labVec);
+ daughters.second.Boost(labVec);
+
+ return daughters;
+}
+
+double TwoBodyVectorMesonDecay::cosTheta(double r, int daughterID)
+{
+ //
+ // cos(theta) is generated according to the distributions derived by H1 et al.
+ // See: Eur. Phys. J. C 6, 603 (1999) and Eur. Phys. J. C 13, 371 (2000)
+ //
+ // To generate a random number we invert the cumulative of this distribution.
+ // It ends up with a (reduced) cubic equation that is solved analytically.
+ // We cover several cases since the solution scheme depends on the matrix
+ // element r (R^04_00).
+ //
+ // If something goes wrong we return 42 instead of -1 < cos(theta) < 1.
+ //
+ double p, q;
+
+ double rnd = mRandom->Uniform(1);
+
+ if (abs(daughterID) == 211 || abs(daughterID) == 321) { // spin 0 decay particle
+ p = (3-3*r)/(3*r-1);
+ q = 1+(3-3*r-4*rnd)/(3*r-1);
+ }
+ else if (abs(daughterID) == 11 || abs(daughterID) == 13) { // spin 1/2 decay particle
+ p = (3+3*r)/(1-3*r);
+ q = (4-8*rnd)/(1-3*r);
+ }
+ else {
+ cout << "TwoBodyVectorMesonDecay::cosTheta(): Cannot handle particle ID=" << daughterID << "." << endl;
+ return 42;
+ }
+
+ double R = sqrt(fabs(p)/3.)*(q < 0 ? -1 : 1);
+ double D = pow(p/3, 3) + q*q/4;
+
+ double result, phi;
+
+ if (p > 0) {
+ phi = asinh(q/(2*R*R*R));
+ result = -2*R*sinh(phi/3);
+ }
+ else if (p < 0) {
+ if (D <= 0) {
+ phi = acos(q/(2*R*R*R));
+ //result = -2*R*cos(phi/3);
+ //result = -2*R*cos(phi/3 + 2*M_PI/3);
+ result = -2*R*cos(phi/3 + 4*M_PI/3);
+ }
+ else {
+ phi = acosh(q/(2*R*R*R));
+ result = -2*R*cosh(phi/3);
+ }
+ }
+ else {
+ if (-q >= 0)
+ result = pow(-q, 1./3.);
+ else {
+ cout << "TwoBodyVectorMesonDecay::cosTheta(): Warning, cannot generate cos(theta)." << endl;
+ return 42;
+ }
+ }
+
+ return result;
+}