ExhumeInterface/interface/CrossSection.h

0001 //-*-C++-*-
0002 //-*-CrossSection.h-*-
0003 //   Written by James Monk and Andrew Pilkington
0004 /////////////////////////////////////////////////////////////////////////////
0005
0006 #ifndef CROSSSECTION_HH
0007 #define CROSSSECTION_HH
0008
0009 #include <cmath>
0010 #include <complex>
0011 #include <cstdlib>
0012 #include <iostream>
0013 #include <map>
0014 #include <string>
0015 #include <utility>
0016 #include <vector>
0017
0018 //#include "CLHEP/config/CLHEP.h"
0019 #include "CLHEP/Vector/LorentzVector.h"
0020 #include "CLHEP/Vector/ThreeVector.h"
0021
0022 #include "GeneratorInterface/ExhumeInterface/interface/I.h"
0023 //#include "GeneratorInterface/ExhumeInterface/interface/PythiaRecord.h"
0024 #include "GeneratorInterface/ExhumeInterface/interface/Particle.h"
0025
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027
0028 //#include "CLHEP/HepMC/include/PythiaWrapper6_2.h"
0029 //#include "CLHEP/HepMC/ConvertHEPEVT.h"
0030 //#include "CLHEP/HepMC/CBhepevt.h"
0031
0032 namespace CLHEP {
0033   class HepRandomEngine;
0034 }
0035
0036 /////////////////////////////////////////////////////////////////////////////
0037 namespace Exhume {
0038
0039   typedef std::pair<double, double> fEntry;
0040   typedef std::pair<const char *, char *> PCharPair;
0041   typedef std::pair<const char *, const void *> PConstVoidPair;
0042
0043   class CrossSection {
0044   public:
0045     CrossSection(const edm::ParameterSet &);
0046     virtual ~CrossSection();
0047     virtual void MaximiseSubParameters() = 0;
0048     virtual void SetPartons() = 0;
0049     virtual void SetSubParameters() = 0;
0050     virtual double SubParameterRange() = 0;
0051     virtual double SubParameterWeight() = 0;
0052
0053     double AlphaS(const double &);
0054
0055     inline void SetRandomEngine(CLHEP::HepRandomEngine *engine) { randomEngine = engine; }
0056
0057     inline double GetRg(const double &x_, const double &Qt) { return (Rg_(x_, Qt)); };
0058
0059     //.....
0060     inline double Differential() {
0061       if (x1 > 1.0 || x2 > 1.0) {
0062         return (0.0);
0063       }
0064
0065       //return( SubProcess() );
0066
0067       //Factor of 2/sqrt(sHat) because we integrate dM not dln(M^2)
0068       return (2.0 * InvSqrtsHat * Lumi() * SubProcess() * exp(B * (t1 + t2)));
0069     };
0070     //.....
0071
0072     inline double GetB() { return (B); };
0073
0074     inline double GetRoot_s() { return (root_s); };
0075
0076     inline std::map<double, double> Getfg2Map() { return (fg2Map); };
0077
0078     inline CLHEP::HepLorentzVector GetProton1() { return (Proton1); };
0079
0080     inline CLHEP::HepLorentzVector GetProton2() { return (Proton2); };
0081
0082     inline std::vector<Particle> GetPartons() { return (Partons); };
0083
0084     inline double Getx1() { return (x1); };
0085
0086     inline double Getx2() { return (x2); };
0087
0088     inline double Gett1() { return (t1); };
0089
0090     inline double Gett2() { return (t2); };
0091
0092     inline double GetsHat() { return (sHat); };
0093
0094     inline double GetSqrtsHat() { return (SqrtsHat); };
0095
0096     inline double GetEta() { return (y); };
0097
0098     inline double GetPhi1() { return (Phi1); };
0099
0100     inline double GetPhi2() { return (Phi2); };
0101
0102     void Hadronise();
0103     void SetKinematics(const double &, const double &, const double &, const double &, const double &, const double &);
0104
0105     inline std::string GetName() { return (Name); };
0106
0107   private:
0108     void AlphaSInit();
0109
0110     double Fg1Fg2(const double &, const double &, const double &);
0111     double Fg1Fg2(const int &, const double &, const double &);
0112
0113     inline double Fg_Qt2(const double &Qt2_, const double &x_, const double &xp_) {
0114       double grad = 5.0 * (Txg(1.1 * Qt2_, x_) - Txg(0.9 * Qt2_, x_)) / Qt2_;
0115
0116       //protect from Qt2 == 0
0117
0118       if (grad != grad) {
0119         return (0.0);
0120       }
0121       return (grad);
0122     };
0123
0124     template <typename T_>
0125     inline void insert(const std::string _name_, const T_ _x_) {
0126       PMap.insert(std::pair<std::string, PConstVoidPair>(_name_, PConstVoidPair(typeid(_x_).name(), _x_)));
0127     }
0128
0129     virtual double SubProcess() = 0;
0130
0131     void LumiInit();
0132     double Lumi();
0133     double Lumi_();
0134     void NoMem();
0135     void ReadCard(const std::string &);
0136     double Splitting(const double &);
0137     double T(const double &);
0138     double TFast(const double &);
0139
0140     bool TBegin, TInterpolate;
0141     std::map<double, std::map<double, double> > TMap2d;
0142     std::map<double, std::map<double, double> >::iterator TMjuHigh, TMjuLow;
0143
0144     double Txg(const double &, const double &);
0145
0146     double Rg1Rg2(const double &);
0147     double Rg_(const double &, double);
0148
0149     std::map<double, std::map<double, double> > RgMap2d;
0150     std::map<double, std::map<double, double> >::iterator RgHigh[2], RgLow[2];
0151     bool RgInterpolate[2], RgBegin;
0152
0153     double B, Freeze, LambdaQCD, Rg, Survive;
0154
0155     double PDF;
0156
0157     double ASFreeze, ASConst;
0158     //.............................
0159     //Used for Lumi function:
0160
0161     double MinQt2, MidQt2, InvMidQt2, InvMidQt4, InvMaxQt2, LumSimpsIncr;
0162     unsigned int LumAccuracy, LumNStart;
0163     int LumNSimps, LumNSimps_1;
0164     double *LumSimpsFunc, *_Qt2, *_Qt, *_KtLow, *_KtHigh, *_AlphaS, *_CfAs_2PIRg, *_NcAs_2PI;
0165     double LumConst;
0166     double Inv2PI, Nc_2PI, Cf_2PIRg;
0167
0168     std::map<std::string, PConstVoidPair> PMap;
0169     std::map<const char *, char *> TypeMap;
0170     std::map<double, double> fg2Map;
0171
0172     //...............................
0173     //Used for T:
0174     double TConst;
0175     int Tn, Tn_1;
0176     double *TFunc;
0177
0178     //...............................
0179     //Used for Splitting:
0180     double Inv3;
0181
0182     int Proton1Id, Proton2Id;
0183
0184   protected:
0185     std::string Name;
0186
0187     std::complex<double> F0(const double &);
0188     std::complex<double> f(const double &);
0189     std::complex<double> Fsf(const double &);
0190
0191     //PPhi is azimuthal angle between protons.
0192     //InvSqrtsHat = 1/sHat
0193     //Others are obvious.
0194     double x1, x1p, x2, x2p, t1, t2, sHat, SqrtsHat, sHat2, InvsHat, InvsHat2, InvSqrtsHat, y, PPhi, Phi1, Phi2, Mju2,
0195         Mju, LnMju2, Pt1, Pt2, Pt1DotPt2, x1x2, ey;
0196
0197     double AlphaEw, gw, HiggsVev, LambdaW;
0198
0199     double BottomMass, CharmMass, StrangeMass, TopMass;
0200
0201     double MuonMass, TauMass;
0202
0203     double HiggsMass, WMass, ZMass;
0204
0205     int FNAL_or_LHC;
0206     double root_s, s, Invs;
0207
0208     CLHEP::HepLorentzVector CentralVector;
0209     CLHEP::HepLorentzVector Proton1, Proton2, P1In, P2In;
0210     /*
0211     struct Particle {
0212       CLHEP::HepLorentzVector p;
0213       int id;
0214     };
0215     */
0216
0217     std::vector<Particle> Partons;
0218     double Invsx1x2, InvV1MinusV2;
0219
0220     double Gev2fb;
0221
0222     std::string lhapdfSetPath_;
0223
0224     CLHEP::HepRandomEngine *randomEngine;
0225   };
0226 }  // namespace Exhume
0227
0228 #endif
0229 ////////////////////////////////////////////////////////////////////////////