LauKinematics.cc

Go to the documentation of this file.

// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

// Authors:
// Thomas Latham
// John Back
// Paul Harrison

#include <iostream>
using std::cout;
using std::cerr;
using std::endl;

#include "TF2.h"
#include "TMath.h"
#include "TRandom.h"

#include "LauConstants.hh"
#include "LauKinematics.hh"
#include "LauRandom.hh"

ClassImp(LauKinematics)


LauKinematics::LauKinematics(Double_t m1, Double_t m2, Double_t m3, Double_t mParent, Bool_t calcSquareDPCoords) :
        m1_(m1), m2_(m2), m3_(m3), mParent_(mParent),
        m1Sq_(m1*m1), m2Sq_(m2*m2), m3Sq_(m3*m3), mParentSq_(mParent*mParent),
        mDTot_(m1 + m2 + m3),
        massInt_(mParent - (m1+m2+m3)),
        mSqDTot_(m1*m1 + m2*m2 + m3*m3),
        m12_(0.0), m23_(0.0), m13_(0.0),
        m12Sq_(0.0), m23Sq_(0.0), m13Sq_(0.0),
        c12_(0.0), c23_(0.0), c13_(0.0),
        mPrime_(0.0), thetaPrime_(0.0),
        qi_(0.0), qk_(0.0),
        p1_12_(0.0), p3_12_(0.0),
        p2_23_(0.0), p1_23_(0.0),
        p1_13_(0.0), p2_13_(0.0),
        p1_Parent_(0.0), p2_Parent_(0.0), p3_Parent_(0.0),
        squareDP_(calcSquareDPCoords), warnings_(kTRUE)
{
        // Constructor
        mass_.clear(); mMin_.clear(); mMax_.clear(); mSqMin_.clear(); mSqMax_.clear();
        mSq_.clear();
        mSqDiff_.clear();
        mDiff_.clear();

        mass_.push_back(m1_);
        mass_.push_back(m2_);
        mass_.push_back(m3_);

        mSq_.push_back(m1Sq_);
        mSq_.push_back(m2Sq_);
        mSq_.push_back(m3Sq_);

        // DP max and min kinematic boundary for circumscribed box
        // (see figure in PDG book).
        for (Int_t i = 0; i < 3; i++) {
                mMin_.push_back(mDTot_ - mass_[i]);
                mMax_.push_back(mParent_ - mass_[i]);
                mSqMin_.push_back(mMin_[i]*mMin_[i]);
                mSqMax_.push_back(mMax_[i]*mMax_[i]);
                mSqDiff_.push_back(mSqMax_[i] - mSqMin_[i]);
                mDiff_.push_back(mMax_[i] - mMin_[i]);
        }

        if (this->squareDP()) {
                cout<<"LauKinematics: squareDP = kTRUE."<<endl;
        } else {
                cout<<"LauKinematics: squareDP = kFALSE."<<endl;
        }
}

LauKinematics::~LauKinematics()
{
        // Destructor
}

void LauKinematics::updateKinematics(Double_t m13Sq, Double_t m23Sq)
{
        // Sets the internal private data members 
        // m13Sq and m23Sq, as well as m13 and m23.
        // Also calculate m12Sq and m12, given mParent defined in the constructor.
        // Lastly, calculate the helicity cosines.

        // Update the 3 mass-squares
        this->updateMassSquares(m13Sq, m23Sq);

        // Now update the helicity cosines
        this->calcHelicities();

        // Also calculate the m', theta' variables
        if (squareDP_) {this->calcSqDPVars();}
}

void LauKinematics::updateSqDPKinematics(Double_t mPrime, Double_t thetaPrime)
{
        // From square DP co-ordinates, calculate remaining kinematic variables
        this->updateSqDPMassSquares(mPrime, thetaPrime);

        // Finally calculate the helicities and track momenta
        this->calcHelicities();
}

void LauKinematics::calcSqDPVars()
{
        // For given m_12 and cos(theta_12) values, calculate m' and theta' for the square Dalitz plot
        Double_t value = (2.0*(m12_ - mMin_[2])/mDiff_[2]) - 1.0;
        mPrime_ = LauConstants::invPi*TMath::ACos(value);
        thetaPrime_ = LauConstants::invPi*TMath::ACos(c12_);
        // Sometimes events are assigned exactly thetaPrime = 0.0 or 1.0
        // which gives problems with efficiency and other histograms
        if (thetaPrime_ == 0.0)
        {
                thetaPrime_ += 1.0e-10;
        } else if (thetaPrime_ == 1.0)
        {
                thetaPrime_ -= 1.0e-10;
        }
}

Double_t LauKinematics::calcSqDPJacobian()
{
        // Calculate the Jacobian for the transformation 
        // m23^2, m13^2 -> m', theta' (square DP)

        Double_t e1Cms12 = (m12Sq_ - m2Sq_ + m1Sq_)/(2.0*m12_);
        Double_t e3Cms12 = (mParentSq_ - m12Sq_ - m3Sq_)/(2.0*m12_);

        Double_t p1Cms12 = this->pCalc(e1Cms12, m1Sq_);
        Double_t p3Cms12 = this->pCalc(e3Cms12, m3Sq_);

        Double_t deriv1 = LauConstants::piBy2*mDiff_[2]*TMath::Sin(LauConstants::pi*mPrime_);
        Double_t deriv2 = LauConstants::pi*TMath::Sin(LauConstants::pi*thetaPrime_);

        Double_t jacobian = 4.0*p1Cms12*p3Cms12*m12_*deriv1*deriv2;

        return jacobian;
}

void LauKinematics::updateMassSquares(Double_t m13Sq, Double_t m23Sq)
{
        m13Sq_ = m13Sq;
        if (m13Sq_ > 0.0) {
                m13_ = TMath::Sqrt(m13Sq_);
        } else {
                m13_ = 0.0;
        }

        m23Sq_ = m23Sq;
        if (m23Sq_ > 0.0) {
                m23_ = TMath::Sqrt(m23Sq_);
        } else {
                m23_ = 0.0;
        }

        // Now calculate m12Sq and m12.
        this->calcm12Sq();

        // Calculate momenta of tracks in parent (B, D etc.) rest-frame
        this->calcParentFrameMomenta();
}

void LauKinematics::updateSqDPMassSquares(Double_t mPrime, Double_t thetaPrime)
{
        // From square DP co-ordinates, calculate only the mass-squares

        // First set the square-DP variables
        mPrime_ = mPrime; thetaPrime_ = thetaPrime;

        // Next calculate m12 and c12
        Double_t m12 = 0.5*mDiff_[2]*(1.0 + TMath::Cos(LauConstants::pi*mPrime)) + mMin_[2];
        Double_t c12 = TMath::Cos(LauConstants::pi*thetaPrime);

        // From these calculate m13Sq and m23Sq
        this->updateMassSq_m12(m12, c12);

        // Calculate momenta of tracks in parent (B, D etc.) rest-frame
        this->calcParentFrameMomenta();
}

void LauKinematics::calcm12Sq()
{
        // Calculate m12Sq from m13Sq and m23Sq.
        m12Sq_ = mParentSq_ + mSqDTot_ - m13Sq_ - m23Sq_;

        // If m12Sq is too low, return lower limit,
        // and modify m13Sq accordingly.
        if (m12Sq_ < mSqMin_[2]) {
                m12Sq_ = mSqMin_[2] + 1.0e-3;
                m13Sq_ = mParentSq_ + mSqDTot_ - m12Sq_ - m23Sq_;
        }

        if (m12Sq_ > 0.0) {
                m12_ = TMath::Sqrt(m12Sq_);
        } else {
                m12_ = 0.0;
        }
}

void LauKinematics::calcParentFrameMomenta()
{
        Double_t e1 = (mParentSq_ + m1Sq_ - m23Sq_) / (2.0*mParent_);
        Double_t e2 = (mParentSq_ + m2Sq_ - m13Sq_) / (2.0*mParent_);
        Double_t e3 = (mParentSq_ + m3Sq_ - m12Sq_) / (2.0*mParent_);

        p1_Parent_ = this->pCalc(e1, m1Sq_);
        p2_Parent_ = this->pCalc(e2, m2Sq_);
        p3_Parent_ = this->pCalc(e3, m3Sq_);
}

void LauKinematics::calcHelicities()
{
        // Calculate helicity angle cosines, given m12Sq, m13Sq and m23Sq.
        // The notation is confusing for the helicity angle cosines:
        // cij is helicity angle for the pair which is made from tracks i and j.
        // It is the angle between tracks i and k in the ij rest frame  
        // Make sure that setMassSqVars is run first.
        Int_t zero(0), one(1), two(2);

        c12_ = cFromM(m12Sq_, m13Sq_, m12_, zero, one, two);
        p1_12_ = qi_; p3_12_ = qk_; // i, j = 12 (rest frame), k = 3
        c23_ = cFromM(m23Sq_, m12Sq_, m23_, one, two, zero);
        p2_23_ = qi_; p1_23_ = qk_; // i, j = 23 (rest frame), k = 1
        c13_ = cFromM(m13Sq_, m23Sq_, m13_, two, zero, one);
        p1_13_ = qi_; p2_13_ = qk_; // i, j = 31 (rest frame), k = 2

}

Double_t LauKinematics::cFromM(Double_t mijSq, Double_t mikSq, Double_t mij, Int_t i, Int_t j, Int_t k) const
{
        // Routine to calculate the cos(helicity) variables from the masses of the particles.
        // cij is helicity angle for the pair which is made from tracks i and j.
        // It is the angle between tracks i and k in the ij rest frame.
        Double_t EiCmsij = (mijSq - mSq_[j] + mSq_[i])/(2.0*mij);
        Double_t EkCmsij = (mParentSq_ - mijSq - mSq_[k])/(2.0*mij);

        if (EiCmsij < mass_[i]) {
                if (warnings_) {
                        cerr<<"Warning: EiCmsij = "<<EiCmsij<<" too small < mass_["<<i<<"] = "<<mass_[i]<<" in cFromM, i, j, k = "<<i<<", "<<j<<", "<<k<<endl;
                        cerr<<"mijSq = "<<mijSq<<"; mij = "<<mij<<"; mSq_["<<j<<"] = "<<mSq_[j]
                                <<"; mSq_["<<i<<"] = "<<mSq_[i]<<endl;
                }
                return 0.0;
        }
        if (EkCmsij < mass_[k]) {
                if (warnings_) {
                        cerr<<"Warning: EkCmsij = "<<EkCmsij<<" too small < mass_["<<k<<"] = "<<mass_[k]<<" in cFromM, i, j, k = "<<i<<", "<<j<<", "<<k<<endl;
                }
                return 0.0;
        }

        // Find track i and k momenta in ij rest frame
        qi_ = this->pCalc(EiCmsij, mSq_[i]);
        qk_ = this->pCalc(EkCmsij, mSq_[k]);

        // Find ij helicity angle
        Double_t cosHel = -(mikSq - mSq_[i] - mSq_[k] - 2.0*EiCmsij*EkCmsij)/(2.0*qi_*qk_);

        if (cosHel > 1.0) {
                cosHel = 1.0;
        } else if (cosHel < -1.0) {
                cosHel = -1.0;
        }

        if (i == 1) {cosHel *= -1.0;}

        return cosHel;

}

Double_t LauKinematics::mFromC(Double_t mijSq, Double_t cij, Double_t mij, Int_t i, Int_t j, Int_t k) const
{
        // Returns the mass-squared for a pair of particles, i,j, given their
        // invariant mass (squared) and the helicity angle.
        // cij is helicity angle for the pair which is made from tracks i and j.
        // It is the angle between tracks i and k in the ij rest frame.

        Double_t EiCmsij = (mijSq - mSq_[j] + mSq_[i])/(2.0*mij);
        Double_t EkCmsij = (mParentSq_ - mijSq - mSq_[k])/(2.0*mij);

        if (TMath::Abs(EiCmsij - mass_[i]) > 1e-6 && EiCmsij < mass_[i]) {
                if (warnings_) {
                        cerr<<"Warning: EiCmsij = "<<EiCmsij<<" < "<<mass_[i]<<" in mFromC, i, j, k = "<<i<<", "<<j<<", "<<k<<endl;
                }
                return 0.0;
        }
        if (TMath::Abs(EkCmsij - mass_[k]) > 1e-6 && EkCmsij < mass_[k]) {
                if (warnings_) {
                        cerr<<"Warning: EkCmsij = "<<EkCmsij<<" < "<<mass_[k]<<" in mFromC, i, j, k = "<<i<<", "<<j<<", "<<k<<endl;
                }
                return 0.0;
        }

        // Find track i and k momenta in ij rest fram
        Double_t qi = this->pCalc(EiCmsij, mSq_[i]);
        Double_t qk = this->pCalc(EkCmsij, mSq_[k]);

        // Find mikSq
        Double_t massSq = mSq_[i] + mSq_[k] + 2.0*EiCmsij*EkCmsij - 2.0*qi*qk*cij;

        if (massSq < mSqMin_[j]) {
                if (warnings_) {
                        cerr<<"Warning: mFromC below bound: i, j, k = "<<i<<", "<<j<<", "<<k<<endl;
                }
                massSq = mSqMin_[j];
        }

        return massSq;

}

void LauKinematics::genFlatPhaseSpace(Double_t& m13Sq, Double_t& m23Sq) const
{
        // Routine to generate flat phase-space events.
        // DP max kinematic boundaries in circumscribed box
        // See DP figure in PDG book.
        // m13max=mbrec-mass(2)
        // m13sqmax=m13max*m13max
        // m23max=mbrec-mass(1)
        // m23sqmax=m23max*m23max

        // Generate m13Sq and m23Sq flat within DP overall rectangular box
        // Loop if the generated point is not within the DP kinematic boundary

        do {
                m13Sq = mSqMin_[1] + LauRandom::randomFun()->Rndm()*mSqDiff_[1];
                m23Sq = mSqMin_[0] + LauRandom::randomFun()->Rndm()*mSqDiff_[0];

        } while ( ! this->withinDPLimits( m13Sq, m23Sq ) );
}

void LauKinematics::genFlatSqDP(Double_t& mPrime, Double_t& thetaPrime) const
{
        // Generate random event in the square Dalitz plot
        mPrime = LauRandom::randomFun()->Rndm();
        thetaPrime = LauRandom::randomFun()->Rndm();
}

Bool_t LauKinematics::withinDPLimits() const
{
        return this->withinDPLimits(m13Sq_,m23Sq_);
}

Bool_t LauKinematics::withinDPLimits(Double_t m13Sq, Double_t m23Sq) const
{
        // Find out whether the point (m13Sq,m23Sq) is within the limits of the
        // Dalitz plot. The limits are specified by the invariant masses
        // of the parent (e.g. B) and its three daughters that were
        // defined in the constructor of this class. Here
        // m_13Sq = square of invariant mass of daughters 1 and 3 
        // m_23Sq = square of invariant mass of daughters 2 and 3.

        Bool_t withinDP = kFALSE;

        // First check that m13Sq is within its absolute min and max
        if (!((m13Sq > mSqMin_[1]) && (m13Sq < mSqMax_[1]))) {
                return kFALSE;
        }

        // Now for the given value of m13Sq calculate the local min and max of m23Sq
        Double_t m13 = TMath::Sqrt(m13Sq);

        Double_t e3Cms13 = (m13Sq - m1Sq_ + m3Sq_)/(2.0*m13);
        Double_t p3Cms13 = this->pCalc(e3Cms13, m3Sq_);

        Double_t e2Cms13 = (mParentSq_ - m13Sq - m2Sq_)/(2.0*m13);
        Double_t p2Cms13 = this->pCalc(e2Cms13, m2Sq_);

        Double_t term = 2.0*e2Cms13*e3Cms13 + m2Sq_ + m3Sq_;

        Double_t m23SqLocMin = term - 2.0*p2Cms13*p3Cms13;
        Double_t m23SqLocMax = term + 2.0*p2Cms13*p3Cms13;

        // Check whether the given value of m23Sq satisfies these bounds
        if (m23Sq > m23SqLocMin && m23Sq < m23SqLocMax) {    
                withinDP = kTRUE;
        }

        return withinDP;
}

Bool_t LauKinematics::withinDPLimits2(Double_t m13Sq, Double_t m23Sq) const
{
        // Same as withinDPLimits, but this time testing whether the m13Sq 
        // variable is within the kinematic boundary for the given m23Sq value

        Bool_t withinDP = kFALSE;

        // First check that m13Sq is within its absolute min and max
        if (!((m23Sq > mSqMin_[0]) && (m23Sq < mSqMax_[0]))) {
                return kFALSE;
        }

        // Now for the given value of m13Sq calculate the local min and max of m23Sq
        Double_t m23 = TMath::Sqrt(m23Sq);

        Double_t e3Cms23 = (m23Sq - m2Sq_ + m3Sq_)/(2.0*m23);
        Double_t p3Cms23 = this->pCalc(e3Cms23, m3Sq_);

        Double_t e1Cms23 = (mParentSq_ - m23Sq - m1Sq_)/(2.0*m23);
        Double_t p1Cms23 = this->pCalc(e1Cms23, m1Sq_);

        Double_t term = 2.0*e3Cms23*e1Cms23 + m1Sq_ + m3Sq_;

        Double_t m13SqLocMin = term - 2.0*p3Cms23*p1Cms23;
        Double_t m13SqLocMax = term + 2.0*p3Cms23*p1Cms23;

        // Check whether the given value of m23Sq satisfies these bounds
        if (m13Sq > m13SqLocMin && m13Sq < m13SqLocMax) {
                withinDP = kTRUE;
        }

        return withinDP;
}

Bool_t LauKinematics::withinSqDPLimits(Double_t mPrime, Double_t thetaPrime) const
{
        // Check whether m' and theta' are between 0 and 1
        Bool_t withinDP(kFALSE);
        if (mPrime > 0.0 && mPrime < 1.0 &&
                        thetaPrime > 0.0 && thetaPrime < 1.0) {
                withinDP = kTRUE;
        }

        return withinDP;
}

Double_t LauKinematics::calcThirdMassSq(Double_t firstMassSq, Double_t secondMassSq) const
{
        // Calculate one massSq from the other two
        return mParentSq_ + mSqDTot_ - firstMassSq - secondMassSq;
}

Double_t LauKinematics::distanceFromDPCentre() const
{
        return this->distanceFromDPCentre(m13Sq_,m23Sq_);
}

Double_t LauKinematics::distanceFromDPCentre(Double_t m13Sq, Double_t m23Sq) const
{
        // DP centre is defined as the point where m12=m13=m23

        // First find the m^2_ij value for the centre
        Double_t centreMijSq = (mParentSq_ + m1Sq_ + m2Sq_ + m3Sq_)/3.0;

        // Then find the difference between this and the two provided co-ords
        Double_t diff13 = m13Sq - centreMijSq;
        Double_t diff23 = m23Sq - centreMijSq;

        // Calculate the total distance
        Double_t distance = TMath::Sqrt(diff13*diff13 + diff23*diff23);
        return distance;
}

Double_t LauKinematics::pCalc(Double_t energy, Double_t massSq) const
{
        // Routine to calculate the momentum of a particle, given its energy and 
        // invariant mass (squared).
        Double_t arg = energy*energy - massSq;

        if (arg < 0.0) {
                //if (warnings_) {
                        //cerr<<"Warning. In pcalc, arg < 0.0: "<<arg<<" for e = "<<energy<<", mSq = "<<massSq<<endl;
                //}
                arg = 0.0;
        }
        Double_t pCalcVal = TMath::Sqrt(arg);
        return pCalcVal;

}

void LauKinematics::flipAndUpdateKinematics()
{
        // Flips the DP variables m13^2 <-> m23^2.
        // Used in the case of symmetrical Dalitz plots (i.e. when two of
        // the daughter tracks are the same type) within the Dynamics()
        // function.
        this->updateKinematics(m23Sq_, m13Sq_);

}

void LauKinematics::updateMassSq_m23(Double_t m23, Double_t c23)
{
        // Update the variables m12Sq_ and m13Sq_ given m23 and c23.
        m23_ = m23; m23Sq_ = m23*m23; c23_ = c23;

        Int_t zero(0), one(1), two(2);
        m13Sq_ = this->mFromC(m23Sq_, -c23_, m23_, two, one, zero);
        m12Sq_ = mParentSq_ - m23Sq_ - m13Sq_ + mSqDTot_;

        m13_ = TMath::Sqrt(m13Sq_);
        m12_ = TMath::Sqrt(m12Sq_);
        //c12_ = cFromM(m12Sq_, m13Sq_, m12_, zero, one, two);
        //c13_ = cFromM(m13Sq_, m23Sq_, m13_, two, zero, one);
}

void LauKinematics::updateMassSq_m13(Double_t m13, Double_t c13)
{
        // Update the variables m12Sq_ and m23Sq_ given m13 and c13.
        m13_ = m13; m13Sq_ = m13*m13; c13_ = c13;

        Int_t zero(0), one(1), two(2);
        m23Sq_ = this->mFromC(m13Sq_, c13_, m13_, two, zero, one);
        m12Sq_ = mParentSq_ - m23Sq_ - m13Sq_ + mSqDTot_;

        m23_ = TMath::Sqrt(m23Sq_);
        m12_ = TMath::Sqrt(m12Sq_);
        //c23_ = cFromM(m23Sq_, m12Sq_, m23_, one, two, zero);
        //c12_ = cFromM(m12Sq_, m13Sq_, m12_, zero, one, two);
}

void LauKinematics::updateMassSq_m12(Double_t m12, Double_t c12)
{
        // Update the variables m23Sq_ and m13Sq_ given m12 and c12.
        m12_ = m12; m12Sq_ = m12*m12; c12_ = c12;

        Int_t zero(0), one(1), two(2);
        m13Sq_ = this->mFromC(m12Sq_, c12_, m12_, zero, one, two);
        m23Sq_ = mParentSq_ - m12Sq_ - m13Sq_ + mSqDTot_;
        m13_ = TMath::Sqrt(m13Sq_);
        m23_ = TMath::Sqrt(m23Sq_);
        //c23_ = cFromM(m23Sq_, m12Sq_, m23_, one, two, zero);
        //c13_ = cFromM(m13Sq_, m23Sq_, m13_, two, zero, one);
}

Double_t LauKinematics::genm13Sq() const
{
        Double_t m13Sq = mSqMin_[1] + LauRandom::randomFun()->Rndm()*mSqDiff_[1];
        return m13Sq;
}

Double_t LauKinematics::genm23Sq() const
{
        Double_t m23Sq = mSqMin_[0] + LauRandom::randomFun()->Rndm()*mSqDiff_[0];
        return m23Sq;
}

Double_t LauKinematics::genm12Sq() const
{
        Double_t m12Sq = mSqMin_[2] + LauRandom::randomFun()->Rndm()*mSqDiff_[2];
        return m12Sq;
}


void LauKinematics::drawDPContour(Int_t orientation, Int_t nbins)
{
        // orientation -
        // 1323 : x = m13, y = m23
        // etc.

        Double_t m1 = this->getm1();
        Double_t m2 = this->getm2();
        Double_t m3 = this->getm3();
        Double_t mParent = this->getmParent();

        Double_t m13SqMin = this->getm13SqMin();
        Double_t m23SqMin = this->getm23SqMin();
        Double_t m12SqMin = this->getm12SqMin();
        Double_t m13SqMax = this->getm13SqMax();
        Double_t m23SqMax = this->getm23SqMax();
        Double_t m12SqMax = this->getm12SqMax();

        Double_t xMin(0.0);
        Double_t xMax(mParent*mParent);
        Double_t yMin(0.0);
        Double_t yMax(mParent*mParent);
        if (orientation == 1323) {
                xMin = m13SqMin-1.0; xMax = m13SqMax+1.0;
                yMin = m23SqMin-1.0; yMax = m23SqMax+1.0;
        } else if (orientation == 2313) {
                xMin = m23SqMin-1.0; xMax = m23SqMax+1.0;
                yMin = m13SqMin-1.0; yMax = m13SqMax+1.0;
        } else if (orientation == 1213) {
                xMin = m12SqMin-1.0; xMax = m12SqMax+1.0;
                yMin = m13SqMin-1.0; yMax = m13SqMax+1.0;
        } else if (orientation == 1312) {
                xMin = m13SqMin-1.0; xMax = m13SqMax+1.0;
                yMin = m12SqMin-1.0; yMax = m12SqMax+1.0;
        } else if (orientation == 1223) {
                xMin = m12SqMin-1.0; xMax = m12SqMax+1.0;
                yMin = m23SqMin-1.0; yMax = m23SqMax+1.0;
        } else if (orientation == 2312) {
                xMin = m23SqMin-1.0; xMax = m23SqMax+1.0;
                yMin = m12SqMin-1.0; yMax = m12SqMax+1.0;
        } else {
                cerr<<"Unrecognised orientation, not drawing contour."<<endl;
                return;
        }

        Int_t npar = 4;
        TF2 * f2 = new TF2("contour", &dal, xMin, xMax, yMin, yMax, npar);

        // Set the parameters
        f2->SetParameter(0,mParent);
        if (orientation == 1323) {
                f2->SetParameter(1,m1);
                f2->SetParameter(2,m2);
                f2->SetParameter(3,m3);
        } else if (orientation == 2313) {
                f2->SetParameter(1,m2);
                f2->SetParameter(2,m1);
                f2->SetParameter(3,m3);
        } else if (orientation == 1213) {
                f2->SetParameter(1,m2);
                f2->SetParameter(2,m3);
                f2->SetParameter(3,m1);
        } else if (orientation == 1312) {
                f2->SetParameter(1,m3);
                f2->SetParameter(2,m2);
                f2->SetParameter(3,m1);
        } else if (orientation == 1223) {
                f2->SetParameter(1,m1);
                f2->SetParameter(2,m3);
                f2->SetParameter(3,m2);
        } else if (orientation == 2312) {
                f2->SetParameter(1,m3);
                f2->SetParameter(2,m1);
                f2->SetParameter(3,m2);
        }

        // Set up the contour to be drawn when the value of the function == 1.0
        Double_t b[]={1.0};
        f2->SetContour(1,b);

        // Set the number of bins for the contour to be sampled over
        f2->SetNpx(nbins);
        f2->SetNpy(nbins);
        // and the line style
        f2->SetLineWidth(3);
        f2->SetLineStyle(kSolid);

        // Draw the contour on top of the histo that should already have been drawn
        f2->DrawCopy("same");

        delete f2;
}

Double_t dal(Double_t* x, Double_t* par)
{
        Double_t mParent = par[0];
        Double_t mi = par[1];
        Double_t mj = par[2];
        Double_t mk = par[3];

        Double_t mikSq=x[0];
        Double_t mjkSq=x[1];
        Double_t mik = TMath::Sqrt(mikSq);
        Double_t mjk = TMath::Sqrt(mjkSq);

        Double_t ejcmsik = (mParent*mParent-mj*mj-mik*mik)/(2.0*mik);
        Double_t ekcmsik = (mik*mik+mk*mk-mi*mi)/(2.0*mik);
        if (ekcmsik<mk || ejcmsik<mj) return 2.0;

        Double_t pj = TMath::Sqrt(ejcmsik*ejcmsik-mj*mj);
        Double_t pk = TMath::Sqrt(ekcmsik*ekcmsik-mk*mk);
        Double_t coshelik = (mjk*mjk - mk*mk - mj*mj - 2.0*ejcmsik*ekcmsik)/(2.0*pj*pk);

        Double_t coshelikSq = coshelik*coshelik;
        return coshelikSq;
}