CMSSW_4_4_3_patch1/src/GeneratorInterface/ExhumeInterface/interface/CrossSection.h

Go to the documentation of this file.

//-*-C++-*-
//-*-CrossSection.h-*-
//   Written by James Monk and Andrew Pilkington

#ifndef CROSSSECTION_HH
#define CROSSSECTION_HH

#include <cmath>
#include <complex>
#include <cstdlib>
#include <iostream>
#include <map>
#include <string>
#include <utility>
#include <vector>

//#include "CLHEP/config/CLHEP.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Vector/ThreeVector.h"

#include "GeneratorInterface/ExhumeInterface/interface/I.h"
#include "GeneratorInterface/ExhumeInterface/interface/PI.h"
//#include "GeneratorInterface/ExhumeInterface/interface/PythiaRecord.h"
#include "GeneratorInterface/ExhumeInterface/interface/Particle.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

//#include "CLHEP/HepMC/include/PythiaWrapper6_2.h"
//#include "CLHEP/HepMC/ConvertHEPEVT.h"
//#include "CLHEP/HepMC/CBhepevt.h"

namespace CLHEP {
class HepRandomEngine;
}

namespace Exhume{
  
  typedef std::pair<double,double> fEntry;
  typedef std::pair<const char*,char*> PCharPair;
  typedef std::pair<const char*,const void*> PConstVoidPair;
  
  class CrossSection{

  public:

    CrossSection(const edm::ParameterSet&);
    virtual ~CrossSection();
    virtual void MaximiseSubParameters()=0;
    virtual void SetPartons()=0;
    virtual void SetSubParameters()=0;
    virtual double SubParameterRange()=0;
    virtual double SubParameterWeight()=0;

    double AlphaS(const double&);

    inline void SetRandomEngine(CLHEP::HepRandomEngine* engine){randomEngine = engine;}

    inline double GetRg(const double &x_, const double &Qt){

      return(Rg_(x_, Qt) );
    };

    //.....
    inline double Differential(){

      if(x1>1.0 || x2 > 1.0){
        return(0.0);
      }

      //return( SubProcess() );

      //Factor of 2/sqrt(sHat) because we integrate dM not dln(M^2)
      return( 2.0 * InvSqrtsHat * Lumi() * SubProcess() * exp(B*(t1 + t2))); 
    };
    //.....

    inline double GetB(){
      return(B);
    };

    inline double GetRoot_s(){
      return(root_s);
    };

    inline std::map<double,double> Getfg2Map(){
      return(fg2Map);
    };

    inline CLHEP::HepLorentzVector GetProton1(){
      return(Proton1);
    };

    inline CLHEP::HepLorentzVector GetProton2(){
      return(Proton2);
    };

    inline std::vector<Particle> GetPartons(){
      return(Partons);
    };

    inline double Getx1(){
      return(x1);
    };
    
    inline double Getx2(){
      return(x2);
    };

    inline double Gett1(){
      return(t1);
    };

    inline double Gett2(){
      return(t2);
    };
    
    inline double GetsHat(){
      return(sHat);
    };

    inline double GetSqrtsHat(){
      return(SqrtsHat);
    };

    inline double GetEta(){
      return(y);
    };

    inline double GetPhi1(){
      return(Phi1);
    };

    inline double GetPhi2(){
      return(Phi2);
    };

    void Hadronise();
    void SetKinematics(const double&, const double&,
                       const double&, const double&, 
                       const double&, const double&);

    inline std::string GetName(){
      return(Name);
    };

  private:

    void AlphaSInit();

    double Fg1Fg2(const double&, const double&, const double&);
    double Fg1Fg2(const int&, const double&, const double&);

    inline double Fg_Qt2(const double &Qt2_, 
                         const double &x_, const double &xp_){

      double grad = 5.0 * (Txg(1.1 * Qt2_, x_) - Txg(0.9 * Qt2_, x_) ) / Qt2_;

      //protect from Qt2 == 0

      if(grad!=grad){

        return(0.0);
      }
      return(grad);    
    };
    
    template <typename T_>
    inline void insert(const std::string _name_, const T_ _x_){
      PMap.insert
        (std::pair<std::string,PConstVoidPair>
         (_name_,PConstVoidPair(typeid(_x_).name(),_x_)));
    }

    virtual double SubProcess() = 0;

    void LumiInit();
    double Lumi();
    double Lumi_();
    void NoMem();
    void   ReadCard(const std::string&);
    double Splitting(const double&);
    double T(const double&);
    double TFast(const double&);

    bool TBegin, TInterpolate;
    std::map<double, std::map<double, double> > TMap2d;
    std::map<double, std::map<double, double> >::iterator TMjuHigh, TMjuLow;

    double Txg(const double&, const double&);

    double Rg1Rg2(const double&);
    double Rg_(const double&, double);

    std::map<double, std::map<double, double> > RgMap2d;
    std::map<double, std::map<double, double> >::iterator RgHigh[2], RgLow[2];
    bool RgInterpolate[2], RgBegin;

    double B,
      Freeze,
      LambdaQCD,
      Rg,
      Survive;

    double PDF;
      
    double ASFreeze, ASConst;
    //.............................
    //Used for Lumi function:
    
    double MinQt2, MidQt2, InvMidQt2, InvMidQt4, InvMaxQt2,
      LumSimpsIncr;
    unsigned int LumAccuracy, LumNStart;
    int LumNSimps, LumNSimps_1;
    double *LumSimpsFunc, *_Qt2, *_Qt, *_KtLow, *_KtHigh, 
      *_AlphaS, *_CfAs_2PIRg, *_NcAs_2PI;
    double LumConst;
    double Inv2PI, Nc_2PI, Cf_2PIRg;

    std::map<std::string,PConstVoidPair> PMap;
    std::map<const char*,char*> TypeMap;
    std::map<double, double> fg2Map;

    //...............................
    //Used for T:
    double TConst;
    int Tn, Tn_1;
    double *TFunc;

    //...............................
    //Used for Splitting:
    double Inv3;

    int Proton1Id, Proton2Id;

  protected:

    std::string Name;

    std::complex<double> F0(const double&);
    std::complex<double> f(const double&);
    std::complex<double> Fsf(const double&);

    //PPhi is azimuthal angle between protons.  
    //InvSqrtsHat = 1/sHat
    //Others are obvious.
    double x1, x1p, x2, x2p, t1, t2, sHat,  SqrtsHat,
      sHat2, InvsHat, InvsHat2, InvSqrtsHat, y, PPhi, Phi1, Phi2,
      Mju2, Mju, LnMju2,
      Pt1, Pt2, Pt1DotPt2,
      x1x2, ey;

    double AlphaEw,
      gw,
      HiggsVev,
      LambdaW;

    double BottomMass,
      CharmMass,
      StrangeMass,
      TopMass;

    double MuonMass,
      TauMass;

    double HiggsMass,
      WMass,
      ZMass;

    int FNAL_or_LHC;
    double root_s,
      s, Invs;

    CLHEP::HepLorentzVector CentralVector;
    CLHEP::HepLorentzVector Proton1, Proton2, P1In, P2In;
    /*
    struct Particle {
      CLHEP::HepLorentzVector p;
      int id;
    };
    */

    std::vector<Particle> Partons;
    double Invsx1x2, InvV1MinusV2;

    double Gev2fb;

    std::string lhapdfSetPath_;

    CLHEP::HepRandomEngine* randomEngine;
  };
}

#endif