G4PomeronCrossSection.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// $Id: G4PomeronCrossSection.cc 79761 2014-03-13 15:14:37Z gcosmo $
//

#include "G4PomeronCrossSection.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "G4Pow.hh"                                      // factorial()

G4PomeronCrossSection::G4PomeronCrossSection() :
   pomeron_Alpha(0), pomeron_Alpha_Hard(0), pomeron_Alphaprime(0),
   pomeron_C(0), pomeron_Gamma(0), pomeron_Gamma_Hard(0),
   pomeron_Rsquare(0), pomeron_S(0)
{}


G4PomeronCrossSection::~G4PomeronCrossSection()
{;}
//**********************************************************************************************

G4PomeronCrossSection::G4PomeronCrossSection(const G4ParticleDefinition * particle)
    {
    G4int Encoding = std::abs(particle->GetPDGEncoding());

    if (std::abs(particle->GetBaryonNumber())!=0)
        InitForNucleon();
    else if (Encoding/100== 3 || Encoding/10 == 3)
        InitForKaon();
    else
        InitForPion();
    }

//**********************************************************************************************

G4PomeronCrossSection::G4PomeronCrossSection(const G4Proton * )
{
        InitForNucleon();
}
G4PomeronCrossSection::G4PomeronCrossSection(const G4Neutron * )
{
        InitForNucleon();
}

G4PomeronCrossSection::G4PomeronCrossSection(const G4PionPlus * )
{
        InitForPion();
}
G4PomeronCrossSection::G4PomeronCrossSection(const G4PionMinus * )
{
        InitForPion();
}
G4PomeronCrossSection::G4PomeronCrossSection(const G4PionZero * )
{
        InitForPion();
}

G4PomeronCrossSection::G4PomeronCrossSection(const G4KaonPlus * )
{
        InitForKaon();
}
G4PomeronCrossSection::G4PomeronCrossSection(const G4KaonMinus * )
{
        InitForKaon();
}
G4PomeronCrossSection::G4PomeronCrossSection(const G4KaonZero * )
{
        InitForKaon();
}
G4PomeronCrossSection::G4PomeronCrossSection(const G4KaonZeroLong * )
{
        InitForKaon();
}
G4PomeronCrossSection::G4PomeronCrossSection(const G4KaonZeroShort * )
{
        InitForKaon();
}
G4PomeronCrossSection::G4PomeronCrossSection(const G4Gamma * )
{
        InitForGamma();
}

G4double G4PomeronCrossSection::GetTotalCrossSection(const G4double S)
{
        G4double FZ2= Expand(Z(S)/2);
        return SigP(S) * FZ2;
}


G4double G4PomeronCrossSection::GetElasticCrossSection(const G4double S)
{
        return SigP(S)/pomeron_C *(Expand(Z(S)/2) - Expand(Z(S)));
}

G4double G4PomeronCrossSection::GetDiffractiveCrossSection(const G4double S)
{
        return ( pomeron_C -1) * GetElasticCrossSection(S);
}

G4double G4PomeronCrossSection::GetInelasticCrossSection(const G4double S)
{
        return GetTotalCrossSection(S) - GetElasticCrossSection(S);
}

//-------------------------Probabilities ----------------------------

G4double G4PomeronCrossSection::GetTotalProbability(const G4double S,
                                  const G4double impactsquare)
{
        return 2/pomeron_C*(1-std::exp(-1*Eikonal(S,impactsquare)));
}

G4double G4PomeronCrossSection::GetDiffractiveProbability(const G4double S,
                                  const G4double impactsquare)
{
        return (pomeron_C-1)/pomeron_C * 
                             (GetTotalProbability(S,impactsquare) -
                              GetNondiffractiveProbability(S,impactsquare));
}

G4double G4PomeronCrossSection::GetNondiffractiveProbability(const G4double S,
                                  const G4double impactsquare)
{
        return (1-std::exp(-2*Eikonal(S,impactsquare)))/pomeron_C;
}

G4double G4PomeronCrossSection::GetElasticProbability(const G4double S,
                                  const G4double impactsquare)
{
        return (GetTotalProbability(S,impactsquare) -
                GetInelasticProbability(S,impactsquare));
}

G4double G4PomeronCrossSection::GetInelasticProbability(const G4double S,
                                  const G4double impactsquare)
{
        return GetNondiffractiveProbability(S,impactsquare) +
                GetDiffractiveProbability(S,impactsquare);
}

G4double G4PomeronCrossSection::GetCutPomeronProbability(const G4double S,
                            const G4double impactsquare, const G4int nPomerons)
{
        G4double factorial=G4Pow::GetInstance()->factorial(nPomerons);
         
        return std::exp(-2*Eikonal(S,impactsquare))/pomeron_C*
                std::pow(2*Eikonal(S,impactsquare),nPomerons)/factorial;
}

// ---------------Temporary --- GF
void G4PomeronCrossSection::Setgamma(const G4double agam)
{
        pomeron_Gamma=agam/GeV/GeV;
}


//-----------------  private/Implementation methods

void G4PomeronCrossSection::InitForNucleon()
{
//      pomeron_S=              3.0*GeV*GeV;
        pomeron_S=              2.7*GeV*GeV;
//      pomeron_Gamma=          2.16/GeV/GeV;
//      pomeron_Gamma=          3.96/GeV/GeV;
        pomeron_Gamma=          (2.6+3.96)/GeV/GeV;
        pomeron_C=              1.4;
        pomeron_Rsquare=        3.56/GeV/GeV;
//      pomeron_Alpha=          1.0808;
        pomeron_Alpha=          0.9808;
        pomeron_Alphaprime=     0.25/GeV/GeV;     
        pomeron_Gamma_Hard =    0.0002/GeV/GeV;  // Note! if pomeron_Gamma_Hard != 0 to fit total pp-crosscection
                                                 // pomeron_Gamma_Soft shold be 2.35/GeV/GeV
        pomeron_Alpha_Hard =    1.47;            
}

void G4PomeronCrossSection::InitForPion()
{
        pomeron_S=              1.5*GeV*GeV;
//      pomeron_Gamma=          1.46/GeV/GeV;
        pomeron_Gamma=          2.17/GeV/GeV;
        pomeron_C=              1.6;
        pomeron_Rsquare=        2.36/GeV/GeV;
        pomeron_Alpha=          1.0808;
        pomeron_Alphaprime=     0.25/GeV/GeV;
        pomeron_Gamma_Hard =    0.0002/GeV/GeV;
        pomeron_Alpha_Hard =    1.47;
}

void G4PomeronCrossSection::InitForKaon()
{
        pomeron_S=              2.3*GeV*GeV;
//      pomeron_Gamma=          1.31/GeV/GeV;
        pomeron_Gamma=          1.92/GeV/GeV;
        pomeron_C=              1.8;
        pomeron_Rsquare=        1.96/GeV/GeV;
        pomeron_Alpha=          1.0808;
        pomeron_Alphaprime=     0.25/GeV/GeV;
        pomeron_Gamma_Hard =    0.0002/GeV/GeV;
        pomeron_Alpha_Hard =    1.47;
}

void G4PomeronCrossSection::InitForGamma()
{
        pomeron_S=              1.7*GeV*GeV;
//      pomeron_Gamma=          1.42/GeV/GeV;
        pomeron_Gamma=          2.07/GeV/GeV;
        pomeron_C=              1.7;
        pomeron_Rsquare=        2.16/GeV/GeV;
        pomeron_Alpha=          1.0808;
        pomeron_Alphaprime=     0.25/GeV/GeV;
        pomeron_Gamma_Hard =    0.0002/GeV/GeV;
        pomeron_Alpha_Hard =    1.47;
}

G4double G4PomeronCrossSection::Expand(G4double z)
{

        G4double sum=1.;
        G4double current=1.;
        for (G4int j=2; j<21; j++ )
        {
                current *= -z *(j-1)/sqr(j);
                sum+=current;
        }
        return sum;
}

G4double G4PomeronCrossSection::Power(const G4double S)
{
        return pomeron_Gamma *std::pow(S/pomeron_S, pomeron_Alpha -1);
}

G4double G4PomeronCrossSection::Z(const G4double S)
{
        return 2*pomeron_C * Power(S) / Lambda(S);
}

G4double G4PomeronCrossSection::Lambda(const G4double S)
{
        return pomeron_Rsquare+pomeron_Alphaprime*std::log(S/pomeron_S);
}

G4double G4PomeronCrossSection::SigP(const G4double S)
{
        return 8 * pi * hbarc_squared * Power(S);
}

G4double G4PomeronCrossSection::Eikonal(const G4double S,
                                               const G4double impactsquare)
{
        return Z(S)/2 * std::exp(-impactsquare/(4*Lambda(S)*hbarc_squared));
}
//*************************************************************************************************
G4double G4PomeronCrossSection::PowerSoft(const G4double S)
    {
    return pomeron_Gamma *std::pow(S/pomeron_S, pomeron_Alpha -1);
    }

G4double G4PomeronCrossSection::PowerHard(const G4double S)
    {
    return pomeron_Gamma_Hard*std::pow(S/pomeron_S, pomeron_Alpha_Hard -1);
    }

G4double G4PomeronCrossSection::LambdaSoft(const G4double S)
    {
    return pomeron_Rsquare+pomeron_Alphaprime*std::log(S/pomeron_S);
    }
    
G4double G4PomeronCrossSection::LambdaHard(const G4double /*S*/)
    {
    return pomeron_Rsquare; //+pomeron_Alphaprime*std::log(s/pomeron_S);
    }

G4double G4PomeronCrossSection::Zsoft(const G4double S)
   {
   return 2*pomeron_C*PowerHard(S) / LambdaSoft(S);
   }

G4double G4PomeronCrossSection::Zhard(const G4double S)
   {
   return 2*pomeron_C*PowerHard(S)/LambdaHard(S);
   }

G4double G4PomeronCrossSection::SoftEikonal(G4double S, G4double impactsquare)
    {
    return Zsoft(S)/2*std::exp(-impactsquare/LambdaSoft(S)/hbarc_squared/4);
    }
     
G4double G4PomeronCrossSection::HardEikonal(G4double S, G4double impactsquare)
    {
    return Zhard(S)/2*std::exp(-impactsquare/LambdaHard(S)/hbarc_squared/4);
    }
    
//*************************************************************************************************