Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:29:40

0001 //

0002 // CMSCGEN.cc   version 3.0    Thomas Hebbeker 2007-05-15

0003 //

0004 // implemented in CMSSW by P. Biallass 2007-05-28

0005 // see header for documentation and CMS internal note 2007 "Improved Parametrization of the Cosmic Muon Flux for the generator CMSCGEN" by Biallass + Hebbeker

0006 //

0007 
0008 #include <CLHEP/Random/RandomEngine.h>
0009 #include <CLHEP/Random/JamesRandom.h>
0010 
0011 #include "GeneratorInterface/CosmicMuonGenerator/interface/CMSCGEN.h"
0012 
0013 CMSCGEN::CMSCGEN() : initialization(0), RanGen2(nullptr), delRanGen(false) {}
0014 
0015 CMSCGEN::~CMSCGEN() {
0016   if (delRanGen)
0017     delete RanGen2;
0018 }
0019 
0020 void CMSCGEN::setRandomEngine(CLHEP::HepRandomEngine *v) {
0021   if (delRanGen)
0022     delete RanGen2;
0023   RanGen2 = v;
0024   delRanGen = false;
0025 }
0026 
0027 int CMSCGEN::initialize(double pmin_in,
0028                         double pmax_in,
0029                         double thetamin_in,
0030                         double thetamax_in,
0031                         CLHEP::HepRandomEngine *rnd,
0032                         bool TIFOnly_constant,
0033                         bool TIFOnly_linear) {
0034   if (delRanGen)
0035     delete RanGen2;
0036   RanGen2 = rnd;
0037   delRanGen = false;
0038 
0039   //set bools for TIFOnly options (E<2GeV with unphysical energy dependence)

0040   TIFOnly_const = TIFOnly_constant;
0041   TIFOnly_lin = TIFOnly_linear;
0042 
0043   // units: GeV

0044 
0045   // WARNING: coordinate system:

0046   //   - to outside world define z axis downwards, i.e.

0047   //          muon coming from above, vertically: cos = 1

0048   //      (used for backward compatibility)

0049   //   - internally use frame with z axis upwards, i.e.

0050   //          muon coming from above, vertically: cos = -1

0051   //     (corresponds to CMS note definition)

0052 
0053   //set cmin and cmax, here convert between coordinate systems:

0054   cmin_in = -TMath::Cos(thetamin_in);  //input angle already converted from Deg to Rad!

0055   cmax_in = -TMath::Cos(thetamax_in);  //input angle already converted from Deg to Rad!

0056 
0057   //allowed energy range

0058   pmin_min = 3.;
0059   //pmin_max = 100.;

0060   pmin_max = 3000.;
0061   pmax = 3000.;
0062   //allowed angular range

0063   //cmax_max = -0.1,

0064   cmax_max = -0.01, cmax_min = -0.9999;
0065 
0066   if (TIFOnly_const == true || TIFOnly_lin == true)
0067     pmin_min = 0.;  //forTIF

0068 
0069   // set pmin

0070   if (pmin_in < pmin_min || pmin_in > pmin_max) {
0071     std::cout << " >>> CMSCGEN.initialize <<< warning: illegal pmin_in =" << pmin_in;
0072     return (-1);
0073   } else if (pmax_in > pmax) {
0074     std::cout << " >>> CMSCGEN.initialize <<< warning: illegal pmax_in =" << pmax_in;
0075     return (-1);
0076   } else {
0077     pmin = pmin_in;
0078     pmax = pmax_in;
0079     xemax = 1. / (pmin * pmin);
0080     xemin = 1. / (pmax * pmax);
0081   }
0082 
0083   // set cmax and cmin

0084   if (cmax_in < cmax_min || cmax_in > cmax_max) {
0085     std::cout << " >>> CMSCGEN.initialize <<< warning: illegal cmax_in =" << cmax_in;
0086     return (-1);
0087   } else {
0088     cmax = cmax_in;
0089     cmin = cmin_in;
0090   }
0091 
0092   initialization = 1;
0093 
0094   if (TIFOnly_const == true || TIFOnly_lin == true)
0095     pmin_min = 3.;  //forTIF

0096 
0097   //  Lmin = log10(pmin_min);

0098   Lmin = log10(pmin);
0099   Lmax = log10(pmax);
0100   Lfac = 100. / (Lmax - Lmin);
0101 
0102   //

0103   // +++ coefficients for energy spectrum

0104   //

0105 
0106   pe[0] = -1.;
0107   pe[1] = 6.22176;
0108   pe[2] = -13.9404;
0109   pe[3] = 18.1643;
0110   pe[4] = -9.22784;
0111   pe[5] = 1.99234;
0112   pe[6] = -0.156434;
0113   pe[7] = 0.;
0114   pe[8] = 0.;
0115 
0116   //

0117   // +++ coefficients for cos theta distribution

0118   //

0119 
0120   b0c[0] = 0.6639;
0121   b0c[1] = -0.9587;
0122   b0c[2] = 0.2772;
0123 
0124   b1c[0] = 5.820;
0125   b1c[1] = -6.864;
0126   b1c[2] = 1.367;
0127 
0128   b2c[0] = 10.39;
0129   b2c[1] = -8.593;
0130   b2c[2] = 1.547;
0131 
0132   //

0133   // +++ calculate correction table for different cos theta dependence!

0134   //     reference range: c1  to  c2

0135   //

0136   //  explanation: the parametrization of the energy spectrum as used above

0137   //    is the integral over the c = cos(zenith angle) range -1...-0.1

0138   //    since the c distribution depends on energy, the integrated energy

0139   //    spectrum depends on this range. Here a correction factor is determined,

0140   //    based on the linear c dependence of the c distribution.

0141   //    The correction is calculated for 100 bins in L = log10(energy).

0142   //

0143   // +++ in same loop calculate integrated flux

0144   //      (integrated over angles and momentum)

0145 
0146   c1 = -1.;
0147   c2 = -0.1;
0148 
0149   double cemax0 = 1.0;
0150   double L, L2;
0151   double s;
0152   double p, p1, p2;
0153   double integral_here, integral_ref;
0154   double c_cut;
0155 
0156   integrated_flux = 0.;
0157 
0158   for (int k = 1; k <= 100; k++) {
0159     L = Lmin + (k - 0.5) / Lfac;
0160     L2 = L * L;
0161     p = pow(10, L);
0162     p1 = pow(10, L - 0.5 / Lfac);
0163     p2 = pow(10, L + 0.5 / Lfac);
0164 
0165     b0 = b0c[0] + b0c[1] * L + b0c[2] * L2;
0166     b1 = b1c[0] + b1c[1] * L + b1c[2] * L2;
0167     b2 = b2c[0] + b2c[1] * L + b2c[2] * L2;
0168 
0169     // cut out explicitly regions of zero flux

0170     // (for low momentum and near horizontal showers)

0171     // since parametrization for z distribution doesn't work here

0172     // (can become negative)

0173 
0174     c_cut = -0.42 + L * 0.35;
0175 
0176     if (c_cut > c2)
0177       c_cut = c2;
0178 
0179     integral_ref =
0180         b0 * (c_cut - c1) + b1 / 2. * (c_cut * c_cut - c1 * c1) + b2 / 3. * (c_cut * c_cut * c_cut - c1 * c1 * c1);
0181 
0182     if (c_cut > cmax)
0183       c_cut = cmax;
0184 
0185     integral_here = b0 * (c_cut - cmin) + b1 / 2. * (c_cut * c_cut - cmin * cmin) +
0186                     b2 / 3. * (c_cut * c_cut * c_cut - cmin * cmin * cmin);
0187 
0188     corr[k] = integral_here / integral_ref;
0189 
0190     s = (((((((pe[8] * L + pe[7]) * L + pe[6]) * L + pe[5]) * L + pe[4]) * L + pe[3]) * L + pe[2]) * L + pe[1]) * L +
0191         pe[0];
0192 
0193     integrated_flux += 1. / pow(p, 3) * s * corr[k] * (p2 - p1);
0194 
0195     /*

0196     std::cout << k << " " 

0197     << corr[k] << " " 

0198     << p << " " 

0199     << s << " " 

0200     << p1 << " " 

0201     << p2 << " " 

0202     << integrated_flux << " " 

0203     << std::endl;

0204       */
0205 
0206     // std::cout << k << " " << corr[k] << " " << std::endl;

0207   }
0208 
0209   integrated_flux *= 1.27E3;
0210   std::cout << " >>> CMSCGEN.initialize <<< "
0211             << " Integrated flux = " << integrated_flux << " /m**2/s " << std::endl;
0212 
0213   //  find approximate peak value, for Monte Carlo sampling

0214   //      peak is near L = 2

0215 
0216   double ce;
0217 
0218   ce = (((((((pe[8] * 2. + pe[7]) * 2. + pe[6]) * 2. + pe[5]) * 2. + pe[4]) * 2. + pe[3]) * 2. + pe[2]) * 2. + pe[1]) *
0219            2. +
0220        pe[0];
0221 
0222   // normalize to 0.5 (not 1) to have some margin if peak is not at L=2

0223   //

0224   ce = 0.5 / ce;
0225 
0226   for (int k = 0; k < 9; k++) {
0227     pe[k] = pe[k] * ce;
0228   }
0229 
0230   cemax = cemax0 * corr[50];
0231 
0232   return initialization;
0233 }
0234 
0235 int CMSCGEN::initialize(double pmin_in,
0236                         double pmax_in,
0237                         double thetamin_in,
0238                         double thetamax_in,
0239                         int RanSeed,
0240                         bool TIFOnly_constant,
0241                         bool TIFOnly_linear) {
0242   CLHEP::HepRandomEngine *rnd = new CLHEP::HepJamesRandom;
0243   //set seed for Random Generator (seed can be controled by config-file), P.Biallass 2006

0244   rnd->setSeed(RanSeed, 0);
0245   delRanGen = true;
0246   return initialize(pmin_in, pmax_in, thetamin_in, thetamax_in, rnd, TIFOnly_constant, TIFOnly_linear);
0247 }
0248 
0249 int CMSCGEN::generate() {
0250   if (initialization == 0) {
0251     std::cout << " >>> CMSCGEN <<< warning: not initialized" << std::endl;
0252     return -1;
0253   }
0254 
0255   // note: use historical notation (fortran version l3cgen.f)

0256 
0257   //

0258   // +++ determine x = 1/e**2

0259   //

0260   //  explanation: the energy distribution

0261   //        dn/d(1/e**2) = dn/de * e**3 = dn/dlog10(e) * e**2

0262   //     is parametrized by a polynomial. accordingly xe = 1/e**2 is sampled

0263   //     and e calculated

0264   //

0265   //     need precise random variable with high precison since for

0266   //     emin = 3 GeV energies around 3000 GeV are very rare!

0267   //

0268 
0269   double r1, r2, r3;
0270   double xe, e, ce, L, L2;
0271   int k;
0272   double prob;
0273 
0274   double c_max;
0275   double z, z_max;
0276 
0277   while (true) {
0278     prob = RanGen2->flat();
0279     r1 = double(prob);
0280     prob = RanGen2->flat();
0281     r2 = double(prob);
0282 
0283     xe = xemin + r1 * (xemax - xemin);
0284 
0285     if ((1. / sqrt(xe) < 3) &&
0286         TIFOnly_const == true) {  //generate constant energy dependence for E<2GeV, only used for TIF

0287       //compute constant to match to CMSCGEN spectrum

0288       e = 3.;
0289       L = log10(e);
0290       L2 = L * L;
0291 
0292       ce = (((((((pe[8] * L + pe[7]) * L + pe[6]) * L + pe[5]) * L + pe[4]) * L + pe[3]) * L + pe[2]) * L + pe[1]) * L +
0293            pe[0];
0294 
0295       k = int((L - Lmin) * Lfac + 1.);
0296       k = TMath::Max(1, TMath::Min(k, 100));
0297       ce = ce * corr[k];
0298 
0299       e = 1. / sqrt(xe);
0300       if (r2 < (e * e * e * ce / (cemax * 3. * 3. * 3.)))
0301         break;
0302 
0303     } else if ((1. / sqrt(xe) < 3) &&
0304                TIFOnly_lin == true) {  //generate linear energy dependence for E<2GeV, only used for TIF

0305       //compute constant to match to CMSCGEN spectrum

0306       e = 3.;
0307       L = log10(e);
0308       L2 = L * L;
0309 
0310       ce = (((((((pe[8] * L + pe[7]) * L + pe[6]) * L + pe[5]) * L + pe[4]) * L + pe[3]) * L + pe[2]) * L + pe[1]) * L +
0311            pe[0];
0312 
0313       k = int((L - Lmin) * Lfac + 1.);
0314       k = TMath::Max(1, TMath::Min(k, 100));
0315       ce = ce * corr[k];
0316 
0317       e = 1. / sqrt(xe);
0318       if (r2 < (e * e * e * e * ce / (cemax * 3. * 3. * 3. * 3.)))
0319         break;
0320 
0321     } else {  //this is real CMSCGEN energy-dependence

0322 
0323       e = 1. / sqrt(xe);
0324       L = log10(e);
0325       L2 = L * L;
0326 
0327       ce = (((((((pe[8] * L + pe[7]) * L + pe[6]) * L + pe[5]) * L + pe[4]) * L + pe[3]) * L + pe[2]) * L + pe[1]) * L +
0328            pe[0];
0329 
0330       k = int((L - Lmin) * Lfac + 1.);
0331       k = TMath::Max(1, TMath::Min(k, 100));
0332       ce = ce * corr[k];
0333 
0334       if (cemax * r2 < ce)
0335         break;
0336 
0337     }  //end of CMSCGEN energy-dependence

0338   }    //end of while

0339 
0340   pq = e;
0341 
0342   //

0343   // +++ charge ratio 1.280

0344   //

0345   prob = RanGen2->flat();
0346   r3 = double(prob);
0347 
0348   double charg = 1.;
0349   if (r3 < 0.439)
0350     charg = -1.;
0351 
0352   pq = pq * charg;
0353 
0354   //

0355   //  +++ determine cos(angle)

0356   //

0357   //  simple trial and rejection method

0358   //

0359 
0360   // first calculate energy dependent coefficients b_i

0361 
0362   if (TIFOnly_const == true && e < 3.) {  //forTIF (when E<2GeV use angles of 2GeV cosmic)

0363     L = log10(3.);
0364     L2 = L * L;
0365   }
0366   if (TIFOnly_lin == true && e < 3.) {  //forTIF (when E<2GeV use angles of 2GeV cosmic)

0367     L = log10(3.);
0368     L2 = L * L;
0369   }
0370 
0371   b0 = b0c[0] + b0c[1] * L + b0c[2] * L2;
0372   b1 = b1c[0] + b1c[1] * L + b1c[2] * L2;
0373   b2 = b2c[0] + b2c[1] * L + b2c[2] * L2;
0374 
0375   //

0376   // need to know the maximum of z(c)

0377   //

0378   // first calculate c for which c distribution z(c) = maximum

0379   //

0380   // (note: maximum of curve is NOT always at c = -1, but never at c = -0.1)

0381   //

0382 
0383   // try extremal value (from z'(c) = 0), but only if z''(c) < 0

0384   //

0385   // z'(c) = b1 + b2 * c   =>  at c_max = - b1 / (2 b_2) is z'(c) = 0

0386   //

0387   // z''(c) = b2

0388 
0389   c_max = -1.;
0390 
0391   if (b2 < 0.) {
0392     c_max = -0.5 * b1 / b2;
0393     if (c_max < -1.)
0394       c_max = -1.;
0395     if (c_max > -0.1)
0396       c_max = -0.1;
0397   }
0398 
0399   z_max = b0 + b1 * c_max + b2 * c_max * c_max;
0400 
0401   // again cut out explicitly regions of zero flux

0402   double c_cut = -0.42 + L * 0.35;
0403   if (c_cut > cmax)
0404     c_cut = cmax;
0405 
0406   // now we throw dice:

0407 
0408   while (true) {
0409     prob = RanGen2->flat();
0410     r1 = double(prob);
0411     prob = RanGen2->flat();
0412     r2 = double(prob);
0413     c = cmin + (c_cut - cmin) * r1;
0414     z = b0 + b1 * c + b2 * c * c;
0415     if (z > z_max * r2)
0416       break;
0417   }
0418 
0419   return 0;
0420 }
0421 
0422 double CMSCGEN::momentum_times_charge() {
0423   if (initialization == 1) {
0424     return pq;
0425   } else {
0426     std::cout << " >>> CMSCGEN <<< warning: not initialized" << std::endl;
0427     return -9999.;
0428   }
0429 }
0430 
0431 double CMSCGEN::cos_theta() {
0432   if (initialization == 1) {
0433     // here convert between coordinate systems:

0434     return -c;
0435   } else {
0436     std::cout << " >>> CMSCGEN <<< warning: not initialized" << std::endl;
0437     return -0.9999;
0438   }
0439 }
0440 
0441 double CMSCGEN::flux() {
0442   if (initialization == 1) {
0443     return integrated_flux;
0444   } else {
0445     std::cout << " >>> CMSCGEN <<< warning: not initialized" << std::endl;
0446     return -0.9999;
0447   }
0448 }
0449 
0450 int CMSCGEN::initializeNuMu(double pmin_in,
0451                             double pmax_in,
0452                             double thetamin_in,
0453                             double thetamax_in,
0454                             double Enumin_in,
0455                             double Enumax_in,
0456                             double Phimin_in,
0457                             double Phimax_in,
0458                             double ProdAlt_in,
0459                             CLHEP::HepRandomEngine *rnd) {
0460   if (delRanGen)
0461     delete RanGen2;
0462   RanGen2 = rnd;
0463   delRanGen = false;
0464 
0465   ProdAlt = ProdAlt_in;
0466 
0467   Rnunubar = 1.2;
0468 
0469   sigma = (0.72 * Rnunubar + 0.09) / (1 + Rnunubar) * 1.e-38;  //cm^2GeV^-1

0470 
0471   AR = (0.69 + 0.06 * Rnunubar) / (0.09 + 0.72 * Rnunubar);
0472 
0473   //set smin and smax, here convert between coordinate systems:

0474   pmin = pmin_in;
0475   pmax = pmax_in;
0476   cmin = TMath::Cos(thetamin_in);                //input angle already converted from Deg to Rad!

0477   cmax = TMath::Cos(thetamax_in);                //input angle already converted from Deg to Rad!

0478   enumin = (Enumin_in < 10.) ? 10. : Enumin_in;  //no nu's below 10GeV

0479   enumax = Enumax_in;
0480 
0481   //do initial run of flux rate to determine Maximum

0482   integrated_flux = 0.;
0483   dNdEmudEnuMax = 0.;
0484   negabs = 0.;
0485   negfrac = 0.;
0486   int trials = 100000;
0487   for (int i = 0; i < trials; ++i) {
0488     double ctheta = cmin + (cmax - cmin) * RanGen2->flat();
0489     double Emu = pmin + (pmax - pmin) * RanGen2->flat();
0490     double Enu = enumin + (enumax - enumin) * RanGen2->flat();
0491     double rate = dNdEmudEnu(Enu, Emu, ctheta);
0492     //std::cout << "trial=" << i << " ctheta=" << ctheta << " Emu=" << Emu << " Enu=" << Enu

0493     //      << " rate=" << rate << std::endl;

0494     //std::cout << "cmin=" << cmin << " cmax=" << cmax

0495     //      << " pmin=" << pmin << " pmax=" << pmax

0496     //      << " enumin=" << enumin << " enumax=" << enumax << std::endl;

0497     if (rate > 0.) {
0498       integrated_flux += rate;
0499       if (rate > dNdEmudEnuMax)
0500         dNdEmudEnuMax = rate;
0501     } else
0502       negabs++;
0503   }
0504   negfrac = negabs / trials;
0505   integrated_flux /= trials;
0506 
0507   std::cout << "CMSCGEN::initializeNuMu: After " << trials << " trials:" << std::endl;
0508   std::cout << "dNdEmudEnuMax=" << dNdEmudEnuMax << std::endl;
0509   std::cout << "negfrac=" << negfrac << std::endl;
0510 
0511   //multiply by phase space boundaries

0512   integrated_flux *= (cmin - cmax);
0513   integrated_flux *= (Phimax_in - Phimin_in);
0514   integrated_flux *= (pmax - pmin);
0515   integrated_flux *= (enumax - enumin);
0516   //remove negative phase space areas which do not contribute anything

0517   integrated_flux *= (1. - negfrac);
0518   std::cout << " >>> CMSCGEN.initializeNuMu <<< "
0519             << " Integrated flux = " << integrated_flux << " units??? " << std::endl;
0520 
0521   initialization = 1;
0522 
0523   return initialization;
0524 }
0525 
0526 int CMSCGEN::initializeNuMu(double pmin_in,
0527                             double pmax_in,
0528                             double thetamin_in,
0529                             double thetamax_in,
0530                             double Enumin_in,
0531                             double Enumax_in,
0532                             double Phimin_in,
0533                             double Phimax_in,
0534                             double ProdAlt_in,
0535                             int RanSeed) {
0536   CLHEP::HepRandomEngine *rnd = new CLHEP::HepJamesRandom;
0537   //set seed for Random Generator (seed can be controled by config-file), P.Biallass 2006

0538   rnd->setSeed(RanSeed, 0);
0539   delRanGen = true;
0540   return initializeNuMu(
0541       pmin_in, pmax_in, thetamin_in, thetamax_in, Enumin_in, Enumax_in, Phimin_in, Phimax_in, ProdAlt_in, rnd);
0542 }
0543 
0544 double CMSCGEN::dNdEmudEnu(double Enu, double Emu, double ctheta) {
0545   double cthetaNu = 1. + ctheta;  //swap cos(theta) from down to up range

0546   double thetas = asin(sin(acos(cthetaNu)) * (Rearth - SurfaceOfEarth) / (Rearth + ProdAlt));
0547   double costhetas = cos(thetas);
0548   double dNdEnudW =
0549       0.0286 * pow(Enu, -2.7) *
0550       (1. / (1. + (6. * Enu * costhetas) / 115.) + 0.213 / (1. + (1.44 * Enu * costhetas) / 850.));  //cm^2*s*sr*GeV

0551   double dNdEmudEnu = N_A * sigma / alpha * dNdEnudW * 1. / (1. + Emu / epsilon) *
0552                       (Enu - Emu + AR / 3 * (Enu * Enu * Enu - Emu * Emu * Emu) / (Enu * Enu));
0553   return dNdEmudEnu;
0554 }
0555 
0556 int CMSCGEN::generateNuMu() {
0557   if (initialization == 0) {
0558     std::cout << " >>> CMSCGEN <<< warning: not initialized" << std::endl;
0559     return -1;
0560   }
0561 
0562   double ctheta, Emu;
0563   while (true) {
0564     ctheta = cmin + (cmax - cmin) * RanGen2->flat();
0565     Emu = pmin + (pmax - pmin) * RanGen2->flat();
0566     double Enu = enumin + (enumax - enumin) * RanGen2->flat();
0567     double rate = dNdEmudEnu(Enu, Emu, ctheta);
0568     if (rate > dNdEmudEnuMax * RanGen2->flat())
0569       break;
0570   }
0571 
0572   c = -ctheta;  //historical sign convention

0573 
0574   pq = Emu;
0575   //

0576   // +++ nu/nubar ratio (~1.2)

0577   //

0578   double charg = 1.;  //nubar -> mu+

0579   if (RanGen2->flat() > Rnunubar / (1. + Rnunubar))
0580     charg = -1.;  //neutrino -> mu-

0581 
0582   pq = pq * charg;
0583 
0584   //int flux += this event rate

0585 
0586   return 1;
0587 }