ExhumeInterface/src/Event.cpp

0001 //-*-c++-*-
0002 //-*-Event.cpp-*-
0003 //   Written by James Monk and Andrew Pilkington
0004 /////////////////////////////////////////////////////////////////////////////
0005
0006 #include "GeneratorInterface/ExhumeInterface/interface/Event.h"
0007 #include "CLHEP/Random/RandomEngine.h"
0008 //////////////////////////////////////////////////////////////////////////////
0009 Exhume::Event::Event(CrossSection &Process_, CLHEP::HepRandomEngine *engine) {
0010   Process = &Process_;
0011   Setx1Max(1.0);
0012   Setx2Max(1.0);
0013   Sett1Max(0.0);
0014   Sett2Max(0.0);
0015   Sett1Min(-10.0);
0016   Sett2Min(-10.0);
0017   SetMassRange(100, 130);
0018   NumberOfEvents = 0;
0019   TotalAttempts = 0;
0020   ymax = 0.0;
0021   ymin = 0.0;
0022   CSMass = 0.0;
0023   Sigmai = 0.0;
0024   yRange = 0.0;
0025   TwoPI = 2.0 * M_PI;
0026   B = Process->GetB();
0027   InvB = 1.0 / B;
0028   InvBlnB = InvB * log(B);
0029   Root_s = Process->GetRoot_s();
0030   InvRoot_s = 1.0 / Root_s;
0031   randomEngine = engine;
0032   std::cout << "   ++ Random number seed   " << randomEngine->getSeed() << " ++" << std::endl << std::endl;
0033   Process->SetRandomEngine(engine);
0034   //srand(Seed_);
0035 }
0036 ///////////////////////////////////////////////////////////////////////////////
0037 Exhume::Event::~Event() {}
0038 ///////////////////////////////////////////////////////////////////////////////
0039 void Exhume::Event::Generate() {
0040   bool pass = false;
0041   do {
0042     SelectValues();
0043     Process->SetKinematics(SqrtsHat, Eta, t1, t2, Phi1, Phi2);
0044     Process->SetSubParameters();
0045     Process->SetPartons();
0046     CSi = Process->Differential();
0047
0048     wgt = GetFunc(SqrtsHat);
0049     CSMass += wgt;
0050     double ebt = exp(B * (t1 + t2 - t1Max - t2Max)) * Process->SubParameterWeight();
0051     wgt = 1.1 * wgt * ebt;
0052
0053     if (CSi > wgt) {  //Don't want this to happen very often - it is slow.
0054                       //Choose your cuts sensibly!
0055       Process->MaximiseSubParameters();
0056       wgt = WeightFunc(SqrtsHat);
0057       AddPoint(SqrtsHat, wgt);
0058       wgt = 1.1 * wgt * ebt;
0059       Process->SetKinematics(SqrtsHat, Eta, t1, t2, Phi1, Phi2);
0060       Process->SetSubParameters();
0061     }
0062
0063     CSi = CSi / wgt;
0064
0065     Sigmai += CSi;
0066
0067     if (CSi > VonNeu) {
0068       if (CSi > 1.0) {  //Then it has cocked it up completely
0069
0070         std::map<double, double> FuncMap_ = GetFuncMap();
0071         std::map<double, double> LineShape_ = GetLineShape();
0072
0073         FILE *dat = fopen("FuncMap.dat", "w");
0074
0075         for (std::map<double, double>::iterator kk = FuncMap_.begin(); kk != FuncMap_.end(); kk++) {
0076           //fprintf(dat,"%le\t%le\n",kk->first, kk->second);
0077           fprintf(dat, "%e\t%e\n", kk->first, kk->second);
0078         }
0079
0080         fclose(dat);
0081
0082         dat = fopen("LineShape.dat", "w");
0083
0084         for (std::map<double, double>::iterator kk = LineShape_.begin(); kk != LineShape_.end(); kk++) {
0085           //fprintf(dat, "%le\t%le\n",kk->first, kk->second);
0086           fprintf(dat, "%e\t%e\n", kk->first, kk->second);
0087         }
0088         fclose(dat);
0089
0090         std::cout << "   This should never happen" << std::endl;
0091         std::cout << "   m   = " << SqrtsHat << std::endl;
0092         std::cout << "   y   = " << Eta << std::endl;
0093         std::cout << "   t1  = " << t1 << std::endl;
0094         std::cout << "   t2  = " << t2 << std::endl;
0095         std::cout << "   x1  = " << Process->Getx1() << std::endl;
0096         std::cout << "   x2  = " << Process->Getx2() << std::endl;
0097         std::cout << "   CSi = " << CSi * wgt << std::endl;
0098         std::cout << "   wgt = " << wgt << std::endl;
0099         std::cout << "   Please e-mail the authors." << std::endl
0100                   << "   Include the files FuncMap.dat and LineShape.dat" << std::endl;
0101         exit(0);
0102       }
0103       pass = true;
0104       NumberOfEvents++;
0105       TotalAttempts++;
0106
0107       Var.push_back(std::pair<double, double>(TotalAttempts * TotalIntegral, Sigmai));
0108
0109     } else {
0110       TotalAttempts++;
0111     }
0112   } while (!pass);
0113 }
0114 //////////////////////////////////////////////////////////////////////////////
0115 void Exhume::Event::SelectValues() {
0116   double MRand = randomEngine->flat();
0117
0118   SqrtsHat = GetValue(MRand);
0119
0120   double tt1 = randomEngine->flat() * (tt1max - tt1min) + tt1min;
0121   double tt2 = randomEngine->flat() * (tt2max - tt2min) + tt2min;
0122
0123   t1 = InvB * log(tt1) + InvBlnB;
0124   t2 = InvB * log(tt2) + InvBlnB;
0125
0126   Phi1 = TwoPI * randomEngine->flat();
0127   Phi2 = TwoPI * randomEngine->flat();
0128
0129   // y = 0.5 log(x1/x2) and at t1 = t2 = 0 x1x2 = sHat/s
0130   // so the max kinematically allowed value of x1 is
0131   // x1max = sHat/ s / x2min
0132   // x2max =  sHat/ s / x1min and so
0133   // ymax = 0.5 log(x1max / x2min) = 0.5 log(x1max^2 * s/sHat)
0134   // ymin = 0.5 log(x1min / x2max) = -0.5 log (x2max^2 * s/sHat)
0135
0136   double yymax = log(x1Max * Root_s / SqrtsHat);
0137   double yymin = log(InvRoot_s * SqrtsHat / x2Max);
0138
0139   yRange += (yymax - yymin);
0140
0141   Eta = (yymax - yymin) * randomEngine->flat() + yymin;
0142
0143   VonNeu = randomEngine->flat();
0144
0145   return;
0146 }
0147 //////////////////////////////////////////////////////////////////////////////
0148 void Exhume::Event::SetParameterSpace() {
0149   tt1max = InvB * exp(t1Max * B);
0150   tt2max = InvB * exp(t2Max * B);
0151   tt1min = InvB * exp(t1Min * B);
0152   tt2min = InvB * exp(t2Min * B);
0153
0154   Process->MaximiseSubParameters();
0155
0156   WeightInit(MinMass, MaxMass);
0157
0158   return;
0159 }
0160 //////////////////////////////////////////////////////////////////////////////
0161 double Exhume::Event::WeightFunc(const double &mm_) {
0162   Process->SetKinematics(mm_, 0.0, t1Max, t2Max, 0.0, 0.0);
0163
0164   if (mm_ < 65.0) {
0165     return (1.00025 * Process->Differential());  //Factor 1.00025 accounts for
0166                                                  //slight kink in lumi function
0167   }
0168
0169   return (Process->Differential());
0170 }
0171 //////////////////////////////////////////////////////////////////////////////
0172 double Exhume::Event::CrossSectionCalculation() {
0173   double InvTotalAttempts = 1.0 / TotalAttempts;
0174   double yy = yRange * InvTotalAttempts;
0175
0176   //\phi dependence is already integrated out of the differential luminosity
0177
0178   return (1.1 * TotalIntegral * (tt1max - tt1min) * (tt2max - tt2min) * yy * fabs(Process->SubParameterRange()) *
0179           Sigmai * InvTotalAttempts);
0180 }
0181 //////////////////////////////////////////////////////////////////////////////