text: // StringLength.cc is a part of the PYTHIA event generator.
  // Copyright (C) 2015 Torbjorn Sjostrand.
  // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
  // Please respect the MCnet Guidelines, see GUIDELINES for details.
  
  // Function definitions (not found in the header) for the
  // StringLength class.
  
  // Calculate the lambda measure for string and junctions.
  
  #include "Pythia8/StringLength.h"
  
  namespace Pythia8 {
  
  //==========================================================================
  
  // The StringLength class.
  
  //--------------------------------------------------------------------------
  
  // Constants: could be changed here if desired, but normally should not.
  // These are of technical nature, as described for each.
  
  // Minimum delta R between two partons. This is to avoid problems
  // with infinities.
  const double StringLength::MINDELTAR = 1e-20;
  const double StringLength::TINY = 1e-20;
  
  void StringLength::init(Info* infoPtrIn, Settings& settings) {
  
    // Save pointers.
    infoPtr = infoPtrIn;
  
    // Store variables.
    m0         = settings.parm("ColourReconnection:m0");
    m0sqr      = pow2(m0);
    juncCorr   = settings.parm("ColourReconnection:junctionCorrection");
    sqrt2      = sqrt(2);
    lambdaForm = settings.mode("ColourReconnection:lambdaForm");
  }
  
  //--------------------------------------------------------------------------
  
  // Calculate string length for two indices in the event record.
  
  double StringLength::getStringLength( Event& event, int i, int j) {
  
    // Find rest frame of particles.
    Vec4 p1 =  event[i].p();
    Vec4 p2 =  event[j].p();
  
    return getStringLength(p1, p2);
  }
  
  //--------------------------------------------------------------------------
  
  // Calculate string length for two particles given their four-momenta.
  
  double StringLength::getStringLength( Vec4 p1, Vec4 p2) {
  
    // Check that particles are not completely paralel.
    if (REtaPhi(p1,p2) < MINDELTAR) {
      return 1e9;
    }
    
    // Check for very small energies.
    if (p1.e() < TINY || p2.e() < TINY) return 1e9;
    
    // Boost to restframe.
    Vec4 pSum = p1 + p2;
    p1.bstback(pSum);
    p2.bstback(pSum);
  
    // Calculate string length.
    Vec4 p0(0,0,0,1.);
  
    return getLength(p1,p0) + getLength(p2,p0);
  }
  
  //--------------------------------------------------------------------------
  
  // Calculate string length of a single particle.
  // The first vector is the 4 vector of the particle.
  // The second vector represents (1,0,0,0) in dipole restframe.
  
  double StringLength::getLength(Vec4 p, Vec4 v, bool isJunc) {
    double m = m0;
    if (isJunc) m *= juncCorr;
  
    if (lambdaForm == 0)
      return log (1. + sqrt2 * v * p/ m );
    else if (lambdaForm == 1)
      return log (1. + 2 * v * p/ m );
    else if (lambdaForm == 2)
      return log (2 * v * p / m);
    else
      return 1e9;
  }
  
  //--------------------------------------------------------------------------
  
  // Calculate the length of a single junction given the 3 entries in the event.
  
  double StringLength::getJuncLength( Event& event, int i, int j, int k) {
    if (i == j || i == k || j == k)
      return 1e9;
  
    Vec4 p1 = event[i].p();
    Vec4 p2 = event[j].p();
    Vec4 p3 = event[k].p();
  
    return getJuncLength(p1, p2, p3);
  }
  //--------------------------------------------------------------------------
  
  // Calculate the length of a single junction given the 3 four-momenta.
  
  double StringLength::getJuncLength(Vec4 p1, Vec4 p2, Vec4 p3) {
  
    if( p1*p2 > 1e-5 || p2*p3 > 1e-5 || p3*p1 > 1e-5) return 1e9;
    // Check for parallel particles.
    if (REtaPhi(p1,p2) < MINDELTAR || REtaPhi(p1,p3) < MINDELTAR ||
        REtaPhi(p2,p3) < MINDELTAR) {
      return 1e9;
    }
    
    // Check for very small energies.
    if (p1.e() < TINY || p2.e() < TINY || p3.e() < TINY) return 1e9;
    
    // Find the junction rest frame.
    RotBstMatrix MfromJRF1 = stringFragmentation.junctionRestFrame(p1,p2,p3);
    MfromJRF1.invert();
    Vec4 v1(0,0,0,1);
    v1.rotbst(MfromJRF1);
  
    // Possible problem when the right system rest frame system is not found.
    if (pow2(p1*v1) - p1*p1 < 0 || pow2(p2*v1) - p2*p2 < 0
      || pow2(p3*v1) - p3*p3 < 0)
      return 1e9;
  
    // Calcualte the junction length.
    return getLength(p1, v1, true) + getLength(p2, v1, true)
      + getLength(p3, v1, true);
  }
  
  //--------------------------------------------------------------------------
  
  // Calculate the length of a double junction given the 4 entries in the event.
  // The first two are expected to be quarks, the second two to be anti quarks.
  
  double StringLength::getJuncLength( Event& event, int i, int j, int k, int l) {
    if (i == j || i == k || i == l || j == k || j == l || k == l)
      return 1e9;
  
    // Simple minimum check of lengths.
    double origLength = getStringLength(event, i, k) +
      getStringLength(event, j, l);
    double minLength  = getStringLength(event, i, j) +
      getStringLength(event, k, l);
  
    if (origLength < minLength)
      return minLength;
  
    Vec4 p1 = event[i].p();
    Vec4 p2 = event[j].p();
    Vec4 p3 = event[k].p();
    Vec4 p4 = event[l].p();
  
    return getJuncLength(p1, p2, p3, p4);
  }
  
  //--------------------------------------------------------------------------
  
  // Calculate the length of a double junction given the 4 four-momenta.
  // The first two are expected to be quarks, the second two to be anti quarks.
  
  double StringLength::getJuncLength(Vec4 p1, Vec4 p2, Vec4 p3, Vec4 p4) {
    // Check for parallel problems.
    if (REtaPhi(p1,p2) < MINDELTAR || REtaPhi(p1,p3) < MINDELTAR ||
        REtaPhi(p1,p4) < MINDELTAR || REtaPhi(p2,p3) < MINDELTAR ||
        REtaPhi(p2,p4) < MINDELTAR || REtaPhi(p3,p4) < MINDELTAR) {
      return 1e9;
    }
    
    // Check for very small energies.
    if (p1.e() < TINY || p2.e() < TINY || p3.e() < TINY || p4.e() < TINY) 
      return 1e9;
  
    // Calculate velocity of first junction.
    Vec4 pSum1 = p3 +p4;
    RotBstMatrix MfromJRF1 = stringFragmentation.junctionRestFrame(p1,p2,pSum1);
    MfromJRF1.invert();
    Vec4 v1(0,0,0,1);
    v1.rotbst(MfromJRF1);
  
    // Calculate velocity of second junction.
    Vec4 pSum2 = p1 + p2;
    RotBstMatrix MfromJRF2 = stringFragmentation.junctionRestFrame(p3,p4,pSum2);
    MfromJRF2.invert();
    Vec4 v2(0,0,0,1);
    v2.rotbst(MfromJRF2);
  
    // This only happens if it is not possible to find the correct rest frame.
    if (pow2(p1*v1) - p1*p1 < 0 || pow2(p2*v1) - p2*p2 < 0 ||
        pow2(p3*v2) - p3*p3 < 0 || pow2(p4*v2) - p4*p4 < 0)
      return 1e9;
  
    double returnValue = getLength(p1, v1, true) + getLength(p2, v1, true)
      + getLength(p3, v2, true) + getLength(p4, v2, true);
    
    if (pow2(v1*v2)-1 > 0)
      returnValue += log(v1*v2 + sqrt(pow2(v1*v2)-1));
    else 
      returnValue += log(v1*v2);
      
    return returnValue;
  }
  
  //==========================================================================
  
  } // end namespace Pythia8
  