Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 //-*-c++-*-
0002 //-*-CrossSection.cpp-*-
0003 //   Written by James Monk and Andrew Pilkington
0004 /////////////////////////////////////////////////////////////////////////////
0005 
0006 #include "GeneratorInterface/ExhumeInterface/interface/CrossSection.h"
0007 #include "GeneratorInterface/ExhumeInterface/interface/PythiaRecord.h"
0008 
0009 #include "HepMC/PythiaWrapper6_4.h"
0010 //#include "CLHEP/HepMC/ConvertHEPEVT.h"
0011 //#include "CLHEP/HepMC/CBhepevt.h"
0012 #include <cstdio>
0013 #include <memory>
0014 #include <cmath>
0015 
0016 // External Fortran routines to link to:
0017 double dsimps_(double *, double *, double *, int *);
0018 
0019 //void call_pyhepc(int);
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 const int pyjets_maxn =4000;
0046 extern struct {
0047     int n, npad, k[5][pyjets_maxn];
0048     double p[5][pyjets_maxn], v[5][pyjets_maxn];
0049 } pyjets_;
0050 #define pyjets pyjets_
0051 
0052 extern struct {
0053     int kchg[4][500];
0054     double pmas[4][500], parf[2000], vckm[4][4];
0055 } pydat2_;
0056 #define pydat2 pydat2_
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   //Put data types into a map and pair with a string
0095   //for formating in/output.
0096   /*TypeMap.insert(PCharPair(typeid(double*).name(),"%lf"));
0097   TypeMap.insert(PCharPair(typeid(float*).name(),"%f"));
0098   TypeMap.insert(PCharPair(typeid(int*).name(),"%d"));*/
0099 
0100   //Associate each variable with a string
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   //(Re-)Compute rest of parameters
0128   Freeze = sqrt(MinQt2);
0129 
0130   /*if(FNAL_or_LHC == 0){
0131     Survive = 0.045;
0132     Rg = 1.4;
0133     root_s = 1960;
0134   }else if(FNAL_or_LHC==1){
0135     Survive = 0.03;
0136     Rg = 1.2;
0137     root_s = 14000;
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   //Now write out all the parameters:
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   //Initialise the PDFs
0163   // setting up lhapdf path name from environment varaible (***)
0164   /*char* lhaPdfs = NULL;
0165   std::cout << std::endl;   
0166   std::cout<<" Trying to find LHAPATH in environment ...";
0167   lhaPdfs = getenv("LHAPATH");
0168   if(lhaPdfs != NULL) {
0169     std::cout<<" done."<<std::endl;
0170     lhapdfSetPath_=std::string(lhaPdfs);
0171     std::cout<<" Using "<< lhapdfSetPath_ << std::endl; 
0172   }
0173   else{
0174     std::cout<<" failed."<<std::endl;
0175     std::cout<<" Using "<< lhapdfSetPath_ << std::endl;
0176   }*/
0177 
0178   std::cout << std::endl
0179             << " ........................................................................." << std::endl
0180             << std::endl
0181             << "  = Initialising PDFs =" << std::endl;
0182   /*char pdfpath[232];
0183   bool dot=false;
0184   for(int i=0; i<232; ++i) {
0185     if(lhapdfSetPath_.c_str()[i]=='\0') dot=true;
0186     if(!dot) pdfpath[i]=lhapdfSetPath_.c_str()[i];
0187     else pdfpath[i]=' ';
0188   }
0189 
0190   setpdfpath(pdfpath);*/
0191 
0192   my_pdfset(PDF);
0193   std::cout << std::endl << " ........................................................................." << std::endl;
0194   //...........................................................................
0195 
0196   //double smass=0.0;//for now
0197   my_pythia_init();
0198   //pydata();
0199   //initpydata();
0200   pyinre();
0201   Proton1Id = 2212;
0202   //Proton2Id = (FNAL_or_LHC==0)?-2212:2212;
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     //std::cout<<"   Beginning Rg calculation"<<std::endl;
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     for(int i=0;i<2;i++){
0276       if(RgInterpolate[i])std::cout<<"   Interpolating Rg"<<std::endl;
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   //if(dlam>1.0) return(1.0);
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   //Mju = 0.5 * SqrtsHat;
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);  //This is the 3 vector between the protons
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   //int one = 1;
0476   int njoin = Partons.size();
0477 
0478   //NOTE: all values initialized in for loop below
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   //std::cout << "--------->DEBUG: MSTU(4) = "<< pydat1.mstu[3] << std::endl;
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;    //Fortran indexing starts at 1 and have already
0510                    //Inserted the 2 protons.
0511     nnc = nn - 1;  //C++ array indexing;
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;   //states there is an undecayed particle/parton
0521     pyjets.k[1][nnc] = id;  //sets the PDG type
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     //HepLorentzVector bb_COM = Partons[1].p.boost(CentralVector.findBoostToCM());
0540     //std::cout<<b_COM.perp()<<"\t"<<bb_COM.perp()<<std::endl;
0541     ps_scale = 2 * b_COM.et();  //Q^2 = 4pt^2 for final state parton showering
0542     pyshow(ip1, ip2, ps_scale);
0543   }
0544 
0545   pyexec();
0546   //pyhepc_(one);
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   //structm returns gl*x:
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   //double SqrtT = T(Qt2_);
0588   double SqrtT = TFast(Qt2_);
0589   double dT = Nc_2PI * AlphaS_ * log(Inv3 * Kt / (Kt + Mju));
0590   //double CfAs_2PIRg = 0.0;//Cf_2PIRg * AlphaS_;
0591 
0592   double Rg1Rg2_ = Rg1Rg2(Kt);
0593   //double Rg1Rg2_ = Rg_(x1, Kt) * Rg_(x2, Kt);
0594   //double Rg1Rg2_ = 1.2*1.2;
0595 
0596   //dT = 0.0;
0597 
0598   //grad1 = 0.0;
0599 
0600   return (Rg1Rg2_ * (SqrtT * (grad1 - dT * gl1)) *  // - CfAs_2PIRg * (upv1 + dnv1) ) *
0601           (SqrtT * (grad2 - dT * gl2)));            //- CfAs_2PIRg * (upv2 + dnv2) ) );
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   //structm returns gl*x:
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   //double SqrtT = T(_Qt2[ii_]);
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   //double Rg1Rg2_ = Rg_(x1, _Qt[ii_]) * Rg_(x2, _Qt[ii_]);
0633   //double Rg1Rg2_ = 1.2*1.2;
0634 
0635   //dT = 0.0;
0636 
0637   //  grad1 = 0.0;
0638 
0639   return (Rg1Rg2_ * (SqrtT * (grad1 - dT * gl1)) *  // - _CfAs_2PIRg[ii_] * (upv1 + dnv1) ) *
0640           (SqrtT * (grad2 - dT * gl2)));            //- _CfAs_2PIRg[ii_] * (upv2 + dnv2)));
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   //Can't calculate Qt at InvQt = 0.0
0677   //but function converges to 0 there:
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   //LumConst = 0.015625 * Survive * Rg * Rg * Rg * Rg * PI * PI;
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   //Do the integral between MinQt and MidQt:
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   //Now do the rest of the integral between MidQt and MaxQt = 1.0/InvMaxQt
0820   //Transform the integration variable to 1/Qt so \inf -> 0
0821   //Have already entered the first element - it is always 0.
0822 
0823   for (int i = 1; i < LumNSimps; i++) {
0824     //LumSimpsFunc[i] = Fg1Fg2(_Qt2[i], x1, x2);
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   //Will modify this to create lookup table for Lumi on the fly - later
0836   //For now just return Lumi_
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       //std::cout<<"Got here"<<std::endl;
0883     }
0884   }
0885 
0886   if (!TInterpolate) {
0887     //std::cout<<"Not interpolating T"<<std::endl;
0888     double Tee = T(Qt2_);
0889     (TMjuHigh->second).insert(std::pair<double, double>(Qt2_, Tee));
0890     return (Tee);
0891   }
0892 
0893   if (TInterpolate) {
0894     //std::cout<<"Interpolating T"<<std::endl;
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       //This is a bit dirty, don't let anyone else screw with this
0913       //'cause I know what I'm doing:
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       //This is a bit dirty, don't let anyone else screw with this
0941       //'cause I know what I'm doing:
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   //N.B. This is really the square root of T.
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   //gl is really gl * x so don't need to multiply by x
0986   return (T(Qt2_) * gl);
0987 }
0988 /////////////////////////////////////////////////////////////////////////////
0989 /*void Exhume::CrossSection::ReadCard(const std::string &filename){
0990       
0991   FILE * card;
0992   card = fopen(filename.c_str(),"r");
0993   if(card==NULL){
0994     std::cout<<std::endl<<"  File "<<filename<<" does not exist"<<
0995       std::endl<<std::endl<<"\t...Exiting"<<std::endl<<std::endl;
0996     exit(0);
0997   }else{
0998     std::cout<<std::endl<<"  Opening file "<<filename<<"..."<<
0999       std::endl<<std::endl; 
1000   }
1001 
1002   char text[50];
1003   int read;
1004   std::pair<const char*,const void*> value;
1005   std::map<std::string,std::pair<const char*,const void*> >
1006     ::iterator val_;
1007   int LinesRead = 0;
1008 
1009 
1010   while(!feof(card)){
1011     read = fscanf(card,"%s",text);
1012     if(read==1){
1013       val_ = PMap.find(text);
1014 
1015       if(val_!=PMap.end()){
1016     value = val_->second;
1017     LinesRead++;
1018     if(LinesRead==1){
1019       std::cout<<"  Reading values for:";
1020       printf("%20s\n",text);
1021     }else{
1022       printf("\t\t\t%17s\n",text);
1023     }
1024     fscanf(card,TypeMap.find(value.first)->second,value.second);
1025       }
1026     }
1027   }
1028 
1029   if(LinesRead==0) std::cout<<"  WARNING: No correct parameters found in "<<
1030     filename<<std::endl<<std::endl<<
1031     "  No values read in   --   using defaults"<<std::endl;
1032   
1033   std::cout<<std::endl<<"  ...Closing file "<<filename<<std::endl<<std::endl;
1034   
1035   fclose(card);
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 //////////////////////////////////////////////////////////////////////////////
[ 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!