nBodyPhaseSpaceGen.h

//
//    Copyright 2010
//
//    This file is part of starlight.
//
//    starlight 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, either version 3 of the License, or
//    (at your option) any later version.
//        
//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
//
//
// File and Version Information:
// $Rev::                             $: revision of last commit
// $Author: jwebb $: author of last commit
// $Date: 2012/11/27 22:27:32 $: date of last commit
//
// Description:
//    calculates n-body phase space (constant matrix element) using various algorithms
//    
//    the n-body decay is split up into (n - 2) successive 2-body decays
//    each 2-body decay is considered in its own center-of-mass frame thereby
//    separating the mass from the (trivial) angular dependence
//    
//    the event is boosted into the same frame in which the n-body system is
//    given
//    
//    based on:
//    GENBOD (CERNLIB W515), see F. James, "Monte Carlo Phase Space", CERN 68-15 (1968)
//    NUPHAZ, see M. M. Block, "Monte Carlo phase space evaluation", Comp. Phys. Commun. 69, 459 (1992)
//    S. U. Chung, "Spin Formalism", CERN Yellow Report
//    S. U. Chung et. al., "Diffractive Dissociation for COMPASS"
//    
//    index convention:
//    - all vectors have the same size (= number of decay daughters)
//    - index i corresponds to the respective value in the (i + 1)-body system: effective mass M, break-up momentum, angles
//    - thus some vector elements are not used like breakupMom[0], theta[0], phi[0], ...
//      this overhead is negligible compared to the ease of notation
//    
//    the following graph illustrates how the n-body decay is decomposed into a sequence of two-body decays
//    
//    n-body       ...                   3-body                 2-body                 single daughter
//    
//    m[n - 1]                           m[2]                   m[1]
//     ^                                  ^                      ^
//     |                                  |                      |
//     |                                  |                      |
//    M[n - 1] --> ... -->               M[2] -->               M[1] -->               M    [0] = m[0]
//    theta[n - 1] ...                   theta[2]               theta[1]               theta[0] = 0 (not used)
//    phi  [n - 1] ...                   phi  [2]               phi  [1]               phi  [0] = 0 (not used)
//    mSum [n - 1] ...                   mSum [2]               mSum [1]               mSum [0] = m[0]
//    = sum_0^(n - 1) m[i]               = m[2] + m[1] + m[0]   = m[1] + m[0]
//    breakUpMom[n - 1] ...              breakUpMom[2]          breakUpMom[1]          breakUpMom[0] = 0 (not used)
//    = q(M[n - 1], m[n - 1], M[n - 2])  = q(M[2], m[2], M[1])  = q(M[1], m[1], m[0])
//
//


#ifndef NBODYPHASESPACEGEN_H
#define NBODYPHASESPACEGEN_H


#include <iostream>
#include <vector>

#include "reportingUtils.h"
#include "lorentzvector.h"
#include "randomgenerator.h"
#include "starlightconstants.h"


// small helper functions
// calculates factorial
inline
unsigned int
factorial(const unsigned int n)
{
        unsigned int fac = 1;
        for (unsigned int i = 1; i <= n; ++i)
                fac *= i;
        return fac;
}


// computes breakup momentum of 2-body decay
inline
double
breakupMomentum(const double M,   // mass of mother particle
                const double m1,  // mass of daughter particle 1
                const double m2)  // mass of daughter particle 2
{
  if (M < m1 + m2)
    return 0;
  return sqrt((M - m1 - m2) * (M + m1 + m2) * (M - m1 + m2) * (M + m1 - m2)) / (2 * M);
}


class nBodyPhaseSpaceGen {

public:

        nBodyPhaseSpaceGen();
        virtual ~nBodyPhaseSpaceGen();
  
        // generator setup
        bool setDecay(const std::vector<double>& daughterMasses);  // daughter particle masses
        bool setDecay(const unsigned int         nmbOfDaughters,   // number of daughter particles
                      const double*              daughterMasses);  // array of daughter particle masses
  
        // random generator
        void   setSeed(const unsigned int seed) { _rnd.SetSeed(seed);  }  
        double random ()                        { return _rnd.Rndom(); }  

        // high-level generator interface
        double generateDecay        (const lorentzVector& nBody);          // Lorentz vector of n-body system in lab frame
        bool   generateDecayAccepted(const lorentzVector& nBody,           // Lorentz vector of n-body system in lab frame
                                     const double         maxWeight = 0);  // if positive, given value is used as maximum weight, otherwise _maxWeight

        void   setMaxWeight          (const double maxWeight) { _maxWeight = maxWeight;    }  
        double maxWeight             () const                 { return _maxWeight;         }  
        double normalization         () const                 { return _norm;              }  
        double eventWeight           () const                 { return _weight;            }  
        double maxWeightObserved     () const                 { return _maxWeightObserved; }  
        void   resetMaxWeightObserved()                       { _maxWeightObserved = 0;    }  

        double estimateMaxWeight(const double       nBodyMass,                 // sic!
                                 const unsigned int nmbOfIterations = 10000);  // number of generated events

        inline bool eventAccepted(const double maxWeight = 0);

        //----------------------------------------------------------------------------
        // trivial accessors
        const lorentzVector&              daughter        (const int index) const { return _daughters[index];  }  
        const std::vector<lorentzVector>& daughters       ()                const { return _daughters;         }  
        unsigned int                      nmbOfDaughters  ()                const { return _n;                 }  
        double                            daughterMass    (const int index) const { return _m[index];          }  
        double                            intermediateMass(const int index) const { return _M[index];          }  
        double                            breakupMom      (const int index) const { return _breakupMom[index]; }  
        double                            cosTheta        (const int index) const { return _cosTheta[index];   }  
        double                            phi             (const int index) const { return _phi[index];        }  


        std::ostream& print(std::ostream& out = std::cout) const;  
        friend std::ostream& operator << (std::ostream&             out,
                                          const nBodyPhaseSpaceGen& gen)
        { return gen.print(out); }

private:

        //----------------------------------------------------------------------------
        // low-level generator interface
        void pickMasses(const double nBodyMass);  // total energy of n-body system in its RF

        double calcWeight();

        inline void pickAngles();
                
        void calcEventKinematics(const lorentzVector& nBody);  // Lorentz vector of n-body system in lab frame

        // external parameters
        std::vector<double> _m;  

        // internal variables
        unsigned int               _n;                  
        std::vector<double>        _M;                  
        std::vector<double>        _cosTheta;           
        std::vector<double>        _phi;                
        std::vector<double>        _mSum;               
        std::vector<double>        _breakupMom;         
        std::vector<lorentzVector> _daughters;          
        double                     _norm;               
        double                     _weight;             
        double                     _maxWeightObserved;  
        double                     _maxWeight;          

        randomGenerator _rnd;  

};


inline
void
nBodyPhaseSpaceGen::pickAngles()
{
        for (unsigned int i = 1; i < _n; ++i) {  // loop over 2- to n-bodies
                _cosTheta[i] = 2 * random() - 1;  // range [-1,    1]
                _phi[i]      = starlightConstants::twoPi * random();  // range [ 0, 2 pi]
        }
}


inline 
bool
nBodyPhaseSpaceGen::eventAccepted(const double maxWeight)  // if maxWeight > 0, given value is used as maximum weight, otherwise _maxWeight
{
        const double max = (maxWeight <= 0) ? _maxWeight : maxWeight;
        if (max <= 0) {
                printWarn << "maximum weight = " << max << " does not make sense. rejecting event." << std::endl;
                return false;
        }
        if ((_weight / max) > random())
                return true;
        return false;
}


#endif  // NBODYPHASESPACEGEN_H