File indexing completed on 2024-04-06 12:13:29
0001
0002
0003
0004
0005
0006 #include "GeneratorInterface/ExhumeInterface/interface/CrossSection.h"
0007 #include "GeneratorInterface/ExhumeInterface/interface/PythiaRecord.h"
0008
0009 #include "HepMC/PythiaWrapper6_4.h"
0010
0011
0012 #include <cstdio>
0013 #include <memory>
0014 #include <cmath>
0015
0016
0017 double dsimps_(double *, double *, double *, int *);
0018
0019
0020
0021 #define my_pdfset my_pdfset_
0022 extern "C" {
0023 double my_pdfset(double &);
0024 }
0025
0026 #define structm structm_
0027 extern "C" {
0028 void structm(
0029 double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *);
0030 }
0031
0032 #define my_pythia_init my_pythia_init_
0033 extern "C" {
0034 void my_pythia_init();
0035 }
0036
0037 extern "C" {
0038 void setpdfpath_(char *);
0039 void setlhaparm_(char *);
0040 }
0041
0042 #define setpdfpath setpdfpath_
0043 #define setlhaparm setlhaparm_
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061 Exhume::CrossSection::CrossSection(const edm::ParameterSet &pset) {
0062 std::cout << std::endl << " ........................................................................." << std::endl;
0063 std::cout << std::endl << " ExHuME version 1.3.2" << std::endl;
0064 std::cout << " Last Updated 17 October 2005" << std::endl << std::endl;
0065
0066 std::cout << " Authors:\tJames Monk\t\tjmonk@hep.man.ac.uk" << std::endl;
0067 std::cout << "\t\tAndrew Pilkington\tpilko@hep.man.ac.uk" << std::endl;
0068
0069 std::cout << std::endl << " ........................................................................." << std::endl;
0070 std::cout << std::endl << " = Initialising CrossSection =" << std::endl;
0071
0072 edm::ParameterSet paramsPSet = pset.getParameter<edm::ParameterSet>("ExhumeParameters");
0073 B = paramsPSet.getParameter<double>("B");
0074 LambdaQCD = paramsPSet.getParameter<double>("LambdaQCD");
0075 Rg = paramsPSet.getParameter<double>("Rg");
0076 Survive = paramsPSet.getParameter<double>("Survive");
0077 PDF = paramsPSet.getParameter<double>("PDF");
0078 MinQt2 = paramsPSet.getParameter<double>("MinQt2");
0079 AlphaEw = paramsPSet.getParameter<double>("AlphaEw");
0080 HiggsVev = paramsPSet.getParameter<double>("HiggsVev");
0081 BottomMass = paramsPSet.getParameter<double>("BottomMass");
0082 CharmMass = paramsPSet.getParameter<double>("CharmMass");
0083 StrangeMass = paramsPSet.getParameter<double>("StrangeMass");
0084 TopMass = paramsPSet.getParameter<double>("TopMass");
0085 MuonMass = paramsPSet.getParameter<double>("MuonMass");
0086 TauMass = paramsPSet.getParameter<double>("TauMass");
0087 HiggsMass = paramsPSet.getParameter<double>("HiggsMass");
0088 WMass = paramsPSet.getParameter<double>("WMass");
0089 ZMass = paramsPSet.getParameter<double>("ZMass");
0090
0091 FNAL_or_LHC = -1;
0092 root_s = pset.getParameter<double>("comEnergy");
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102 insert("AlphaEw", &AlphaEw);
0103 insert("WMass", &WMass);
0104 insert("ZMass", &ZMass);
0105 insert("HiggsMass", &HiggsMass);
0106 insert("HiggsVev", &HiggsVev);
0107 insert("MinQt2", &MinQt2);
0108
0109 insert("TopMass", &TopMass);
0110 insert("BottomMass", &BottomMass);
0111 insert("CharmMass", &CharmMass);
0112 insert("StrangeMass", &StrangeMass);
0113 insert("TauMass", &TauMass);
0114 insert("MuonMass", &MuonMass);
0115 insert("LambdaQCD", &LambdaQCD);
0116 insert("Freeze", &Freeze);
0117 insert("B", &B);
0118 insert("gw", &gw);
0119 insert("LambdaW", &LambdaW);
0120 insert("FNAL_or_LHC", &FNAL_or_LHC);
0121 insert("s", &s);
0122 insert("root_s", &root_s);
0123 insert("Rg", &Rg);
0124 insert("Survive", &Survive);
0125 insert("PDF", &PDF);
0126
0127
0128 Freeze = sqrt(MinQt2);
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140 s = root_s * root_s;
0141 Invs = 1.0 / s;
0142
0143 gw = sqrt(4.0 * M_PI * AlphaEw / (1.0 - WMass * WMass / (ZMass * ZMass)));
0144
0145 LambdaW = 0.5 * gw / WMass;
0146
0147 Gev2fb = 3.88 * pow(10.0, 11);
0148
0149
0150
0151 for (std::map<std::string, PConstVoidPair>::iterator ii = PMap.begin(); ii != PMap.end(); ii++) {
0152 if (typeid(double *).name() == (ii->second).first) {
0153 if (strlen((ii->first).c_str()) < 6) {
0154 printf(" %s\t\t\t%17g\n", (ii->first).c_str(), *(double *)((ii->second).second));
0155 } else {
0156 printf(" %s\t\t%17g\n", (ii->first).c_str(), *(double *)((ii->second).second));
0157 }
0158 }
0159 }
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178 std::cout << std::endl
0179 << " ........................................................................." << std::endl
0180 << std::endl
0181 << " = Initialising PDFs =" << std::endl;
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192 my_pdfset(PDF);
0193 std::cout << std::endl << " ........................................................................." << std::endl;
0194
0195
0196
0197 my_pythia_init();
0198
0199
0200 pyinre();
0201 Proton1Id = 2212;
0202
0203 Proton2Id = 2212;
0204
0205 pydat2.pmas[0][24] = HiggsMass;
0206 pydat2.pmas[0][22] = ZMass;
0207 pydat2.pmas[0][23] = WMass;
0208 pydat2.pmas[0][5] = TopMass;
0209 pydat2.pmas[0][12] = MuonMass;
0210 pydat2.pmas[0][14] = TauMass;
0211 pydat2.pmas[0][4] = BottomMass;
0212 pydat2.pmas[0][3] = CharmMass;
0213 pydat2.pmas[0][2] = StrangeMass;
0214 pydat2.pmas[0][1] = 0.0;
0215 pydat2.pmas[0][0] = 0.0;
0216
0217 AlphaSInit();
0218 LumiInit();
0219 }
0220
0221 Exhume::CrossSection::~CrossSection() {
0222 free(LumSimpsFunc);
0223 free(_Qt2);
0224 free(_Qt);
0225 free(_KtLow);
0226 free(_KtHigh);
0227 free(_AlphaS);
0228 free(_CfAs_2PIRg);
0229 free(_NcAs_2PI);
0230 free(TFunc);
0231 }
0232
0233 double Exhume::CrossSection::AlphaS(const double &scale_) {
0234 if (scale_ > Freeze) {
0235 return (ASConst / log(scale_ / LambdaQCD));
0236 } else {
0237 return (ASFreeze);
0238 }
0239 return (0.0);
0240 }
0241
0242 double Exhume::CrossSection::Rg1Rg2(const double &scale_) {
0243 double Rg1 = 0.;
0244 double Rg2 = 0.;
0245
0246 if (RgBegin) {
0247 RgBegin = false;
0248
0249 RgHigh[0] = RgMap2d.upper_bound(x1);
0250 RgHigh[1] = RgMap2d.upper_bound(x2);
0251
0252 RgLow[0] = RgHigh[0];
0253 RgLow[0]--;
0254 RgLow[1] = RgHigh[1];
0255 RgLow[1]--;
0256
0257
0258
0259 if (fabs(RgLow[0]->first - RgHigh[0]->first) > 0.05 * x1 || RgHigh[0] == RgMap2d.end() ||
0260 RgHigh[0] == RgMap2d.begin()) {
0261 std::map<double, double> Rg1Temp;
0262 RgMap2d.insert(std::pair<double, std::map<double, double> >(x1, Rg1Temp));
0263 RgHigh[0] = RgMap2d.upper_bound(0.9999 * x1);
0264 RgInterpolate[0] = false;
0265 }
0266
0267 if (fabs(RgLow[1]->first - RgHigh[1]->first) > 0.05 * x2 || RgHigh[1] == RgMap2d.end() ||
0268 RgHigh[1] == RgMap2d.begin()) {
0269 std::map<double, double> Rg2Temp;
0270 RgMap2d.insert(std::pair<double, std::map<double, double> >(x2, Rg2Temp));
0271 RgHigh[1] = RgMap2d.upper_bound(0.9999 * x2);
0272 RgInterpolate[1] = false;
0273 }
0274
0275
0276
0277
0278
0279 }
0280
0281 if (!RgInterpolate[0]) {
0282 Rg1 = Rg_(x1, scale_);
0283 (RgHigh[0]->second).insert(std::pair<double, double>(scale_, Rg1));
0284 }
0285
0286 if (!RgInterpolate[1]) {
0287 Rg2 = Rg_(x2, scale_);
0288 (RgHigh[1]->second).insert(std::pair<double, double>(scale_, Rg2));
0289 }
0290
0291 if (RgInterpolate[0]) {
0292 double RgxHigh, RgxLow;
0293
0294 std::map<double, double>::iterator high_, low_;
0295 high_ = (RgHigh[0]->second).upper_bound(scale_);
0296 bool CalculateRg = true;
0297 if (high_ != (RgHigh[0]->second).end() && high_ != (RgHigh[0]->second).begin()) {
0298 low_ = high_;
0299 low_--;
0300 if (high_->first - low_->first < 0.05 * scale_)
0301 CalculateRg = false;
0302 }
0303
0304 if (CalculateRg) {
0305 RgxHigh = Rg_(RgHigh[0]->first, scale_);
0306 (RgHigh[0]->second).insert(std::pair<double, double>(scale_, RgxHigh));
0307 } else {
0308 RgxHigh = low_->second + (scale_ - low_->first) * (high_->second - low_->second) / (high_->first - low_->first);
0309 }
0310
0311 CalculateRg = true;
0312
0313 high_ = (RgLow[0]->second).upper_bound(scale_);
0314 if (high_ != (RgLow[0]->second).end() && high_ != (RgLow[0]->second).begin()) {
0315 low_ = high_;
0316 low_--;
0317 if (high_->first - low_->first < 0.05 * scale_)
0318 CalculateRg = false;
0319 }
0320
0321 if (CalculateRg) {
0322 RgxLow = Rg_(RgLow[0]->first, scale_);
0323 (RgLow[0]->second).insert(std::pair<double, double>(scale_, RgxLow));
0324 } else {
0325 RgxLow = low_->second + (scale_ - low_->first) * (high_->second - low_->second) / (high_->first - low_->first);
0326 }
0327
0328 Rg1 = RgxLow + (x1 - RgLow[0]->first) * (RgxHigh - RgxLow) / (RgHigh[0]->first - RgLow[0]->first);
0329 }
0330
0331 if (RgInterpolate[1]) {
0332 double RgxHigh, RgxLow;
0333
0334 std::map<double, double>::iterator high_, low_;
0335 high_ = (RgHigh[1]->second).upper_bound(scale_);
0336
0337 bool CalculateRg = true;
0338 if (high_ != (RgHigh[1]->second).end() && high_ != (RgHigh[1]->second).begin()) {
0339 low_ = high_;
0340 low_--;
0341 if (high_->first - low_->first < 0.05 * scale_)
0342 CalculateRg = false;
0343 }
0344
0345 if (CalculateRg) {
0346 RgxHigh = Rg_(RgHigh[1]->first, scale_);
0347 (RgHigh[1]->second).insert(std::pair<double, double>(scale_, RgxHigh));
0348 } else {
0349 RgxHigh = low_->second + (scale_ - low_->first) * (high_->second - low_->second) / (high_->first - low_->first);
0350 }
0351
0352 CalculateRg = true;
0353 high_ = (RgLow[1]->second).upper_bound(scale_);
0354
0355 if (high_ != (RgLow[1]->second).end() && high_ != (RgLow[1]->second).begin()) {
0356 low_ = high_;
0357 low_--;
0358 if (high_->first - low_->first < 0.05 * scale_)
0359 CalculateRg = false;
0360 }
0361
0362 if (CalculateRg) {
0363 RgxLow = Rg_(RgLow[1]->first, scale_);
0364 (RgLow[1]->second).insert(std::pair<double, double>(scale_, RgxLow));
0365 } else {
0366 RgxLow = low_->second + (scale_ - low_->first) * (high_->second - low_->second) / (high_->first - low_->first);
0367 }
0368
0369 Rg2 = RgxLow + (x2 - RgLow[1]->first) * (RgxHigh - RgxLow) / (RgHigh[1]->first - RgLow[1]->first);
0370 }
0371
0372 return (Rg1 * Rg2);
0373 }
0374
0375 double Exhume::CrossSection::Rg_(const double &x_, double scale_) {
0376 if (x_ < 0.002 && scale_ < 1.5625)
0377 return (1.0);
0378
0379 double upv, dnv, usea, dsea, str, chm, bot, top, gl1, gl2;
0380
0381 double x1_, x2_;
0382
0383 x1_ = 0.995 * x_;
0384 x2_ = 1.005 * x_;
0385
0386 if (x2_ > 1.0) {
0387 x1_ = 0.99 * x_;
0388 x2_ = x_;
0389 }
0390
0391 structm(&x1_, &scale_, &upv, &dnv, &usea, &dsea, &str, &chm, &bot, &top, &gl1);
0392 structm(&x2_, &scale_, &upv, &dnv, &usea, &dsea, &str, &chm, &bot, &top, &gl2);
0393
0394 double dlam = (gl1 - gl2) * 100.0 / gl1;
0395 if (dlam > 1.0)
0396 dlam = 1.0;
0397
0398 if (dlam < 0.0)
0399 dlam = 0.0;
0400
0401 return (1.0 + dlam * (0.82 + 0.56 * dlam));
0402 }
0403
0404 void Exhume::CrossSection::SetKinematics(const double &SqrtsHat_,
0405 const double &y_,
0406 const double &t1_,
0407 const double &t2_,
0408 const double &Phi1_,
0409 const double &Phi2_) {
0410 RgBegin = true;
0411 RgInterpolate[0] = true;
0412 RgInterpolate[1] = true;
0413
0414 TBegin = true;
0415 TInterpolate = true;
0416
0417 SqrtsHat = SqrtsHat_;
0418 sHat = SqrtsHat * SqrtsHat;
0419 sHat2 = sHat * sHat;
0420 InvSqrtsHat = 1.0 / SqrtsHat;
0421 InvsHat = InvSqrtsHat * InvSqrtsHat;
0422 InvsHat2 = InvsHat * InvsHat;
0423 t1 = t1_;
0424 t2 = t2_;
0425 Phi1 = Phi1_;
0426 Phi2 = Phi2_;
0427 PPhi = Phi1 - Phi2;
0428 y = y_;
0429
0430
0431 Mju = 0.618 * SqrtsHat;
0432 Mju2 = Mju * Mju;
0433 LnMju2 = log(Mju2);
0434
0435 Pt1 = sqrt(-t1);
0436 Pt2 = sqrt(-t2);
0437
0438 Pt1DotPt2 = Pt1 * Pt2 * cos(PPhi);
0439
0440 x1x2 = sqrt((sHat - t1 - t2 + 2.0 * Pt1DotPt2) * Invs);
0441 Invsx1x2 = Invs / x1x2;
0442
0443 ey = exp(y_);
0444
0445 x1 = x1x2 * ey;
0446 x2 = x1x2 / ey;
0447
0448 CLHEP::HepLorentzVector Glu1, Glu2;
0449
0450 Glu1.setE(x1 * P1In.e());
0451 Glu2.setE(x2 * P2In.e());
0452 Glu1.setPz(x1 * P1In.pz());
0453 Glu2.setPz(x2 * P2In.pz());
0454 Glu1.setPx(-Pt1 * cos(Phi1));
0455 Glu2.setPx(-Pt2 * cos(Phi2));
0456 Glu1.setPy(-Pt1 * sin(Phi1));
0457 Glu2.setPy(-Pt2 * sin(Phi2));
0458
0459 Proton1 = P1In - Glu1;
0460 Proton2 = P2In - Glu2;
0461
0462 CentralVector = Glu1 + Glu2;
0463
0464 double V12 = 1.0 - 4.0 * t1 * Invs / (x1 * x1);
0465 double V22 = 1.0 - 4.0 * t2 * Invs / (x2 * x2);
0466
0467 double V1MinusV2 = sqrt(V12 + V22 + 2.0 - 8.0 * Pt1DotPt2 * Invsx1x2);
0468
0469 InvV1MinusV2 = 1.0 / V1MinusV2;
0470
0471 return;
0472 }
0473
0474 void Exhume::CrossSection::Hadronise() {
0475
0476 int njoin = Partons.size();
0477
0478
0479 std::unique_ptr<int[]> ijoin(new int[njoin]);
0480
0481 double e_, theta_, phi_;
0482 int id, nn, nnc;
0483
0484 double ps_scale;
0485
0486
0487
0488 e_ = Proton1.e();
0489 theta_ = Proton1.theta();
0490 phi_ = Proton1.phi();
0491 nn = 1;
0492 py1ent(nn, Proton1Id, e_, theta_, phi_);
0493
0494 e_ = Proton2.e();
0495 theta_ = Proton2.theta();
0496 phi_ = Proton2.phi();
0497 nn = 2;
0498 py1ent(nn, Proton2Id, e_, theta_, phi_);
0499
0500 double px_, py_, pz_, mm_;
0501
0502 for (int i = 0; i < njoin; i++) {
0503 e_ = Partons[i].p.e();
0504 px_ = Partons[i].p.px();
0505 py_ = Partons[i].p.py();
0506 pz_ = Partons[i].p.pz();
0507 mm_ = Partons[i].p.m();
0508
0509 nn = i + 3;
0510
0511 nnc = nn - 1;
0512 id = Partons[i].id;
0513
0514 for (int j = 0; j < 5; j++) {
0515 pyjets.k[j][nnc] = 0;
0516 pyjets.p[j][nnc] = 0.0;
0517 pyjets.v[j][nnc] = 0.0;
0518 }
0519
0520 pyjets.k[0][nnc] = 1;
0521 pyjets.k[1][nnc] = id;
0522
0523 pyjets.p[0][nnc] = px_;
0524 pyjets.p[1][nnc] = py_;
0525 pyjets.p[2][nnc] = pz_;
0526 pyjets.p[3][nnc] = e_;
0527 pyjets.p[4][nnc] = mm_;
0528
0529 pyjets.n = nn + 1;
0530
0531 ijoin[i] = nn;
0532 }
0533
0534 if ((njoin > 1) && (Name != "di-photon")) {
0535 pyjoin(njoin, ijoin.get());
0536 int ip1 = 3;
0537 int ip2 = 4;
0538 CLHEP::HepLorentzVector b_COM = Partons[0].p.boost(CentralVector.findBoostToCM());
0539
0540
0541 ps_scale = 2 * b_COM.et();
0542 pyshow(ip1, ip2, ps_scale);
0543 }
0544
0545 pyexec();
0546
0547 call_pyhepc(1);
0548
0549 return;
0550 }
0551
0552 void Exhume::CrossSection::AlphaSInit() {
0553 ASConst = 12.0 * M_PI / 22.0 / 2.0;
0554 LambdaQCD = 0.001 * LambdaQCD;
0555 ASFreeze = ASConst / (log(Freeze / LambdaQCD));
0556
0557 return;
0558 }
0559
0560 double Exhume::CrossSection::Fg1Fg2(const double &Qt2_, const double &x1_, const double &x2_) {
0561 double Kt = sqrt(Qt2_);
0562 double grad1, grad2;
0563 double KtLow, KtHigh, upv1, dnv1, upv2, dnv2, usea, dsea, str, chm, bot, top, glLow, glHigh, gl1, gl2;
0564 double AlphaS_ = AlphaS(Kt);
0565
0566 KtLow = 0.9 * Kt;
0567 KtHigh = 1.1 * Kt;
0568 double x_ = x1_;
0569
0570
0571 structm(&x_, &KtLow, &upv1, &dnv1, &usea, &dsea, &str, &chm, &bot, &top, &glLow);
0572 structm(&x_, &KtHigh, &upv1, &dnv1, &usea, &dsea, &str, &chm, &bot, &top, &glHigh);
0573
0574 grad1 = 2.5 * (glHigh - glLow);
0575
0576 structm(&x_, &Kt, &upv1, &dnv1, &usea, &dsea, &str, &chm, &bot, &top, &gl1);
0577
0578 x_ = x2_;
0579
0580 structm(&x_, &KtLow, &upv2, &dnv2, &usea, &dsea, &str, &chm, &bot, &top, &glLow);
0581 structm(&x_, &KtHigh, &upv2, &dnv2, &usea, &dsea, &str, &chm, &bot, &top, &glHigh);
0582
0583 grad2 = 2.5 * (glHigh - glLow);
0584
0585 structm(&x_, &Kt, &upv2, &dnv2, &usea, &dsea, &str, &chm, &bot, &top, &gl2);
0586
0587
0588 double SqrtT = TFast(Qt2_);
0589 double dT = Nc_2PI * AlphaS_ * log(Inv3 * Kt / (Kt + Mju));
0590
0591
0592 double Rg1Rg2_ = Rg1Rg2(Kt);
0593
0594
0595
0596
0597
0598
0599
0600 return (Rg1Rg2_ * (SqrtT * (grad1 - dT * gl1)) *
0601 (SqrtT * (grad2 - dT * gl2)));
0602 }
0603
0604 double Exhume::CrossSection::Fg1Fg2(const int &ii_, const double &x1_, const double &x2_) {
0605 double grad1, grad2;
0606 double upv1, dnv1, upv2, dnv2, usea, dsea, str, chm, bot, top, glLow, glHigh, gl1, gl2;
0607
0608 double x_ = x1_;
0609
0610
0611 structm(&x_, &_KtLow[ii_], &upv1, &dnv1, &usea, &dsea, &str, &chm, &bot, &top, &glLow);
0612 structm(&x_, &_KtHigh[ii_], &upv1, &dnv1, &usea, &dsea, &str, &chm, &bot, &top, &glHigh);
0613
0614 grad1 = 2.5 * (glHigh - glLow);
0615
0616 structm(&x_, &_Qt[ii_], &upv1, &dnv1, &usea, &dsea, &str, &chm, &bot, &top, &gl1);
0617
0618 x_ = x2_;
0619
0620 structm(&x_, &_KtLow[ii_], &upv2, &dnv2, &usea, &dsea, &str, &chm, &bot, &top, &glLow);
0621 structm(&x_, &_KtHigh[ii_], &upv2, &dnv2, &usea, &dsea, &str, &chm, &bot, &top, &glHigh);
0622
0623 grad2 = 2.5 * (glHigh - glLow);
0624
0625 structm(&x_, &_Qt[ii_], &upv2, &dnv2, &usea, &dsea, &str, &chm, &bot, &top, &gl2);
0626
0627
0628 double SqrtT = TFast(_Qt2[ii_]);
0629 double dT = _NcAs_2PI[ii_] * log(Inv3 * _Qt[ii_] / (_Qt[ii_] + Mju));
0630
0631 double Rg1Rg2_ = Rg1Rg2(_Qt[ii_]);
0632
0633
0634
0635
0636
0637
0638
0639 return (Rg1Rg2_ * (SqrtT * (grad1 - dT * gl1)) *
0640 (SqrtT * (grad2 - dT * gl2)));
0641 }
0642
0643 void Exhume::CrossSection::LumiInit() {
0644 MidQt2 = 4.0;
0645 InvMidQt2 = 1.0 / MidQt2;
0646 InvMidQt4 = InvMidQt2 * InvMidQt2;
0647 InvMaxQt2 = 0.0;
0648 LumNStart = 10;
0649 LumAccuracy = 100;
0650 LumNSimps = 51;
0651 LumNSimps_1 = LumNSimps - 1;
0652 TConst = -0.25 / M_PI;
0653
0654 Inv2PI = 0.5 / M_PI;
0655 Nc_2PI = 3.0 * Inv2PI;
0656
0657 Inv3 = 1.0 / 3.0;
0658
0659 Cf_2PIRg = 4.0 * Inv3 * Inv2PI / Rg;
0660
0661 LumSimpsIncr = (InvMidQt2 - InvMaxQt2) / ((double)LumNSimps);
0662 LumSimpsFunc = (double *)calloc(LumNSimps, sizeof(double));
0663 _Qt2 = (double *)calloc(LumNSimps, sizeof(double));
0664 _Qt = (double *)calloc(LumNSimps, sizeof(double));
0665 _KtLow = (double *)calloc(LumNSimps, sizeof(double));
0666 _KtHigh = (double *)calloc(LumNSimps, sizeof(double));
0667 _AlphaS = (double *)calloc(LumNSimps, sizeof(double));
0668 _CfAs_2PIRg = (double *)calloc(LumNSimps, sizeof(double));
0669 _NcAs_2PI = (double *)calloc(LumNSimps, sizeof(double));
0670
0671 if (LumSimpsFunc == nullptr || _Qt == nullptr || _Qt2 == nullptr || _KtLow == nullptr || _KtHigh == nullptr ||
0672 _AlphaS == nullptr || _CfAs_2PIRg == nullptr || _NcAs_2PI == nullptr) {
0673 NoMem();
0674 }
0675
0676
0677
0678 LumSimpsFunc[0] = 0.0;
0679 double InvQt2 = InvMaxQt2 + LumSimpsIncr;
0680 for (int i = 1; i < LumNSimps; i++) {
0681 _Qt2[i] = 1.0 / InvQt2;
0682 _Qt[i] = sqrt(_Qt2[i]);
0683 _KtLow[i] = 0.9 * _Qt[i];
0684 _KtHigh[i] = 1.1 * _Qt[i];
0685 _AlphaS[i] = AlphaS(_Qt[i]);
0686 _CfAs_2PIRg[i] = Cf_2PIRg * _AlphaS[i];
0687 _NcAs_2PI[i] = Nc_2PI * _AlphaS[i];
0688 InvQt2 = InvQt2 + LumSimpsIncr;
0689 }
0690
0691
0692 LumConst = 0.015625 * Survive * M_PI * M_PI;
0693
0694 Tn = 81;
0695 Tn_1 = Tn - 1;
0696 TFunc = (double *)calloc(Tn, sizeof(double));
0697
0698 if (TFunc == nullptr) {
0699 NoMem();
0700 }
0701
0702 P1In.setE(0.5 * root_s);
0703 P2In.setE(0.5 * root_s);
0704 P1In.setPz(0.5 * root_s);
0705 P2In.setPz(-0.5 * root_s);
0706 P1In.setPx(0.0);
0707 P2In.setPx(0.0);
0708 P1In.setPy(0.0);
0709 P2In.setPy(0.0);
0710
0711 std::map<double, double> TempMap;
0712 double QtMin = sqrt(MinQt2);
0713 TempMap.insert(std::pair<double, double>(QtMin, 1.0));
0714 TempMap.insert(std::pair<double, double>(100.0, 1.0));
0715
0716 RgMap2d.insert(std::pair<double, std::map<double, double> >(0.0, TempMap));
0717
0718 TempMap.clear();
0719 TempMap.insert(std::pair<double, double>(QtMin, Rg_(0.999, QtMin)));
0720 TempMap.insert(std::pair<double, double>(100.0, Rg_(0.999, 100.0)));
0721
0722 RgMap2d.insert(std::pair<double, std::map<double, double> >(0.999, TempMap));
0723
0724 TempMap.clear();
0725
0726 Mju = 10.0;
0727 Mju2 = Mju * Mju;
0728 LnMju2 = log(Mju2);
0729
0730 TempMap.insert(std::pair<double, double>(MinQt2, T(MinQt2)));
0731 TempMap.insert(std::pair<double, double>(10000.0, T(10000.0)));
0732
0733 TMap2d.insert(std::pair<double, std::map<double, double> >(Mju, TempMap));
0734
0735 TempMap.clear();
0736
0737 Mju = 1000.0;
0738 Mju2 = Mju * Mju;
0739 LnMju2 = log(Mju2);
0740
0741 TempMap.insert(std::pair<double, double>(MinQt2, T(MinQt2)));
0742 TempMap.insert(std::pair<double, double>(10000.0, T(10000.0)));
0743
0744 TMap2d.insert(std::pair<double, std::map<double, double> >(Mju, TempMap));
0745
0746 TempMap.clear();
0747
0748 return;
0749 }
0750
0751 double Exhume::CrossSection::Lumi_() {
0752 double lum;
0753 fg2Map.clear();
0754 std::map<double, double> NextM, M;
0755 std::map<double, double>::iterator high, low, q1, q2;
0756 double LumIncr;
0757
0758 double fg2_Qt4;
0759 double Qt2;
0760
0761
0762
0763
0764 for (Qt2 = MinQt2; Qt2 < 0.0000004; Qt2 += 0.0000001) {
0765 fg2_Qt4 = Fg1Fg2(Qt2, x1, x2) / (Qt2 * Qt2);
0766
0767 if (fg2_Qt4 != fg2_Qt4)
0768 fg2_Qt4 = 0.0;
0769
0770 fg2Map.insert(fEntry(Qt2, fg2_Qt4));
0771 }
0772
0773 fg2_Qt4 = Fg1Fg2(MidQt2, x1, x2) * InvMidQt4;
0774
0775 fg2Map.insert(fEntry(MidQt2, fg2_Qt4));
0776
0777 lum = 1.0;
0778
0779 NextM.insert(fEntry(0.0, MinQt2));
0780 NextM.insert(fEntry(lum, MidQt2));
0781
0782 unsigned int N = LumNStart;
0783 unsigned int NTot = LumNStart;
0784
0785 while (NTot < LumAccuracy) {
0786 N = N * 2;
0787 NTot = NTot + N;
0788 M = NextM;
0789 NextM.clear();
0790 NextM.insert(fEntry(0.0, MinQt2));
0791 LumIncr = lum / (N - 1.0);
0792
0793 for (double q_ = LumIncr; q_ < lum; q_ += LumIncr) {
0794 high = M.upper_bound(q_);
0795 if (high == M.end())
0796 high--;
0797 low = high;
0798 low--;
0799
0800 Qt2 = low->second + (high->second - low->second) * (q_ - low->first) / (high->first - low->first);
0801
0802 fg2_Qt4 = Fg1Fg2(Qt2, x1, x2) / (Qt2 * Qt2);
0803 fg2Map.insert(fEntry(Qt2, fg2_Qt4));
0804 }
0805
0806 q1 = fg2Map.begin();
0807 q2 = q1;
0808
0809 for (q2++; q2 != fg2Map.end(); q2++) {
0810 NextM.insert(
0811 fEntry(NextM.rbegin()->first + 0.5 * (q2->first - q1->first) * (q2->second + q1->second), q2->first));
0812 q1 = q2;
0813 }
0814 lum = NextM.rbegin()->first;
0815 }
0816
0817
0818
0819
0820
0821
0822
0823 for (int i = 1; i < LumNSimps; i++) {
0824
0825 LumSimpsFunc[i] = Fg1Fg2(i, x1, x2);
0826 }
0827 lum = lum + dsimps_(LumSimpsFunc, &InvMaxQt2, &InvMidQt2, &LumNSimps_1);
0828
0829
0830
0831 return (LumConst * lum * lum);
0832 }
0833
0834 double Exhume::CrossSection::Lumi() {
0835
0836
0837
0838 return (Lumi_());
0839 }
0840
0841 void Exhume::CrossSection::NoMem() {
0842 std::cout << " Your computer cannot allocate enough memory." << std::endl
0843 << " Try replacing it with one that wouldn't be more useful as a doorstop." << std::endl
0844 << " ...Exiting" << std::endl;
0845
0846 exit(0);
0847
0848 return;
0849 }
0850
0851 double Exhume::CrossSection::Splitting(const double &Kt) {
0852 double D = Mju / (Mju + Kt);
0853 double D2 = D * D;
0854 double D3 = D2 * D;
0855
0856 unsigned int Nf;
0857 if (Kt < CharmMass) {
0858 Nf = 3;
0859 } else if (Kt < BottomMass) {
0860 Nf = 4;
0861 } else {
0862 Nf = 5;
0863 }
0864
0865 return (6.0 * (-D2 + Inv3 * D3 - 0.25 * D3 * D - log(1 - D)) + Nf * (0.5 * (D - D2) + Inv3 * D3));
0866 }
0867
0868 double Exhume::CrossSection::TFast(const double &Qt2_) {
0869 if (TBegin) {
0870 TBegin = false;
0871
0872 TMjuHigh = TMap2d.upper_bound(Mju);
0873 TMjuLow = TMjuHigh;
0874 TMjuLow--;
0875
0876 if (TMjuHigh->first - TMjuLow->first > 0.01 * Mju || TMjuHigh == TMap2d.end() || TMjuHigh == TMap2d.begin()) {
0877 std::map<double, double> TTemp;
0878
0879 TMap2d.insert(std::pair<double, std::map<double, double> >(Mju, TTemp));
0880 TMjuHigh = TMap2d.upper_bound(0.9999 * Mju);
0881 TInterpolate = false;
0882
0883 }
0884 }
0885
0886 if (!TInterpolate) {
0887
0888 double Tee = T(Qt2_);
0889 (TMjuHigh->second).insert(std::pair<double, double>(Qt2_, Tee));
0890 return (Tee);
0891 }
0892
0893 if (TInterpolate) {
0894
0895
0896 double THigh, TLow;
0897 std::map<double, double>::iterator high_, low_;
0898 high_ = (TMjuHigh->second).upper_bound(Qt2_);
0899 bool CalculateT = true;
0900
0901 if (high_ != (TMjuHigh->second).end() && high_ != (TMjuHigh->second).begin()) {
0902 low_ = high_;
0903 low_--;
0904 if (high_->first - low_->first < 0.01 * Qt2_)
0905 CalculateT = false;
0906 }
0907
0908 if (CalculateT) {
0909 double MjuOld = Mju;
0910 double Mju2Old = Mju2;
0911 double LnMju2Old = LnMju2;
0912
0913
0914 Mju = TMjuHigh->first;
0915 Mju2 = Mju * Mju;
0916 LnMju2 = log(Mju2);
0917 THigh = T(Qt2_);
0918 Mju = MjuOld;
0919 Mju2 = Mju2Old;
0920 LnMju2 = LnMju2Old;
0921 (TMjuHigh->second).insert(std::pair<double, double>(Qt2_, THigh));
0922 } else {
0923 THigh = low_->second + (Qt2_ - low_->first) * (high_->second - low_->second) / (high_->first - low_->first);
0924 }
0925
0926 CalculateT = true;
0927 high_ = (TMjuLow->second).upper_bound(Qt2_);
0928
0929 if (high_ != (TMjuLow->second).end() && high_ != (TMjuLow->second).begin()) {
0930 low_ = high_;
0931 low_--;
0932 if (high_->first - low_->first < 0.01 * Qt2_)
0933 CalculateT = false;
0934 }
0935
0936 if (CalculateT) {
0937 double MjuOld = Mju;
0938 double Mju2Old = Mju2;
0939 double LnMju2Old = LnMju2;
0940
0941
0942 Mju = TMjuLow->first;
0943 Mju2 = Mju * Mju;
0944 LnMju2 = log(Mju2);
0945 TLow = T(Qt2_);
0946 Mju = MjuOld;
0947 Mju2 = Mju2Old;
0948 LnMju2 = LnMju2Old;
0949 (TMjuLow->second).insert(std::pair<double, double>(Qt2_, TLow));
0950 } else {
0951 TLow = low_->second + (Qt2_ - low_->first) * (high_->second - low_->second) / (high_->first - low_->first);
0952 }
0953
0954 return (TLow + (Mju - TMjuLow->first) * (THigh - TLow) / (TMjuHigh->first - TMjuLow->first));
0955 }
0956
0957 return (0.0);
0958 }
0959
0960 double Exhume::CrossSection::T(const double &Qt2_) {
0961 double Kt;
0962 double LnQt2 = log(Qt2_);
0963 double TIncr = (LnMju2 - LnQt2) / (Tn - 1.0);
0964 double TLnKt2 = LnQt2;
0965
0966 for (int i = 0; i < Tn; i++) {
0967 Kt = sqrt(exp(TLnKt2));
0968 TFunc[i] = Splitting(Kt) * AlphaS(Kt);
0969 TLnKt2 = TLnKt2 + TIncr;
0970 }
0971
0972
0973 if (Qt2_ > Mju2) {
0974 return (exp(-TConst * dsimps_(TFunc, &LnQt2, &LnMju2, &Tn_1)));
0975 }
0976 return (exp(TConst * dsimps_(TFunc, &LnQt2, &LnMju2, &Tn_1)));
0977 }
0978
0979 double Exhume::CrossSection::Txg(const double &Qt2_, const double &x_) {
0980 double x__ = x_;
0981 double Qt = sqrt(Qt2_);
0982 double upv, dnv, usea, dsea, str, chm, bot, top, gl;
0983 structm(&x__, &Qt, &upv, &dnv, &usea, &dsea, &str, &chm, &bot, &top, &gl);
0984
0985
0986 return (T(Qt2_) * gl);
0987 }
0988
0989
0990
0991
0992
0993
0994
0995
0996
0997
0998
0999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039 std::complex<double> Exhume::CrossSection::Fsf(const double &Tau_) {
1040 double InvTau = 1.0 / Tau_;
1041 return (InvTau * (1.0 + (1.0 - InvTau) * f(Tau_)));
1042 }
1043
1044 std::complex<double> Exhume::CrossSection::F0(const double &Tau_) {
1045 double InvTau = 1.0 / Tau_;
1046 return (InvTau * (-1.0 + InvTau * f(Tau_)));
1047 }
1048
1049 std::complex<double> Exhume::CrossSection::f(const double &Tau_) {
1050 std::complex<double> Sqrtf;
1051 std::complex<double> f;
1052
1053 if (Tau_ <= 1.0) {
1054 Sqrtf = asin(sqrt(Tau_));
1055 f = Sqrtf * Sqrtf;
1056 }
1057 if (Tau_ > 1.0) {
1058 std::complex<double> SqrtTau = sqrt(Tau_), SqrtTau_1 = sqrt(Tau_ - 1.0);
1059 Sqrtf = log((SqrtTau + SqrtTau_1) / (SqrtTau - SqrtTau_1)) - I * M_PI;
1060 f = -0.25 * Sqrtf * Sqrtf;
1061 }
1062 return (f);
1063 }
1064