Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:13:30 

0001 //-*-C++-*-
0002 //-*-Higgs.cpp-*-
0003 //   Written by James Monk and Andrew Pilkington
0004 /////////////////////////////////////////////////////////////////////////////
0005 #include "GeneratorInterface/ExhumeInterface/interface/Higgs.h"
0006 #include "GeneratorInterface/ExhumeInterface/interface/hdecay.h"
0007 
0008 //#include "CLHEP/HepMC/include/PythiaWrapper6_2.h"
0009 
0010 extern "C" {
0011 extern struct {
0012   int mdcy[3][500], mdme[2][8000];
0013   double brat[8000];
0014   int kfdp[5][8000];
0015 } pydat3_;
0016 }
0017 #define pydat3 pydat3_
0018 
0019 /////////////////////////////////////////////////////////////////////////////
0020 Exhume::Higgs::Higgs(const edm::ParameterSet& pset) : CrossSection(pset) {
0021   std::cout << std::endl << "   = Higgs Subprocess Selected =" << std::endl;
0022 
0023   One = 1.0;
0024   BR = &One;
0025   HiggsMass2 = HiggsMass * HiggsMass;
0026   HiggsWidth = HiggsWidth_();
0027   TotWidth = HiggsWidth;
0028 
0029   SetC();
0030 
0031   /*
0032     double Hm2pGam2 = HiggsMass2 + HiggsWidth * HiggsWidth;
0033 
0034     double SqrtHm2pGam2 = sqrt(Hm2pGam2);
0035     double HmpSqrtHm2pGam2 = HiggsMass + SqrtHm2pGam2;
0036     //C normalises so that integral of |prop|^2 d(m^2) = 1
0037 
0038     C = (2.0/HiggsMass) * M_PI *SqrtHm2pGam2 *
0039       (HiggsWidth*HiggsWidth + HiggsMass * HmpSqrtHm2pGam2) / 
0040       (2.0 * pow(HiggsMass, 4.5) * HiggsWidth * sqrt(2 * HmpSqrtHm2pGam2));
0041 
0042     C = 1.0/sqrt(C);
0043     */
0044   //GGHConst = sqrt((*BR)) * 0.25 * I * LambdaW / M_PI;
0045   FsfTop = 0.25 / (TopMass * TopMass);
0046   FsfBottom = 0.25 / (BottomMass * BottomMass);
0047 
0048   //Gev2fb = 3.88 * pow(10.0,11);
0049   NLOConst = M_PI + 5.5 / M_PI;
0050 
0051   Name = "Higgs";
0052   Partons.resize(1);
0053   Partons[0].id = 25;
0054 }
0055 //////////////////////////////////////////////////////////////////////////////
0056 inline void Exhume::Higgs::SetC() {
0057   double Hm2pGam2 = HiggsMass2 + HiggsWidth * HiggsWidth;
0058   double SqrtHm2pGam2 = sqrt(Hm2pGam2);
0059   double HmpSqrtHm2pGam2 = HiggsMass + SqrtHm2pGam2;
0060 
0061   C = (2.0 / HiggsMass) * M_PI * SqrtHm2pGam2 * (HiggsWidth * HiggsWidth + HiggsMass * HmpSqrtHm2pGam2) /
0062       (2.0 * pow(HiggsMass, 4.5) * HiggsWidth * sqrt(2 * HmpSqrtHm2pGam2));
0063 
0064   C = sqrt(1.0 / C);
0065 
0066   GGHConst = sqrt((*BR)) * 0.25 * I * LambdaW / M_PI;
0067 }
0068 //////////////////////////////////////////////////////////////////////////////
0069 void Exhume::Higgs::SetHiggsMass(const double& HM_) {
0070   HiggsMass = HM_;
0071   HiggsMass2 = HiggsMass * HiggsMass;
0072   HiggsWidth = HiggsWidth_();
0073   TotWidth = HiggsWidth;
0074   HiggsWidth = TotWidth * (*BR);
0075 
0076   SetC();
0077 
0078   /*
0079   double Hm2pGam2 = HiggsMass2 + HiggsWidth * HiggsWidth;
0080   
0081   C = 0.25*M_PI*pow(Hm2pGam2,0.25)*
0082     (2.0*HiggsWidth*HiggsWidth + 3.0*HiggsMass2 - 
0083      HiggsMass*sqrt(Hm2pGam2))/(HiggsWidth*pow(HiggsMass,4.5));
0084   
0085   C = 1.0/sqrt(C);
0086   */
0087   return;
0088 }
0089 /////////////////////////////////////////////////////////////////////////////
0090 double Exhume::Higgs::SubProcess() {
0091   AlphaS_ = AlphaS(SqrtsHat);
0092 
0093   double ModAmp = abs(GluGlu2HiggsAmp() * Propagator() * C);
0094 
0095   return ((1.0 + AlphaS_ * NLOConst) * (Gev2fb * ModAmp * ModAmp * M_PI * InvsHat2));
0096 }
0097 /////////////////////////////////////////////////////////////////////////////
0098 double Exhume::Higgs::HiggsWidth_() {
0099   //  std::cout<<std::endl<<
0100   //"  Using HDECAY to calculate S.M. Higgs Width"<<std::endl;
0101 
0102   flag_.ihiggs = 0;
0103   flag_.nnlo = 0;
0104   flag_.ipole = 1;
0105   oldfash_.nfgg = 5;
0106   onshell_.ionsh = 0;
0107   onshell_.ionwz = 0;
0108   param_.gf = 1.16639e-5;
0109   param_.alph = AlphaEw;  //1.0/137.0;
0110   param_.amtau = 1.7771;  //tau_mass;//1.7771;
0111   param_.ammuon = MuonMass;
0112   param_.amz = ZMass;  //91.187;
0113   param_.amw = WMass;  //80.33;
0114   ckmpar_.vus = 0.2205;
0115   ckmpar_.vcb = 0.04;
0116   ckmpar_.vub = 0.08 * ckmpar_.vcb;
0117   masses_.ams = StrangeMass;
0118   masses_.amc = 1.42;        //charm_mass;//1.42;
0119   masses_.amb = BottomMass;  //4.62;
0120   masses_.amt = TopMass;     //175.0;
0121   strange_.amsb = masses_.ams;
0122   als_.amc0 = masses_.amc;
0123   als_.amb0 = masses_.amb;
0124   als_.amt0 = masses_.amt;
0125   double acc = 1.0e-8;
0126   int nloop = 2;
0127   double alsmz = 0.118;
0128 
0129   als_.xlambda = xitla_(&nloop, &alsmz, &acc);
0130   als_.n0 = 5;
0131   alsini_(&acc);
0132   wzwdth_.gamw = 2.08;
0133   wzwdth_.gamz = 2.49;
0134 
0135   int nber = 18;
0136   bernini_(&nber);
0137   hmass_.amsm = HiggsMass;
0138   double GfAmt = param_.gf * masses_.amt;
0139   double Mw_Mt2 = (param_.amw * param_.amw) / (masses_.amt * masses_.amt);
0140   double one_Mw_Mt2 = 1.0 - Mw_Mt2;
0141   wzwdth_.gamt0 = GfAmt * GfAmt * GfAmt * 0.125 / sqrt(2.0) / M_PI * one_Mw_Mt2 * one_Mw_Mt2 * (1 + 2 * Mw_Mt2);
0142   wzwdth_.gamt1 = wzwdth_.gamt0;
0143   //double tgbet = 1.0;
0144   //hdec_(&tgbet);
0145   hdec_();
0146 
0147   //std::cout<<std::endl<<"  S.M. Higgs width = "<<widthsm_.smwdth<<std::endl;
0148 
0149   return (widthsm_.smwdth);
0150 }
0151 /////////////////////////////////////////////////////////////////////////////
0152 void Exhume::Higgs::SetPartons() {
0153   Partons[0].p = CentralVector;
0154 
0155   return;
0156 }
0157 /////////////////////////////////////////////////////////////////////////////
0158 inline double Exhume::Higgs::SubParameterWeight() { return (1.0); }
0159 /////////////////////////////////////////////////////////////////////////////
0160 inline double Exhume::Higgs::SubParameterRange() { return (1.0); }
0161 /////////////////////////////////////////////////////////////////////////////
0162 void Exhume::Higgs::MaximiseSubParameters() { return; }
0163 /////////////////////////////////////////////////////////////////////////////
0164 void Exhume::Higgs::SetSubParameters() { return; }
0165 /////////////////////////////////////////////////////////////////////////////
0166 void Exhume::Higgs::SetHiggsDecay(const int& id_) {
0167   pydat3.mdme[0][209] = 0;
0168   pydat3.mdme[0][210] = 0;
0169   pydat3.mdme[0][211] = 0;
0170   pydat3.mdme[0][212] = 0;
0171   pydat3.mdme[0][213] = 0;
0172   pydat3.mdme[0][214] = 0;
0173   pydat3.mdme[0][215] = 0;
0174   pydat3.mdme[0][216] = 0;
0175   pydat3.mdme[0][217] = 0;
0176   pydat3.mdme[0][218] = 0;
0177   pydat3.mdme[0][219] = 0;
0178   pydat3.mdme[0][220] = 0;
0179   pydat3.mdme[0][221] = 0;
0180   pydat3.mdme[0][222] = 0;
0181   pydat3.mdme[0][223] = 0;
0182   pydat3.mdme[0][224] = 0;
0183   pydat3.mdme[0][225] = 0;
0184 
0185   //double BR;
0186 
0187   switch (id_) {
0188     case 3:
0189       //strange
0190       pydat3.mdme[0][211] = 1;
0191       BR = &widthsm_.smbrs;
0192       break;
0193     case 4:
0194       //charm
0195       pydat3.mdme[0][212] = 1;
0196       BR = &widthsm_.smbrc;
0197       break;
0198     case 5:
0199       //bottom
0200       pydat3.mdme[0][213] = 1;
0201       BR = &widthsm_.smbrb;
0202       break;
0203     case 6:
0204       //top
0205       pydat3.mdme[0][214] = 1;
0206       BR = &widthsm_.smbrt;
0207       break;
0208     case 13:
0209       //muon
0210       pydat3.mdme[0][218] = 1;
0211       BR = &widthsm_.smbrm;
0212       break;
0213     case 15:
0214       //tau
0215       pydat3.mdme[0][219] = 1;
0216       BR = &widthsm_.smbrl;
0217       break;
0218     case 21:
0219       //gluon
0220       pydat3.mdme[0][221] = 1;
0221       BR = &widthsm_.smbrg;
0222       break;
0223     case 22:
0224       //photon
0225       pydat3.mdme[0][222] = 1;
0226       BR = &widthsm_.smbrga;
0227       break;
0228     case 23:
0229       //Z
0230       pydat3.mdme[0][224] = 1;
0231       BR = &widthsm_.smbrz;
0232       break;
0233     case 24:
0234       //W
0235       pydat3.mdme[0][225] = 1;
0236       BR = &widthsm_.smbrw;
0237       break;
0238 
0239     default:
0240       std::cout << "   Unrecognised PDG code for Higgs Decay" << std::endl;
0241       for (int alldecay = 209; alldecay < 226; alldecay++) {
0242         pydat3.mdme[0][alldecay] = 1;
0243       }
0244       BR = &One;  //1.0;
0245       break;
0246   }
0247 
0248   //HiggsWidth = TotWidth * (*BR);
0249 
0250   SetC();
0251 
0252   /*
0253 
0254   double Hm2pGam2 = HiggsMass2 + HiggsWidth * HiggsWidth;
0255   double SqrtHm2pGam2 = sqrt(Hm2pGam2);
0256   double HmpSqrtHm2pGam2 = HiggsMass + SqrtHm2pGam2;
0257 
0258   C = (2.0/HiggsMass) * M_PI *SqrtHm2pGam2 *
0259     (HiggsWidth*HiggsWidth + HiggsMass * HmpSqrtHm2pGam2) / 
0260     (2.0 * pow(HiggsMass, 4.5) * HiggsWidth * sqrt(2 * HmpSqrtHm2pGam2));
0261 
0262   C = sqrt((*BR)/C);
0263   */
0264   return;
0265 }
0266 /////////////////////////////////////////////////////////////////////////////
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!