File indexing completed on 2024-04-06 12:13:30
0001
0002
0003
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
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) {
0054
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) {
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
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
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
0130
0131
0132
0133
0134
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());
0166
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
0177
0178 return (1.1 * TotalIntegral * (tt1max - tt1min) * (tt2max - tt2min) * yy * fabs(Process->SubParameterRange()) *
0179 Sigmai * InvTotalAttempts);
0180 }
0181