---

CMSCGEN Class Reference

#include <GeneratorInterface/CosmicMuonGenerator/interface/CMSCGEN.h>

List of all members.

Public Member Functions
	CMSCGEN ()
double	cos_theta ()
double	flux ()
int	generate ()
int	initialize (double, double, double, double, int, bool, bool)
double	momentum_times_charge ()
	~CMSCGEN ()
Private Attributes
double	b0
double	b0c [3]
double	b1
double	b1c [3]
double	b2
double	b2c [3]
double	c
double	c1
double	c2
double	cemax
double	cmax
double	cmax_in
double	cmax_max
double	cmax_min
double	cmin
double	cmin_in
double	corr [101]
int	initialization
double	integrated_flux
double	Lfac
double	Lmax
double	Lmin
double	pe [9]
double	pmax
double	pmin
double	pmin_max
double	pmin_min
double	pq
TRandom2	RanGen2
bool	TIFOnly_const
bool	TIFOnly_lin
double	xemax
double	xemin

---

Detailed Description

Definition at line 29 of file CMSCGEN.h.

---

Constructor & Destructor Documentation

CMSCGEN::CMSCGEN

(

)

[inline]

Definition at line 86 of file CMSCGEN.h.

References initialization.

           {
  initialization = 0;
}

CMSCGEN::~CMSCGEN

(

)

[inline]

Definition at line 91 of file CMSCGEN.h.

References initialization.

            {
  initialization = 0;
}

---

Member Function Documentation

double CMSCGEN::cos_theta

(

)

Definition at line 441 of file CMSCGEN.cc.

References c, GenMuonPlsPt100GeV_cfg::cout, lat::endl(), and initialization.

Referenced by CosmicMuonGenerator::nextEvent().

{
  
  if(initialization==1)
    { 
      // here convert between coordinate systems:
      return -c;
    }
  else
    {
      std::cout << " >>> CMSCGEN <<< warning: not initialized" << std::endl;
      return -0.9999;
    }

}

double CMSCGEN::flux

(

)

Definition at line 457 of file CMSCGEN.cc.

References GenMuonPlsPt100GeV_cfg::cout, lat::endl(), initialization, and integrated_flux.

{
  
  if(initialization==1)
    { 
      return integrated_flux;
    }
  else
    {
      std::cout << " >>> CMSCGEN <<< warning: not initialized" << std::endl;
      return -0.9999;
    }
  
}

int CMSCGEN::generate

(

)

Definition at line 235 of file CMSCGEN.cc.

References b0, b0c, b1, b1c, b2, b2c, c, ce, cemax, cmax, cmin, corr, GenMuonPlsPt100GeV_cfg::cout, e, lat::endl(), initialization, int, k, L, Lfac, Lmin, pe, pq, r1, r2, r3, RanGen2, funct::sqrt(), TIFOnly_const, TIFOnly_lin, xemax, xemin, and z.

Referenced by CosmicMuonGenerator::nextEvent().

{

  if(initialization==0)
    {
      std::cout << " >>> CMSCGEN <<< warning: not initialized" << std::endl;
      return -1;
    }

  // note: use historical notation (fortran version l3cgen.f)

  //
  // +++ determine x = 1/e**2
  //
  //  explanation: the energy distribution 
  //        dn/d(1/e**2) = dn/de * e**3 = dn/dlog10(e) * e**2
  //     is parametrized by a polynomial. accordingly xe = 1/e**2 is sampled
  //     and e calculated 
  //
  //     need precise random variable with high precison since for 
  //     emin = 3 GeV energies around 3000 GeV are very rare!
  //

  double r1, r2, r3;
  double xe, e, ce, L, L2;
  int k;    
  double prob; 
  
  double c_max;
  double z, z_max;
  
  while (1)
    {
      prob = RanGen2.Rndm();
      r1 = double(prob);
      prob = RanGen2.Rndm();
      r2 = double(prob);
      
      xe = xemin+r1*(xemax-xemin);
      
      if( (1./sqrt(xe)<3) && TIFOnly_const == true) { //generate constant energy dependence for E<2GeV, only used for TIF
        //compute constant to match to CMSCGEN spectrum
        e=3.;      
        L = log10(e);
        L2 = L*L;
        
        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
          +pe[0];
        
        k = int ((L-Lmin)*Lfac+1.);
        k = TMath::Max(1,TMath::Min(k,100));
        ce = ce * corr[k];
        
        e = 1./sqrt(xe);  
        if(r2 < ( e*e*e*ce/(cemax*3.*3.*3.) )) break;
      
      }else if( (1./sqrt(xe)<3) && TIFOnly_lin == true) { //generate linear energy dependence for E<2GeV, only used for TIF
        //compute constant to match to CMSCGEN spectrum
        e=3.;      
        L = log10(e);
        L2 = L*L;
        
        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
          +pe[0];
      
        k = int ((L-Lmin)*Lfac+1.);
        k = TMath::Max(1,TMath::Min(k,100));
        ce = ce * corr[k];
      
        e = 1./sqrt(xe);  
        if(r2 < ( e*e*e*e*ce/(cemax*3.*3.*3.*3.) )) break;
      
      }else{ //this is real CMSCGEN energy-dependence

        e = 1./sqrt(xe);       
        L = log10(e);
        L2 = L*L;
        
        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
          +pe[0];
        
        k = int ((L-Lmin)*Lfac+1.);
        k = TMath::Max(1,TMath::Min(k,100));
        ce = ce * corr[k];
        
        if(cemax*r2 < ce) break;

      } //end of CMSCGEN energy-dependence
    } //end of while
  
  pq = e;
  
  //
  // +++ charge ratio 1.280
  //
  prob = RanGen2.Rndm();
  r3 = double(prob);
  
  double charg = 1.;
  if(r3 < 0.439) charg=-1.;
  
  pq = pq*charg;
  
  //
  //  +++ determine cos(angle)
  //
  //  simple trial and rejection method
  //
  
  // first calculate energy dependent coefficients b_i

  if(TIFOnly_const == true && e<3.){ //forTIF (when E<2GeV use angles of 2GeV cosmic)
    L = log10(3.);
    L2 = L*L;
  }
  if(TIFOnly_lin == true && e<3.){ //forTIF (when E<2GeV use angles of 2GeV cosmic)
    L = log10(3.);
    L2 = L*L; 
  }
  
  b0 = b0c[0] + b0c[1] * L + b0c[2]* L2;
  b1 = b1c[0] + b1c[1] * L + b1c[2]* L2;
  b2 = b2c[0] + b2c[1] * L + b2c[2]* L2;
  
  //
  // need to know the maximum of z(c)
  //
  // first calculate c for which c distribution z(c) = maximum 
  // 
  // (note: maximum of curve is NOT always at c = -1, but never at c = -0.1)   
  //
  
  // try extremal value (from z'(c) = 0), but only if z''(c) < 0  
  //
  // z'(c) = b1 + b2 * c   =>  at c_max = - b1 / (2 b_2) is z'(c) = 0
  //
  // z''(c) = b2 
  
  c_max = -1.;
  
  if(b2<0.) {
    c_max = - 0.5 * b1/b2;
    if(c_max < -1.) c_max = -1.;
    if(c_max > -0.1) c_max = -0.1;
  } 
  
  z_max = b0 + b1 * c_max + b2 * c_max * c_max;

  // again cut out explicitly regions of zero flux 
  double c_cut = -0.42 + L*0.35;
  if (c_cut > cmax) c_cut = cmax; 
  
  // now we throw dice:
  
  while (1)
    {
      prob = RanGen2.Rndm();
      r1 = double(prob);
      prob = RanGen2.Rndm();
      r2 = double(prob);
      c = cmin + (c_cut-cmin)*r1;    
      z = b0 + b1 * c + b2 * c*c;
      if (z > z_max*r2) break;
    }    
  
  return 0;
  
}

int CMSCGEN::initialize	(	double	pmin_in,
		double	pmax_in,
		double	thetamin_in,
		double	thetamax_in,
		int	RanSeed,
		bool	TIFOnly_constant,
		bool	TIFOnly_linear
	)

Definition at line 11 of file CMSCGEN.cc.

References b0, b0c, b1, b1c, b2, b2c, c1, c2, ce, cemax, cmax, cmax_in, cmax_max, cmax_min, cmin, cmin_in, corr, GenMuonPlsPt100GeV_cfg::cout, lat::endl(), initialization, integrated_flux, k, L, Lfac, Lmax, Lmin, p, p1, p2, pe, pmax, pmin, pmin_max, pmin_min, funct::pow(), RanGen2, s, TIFOnly_const, TIFOnly_lin, xemax, and xemin.

Referenced by CosmicMuonGenerator::initialize().

{

  //set seed for Random Generator (seed can be controled by config-file), P.Biallass 2006
  RanGen2.SetSeed(RanSeed);

  //set bools for TIFOnly options (E<2GeV with unphysical energy dependence)
  TIFOnly_const = TIFOnly_constant;
  TIFOnly_lin = TIFOnly_linear;


  // units: GeV

  // WARNING: coordinate system: 
  //   - to outside world define z axis downwards, i.e.
  //          muon coming from above, vertically: cos = 1
  //      (used for backward compatibility)
  //   - internally use frame with z axis upwards, i.e.
  //          muon coming from above, vertically: cos = -1
  //     (corresponds to CMS note definition)

  //set cmin and cmax, here convert between coordinate systems:
  cmin_in = - TMath::Cos(thetamin_in);//input angle already converted from Deg to Rad!
  cmax_in = - TMath::Cos(thetamax_in);//input angle already converted from Deg to Rad!


  //allowed energy range
  pmin_min = 3.;
  //pmin_max = 100.;
  pmin_max = 3000.;
  pmax = 3000.;
  //allowed angular range
  //cmax_max = -0.1,
  cmax_max = -0.01,
  cmax_min = -0.9999;

 if(TIFOnly_const == true || TIFOnly_lin == true) pmin_min = 0.; //forTIF

  // set pmin  
  if(pmin_in < pmin_min || pmin_in > pmin_max){ 
    std::cout << " >>> CMSCGEN.initialize <<< warning: illegal pmin_in =" << pmin_in;
    return(-1);
  } else if(pmax_in > pmax ){
    std::cout << " >>> CMSCGEN.initialize <<< warning: illegal pmax_in =" << pmax_in;
    return(-1);
  }else{     
    pmin = pmin_in;
    pmax = pmax_in;
    xemax = 1./(pmin*pmin);
    xemin = 1./(pmax*pmax); 
  }


  // set cmax and cmin
  if(cmax_in < cmax_min || cmax_in > cmax_max)
    { 
      std::cout << " >>> CMSCGEN.initialize <<< warning: illegal cmax_in =" << cmax_in;
      return(-1);
    }
  else 
    {
      cmax = cmax_in;
      cmin = cmin_in;
    }


  initialization = 1;

  if(TIFOnly_const == true || TIFOnly_lin == true) pmin_min = 3.; //forTIF

  //  Lmin = log10(pmin_min);
  Lmin = log10(pmin);
  Lmax = log10(pmax);
  Lfac = 100./(Lmax-Lmin);

  //
  // +++ coefficients for energy spectrum
  //

  pe[0] = -1.;
  pe[1] = 6.22176;
  pe[2] = -13.9404;
  pe[3] = 18.1643;
  pe[4] = -9.22784;
  pe[5] = 1.99234;
  pe[6] = -0.156434;
  pe[7] = 0.;
  pe[8] = 0.;

  //
  // +++ coefficients for cos theta distribution
  //

  b0c[0] = 0.6639;
  b0c[1] = -0.9587;
  b0c[2] = 0.2772;
  
  b1c[0] = 5.820;
  b1c[1] = -6.864;
  b1c[2] = 1.367;
  
  b2c[0] = 10.39;
  b2c[1] = -8.593;
  b2c[2] = 1.547;
  
  //
  // +++ calculate correction table for different cos theta dependence!
  //     reference range: c1  to  c2 
  //
  //  explanation: the parametrization of the energy spectrum as used above
  //    is the integral over the c = cos(zenith angle) range -1...-0.1 
  //    since the c distribution depends on energy, the integrated energy
  //    spectrum depends on this range. Here a correction factor is determined,
  //    based on the linear c dependence of the c distribution.
  //    The correction is calculated for 100 bins in L = log10(energy).
  //
  // +++ in same loop calculate integrated flux 
  //      (integrated over angles and momentum)

  c1 = -1.;
  c2 = -0.1;
  
  double cemax0 = 1.0;
  double L, L2; 
  double s;
  double p, p1, p2;
  double integral_here, integral_ref;
  double c_cut;

  integrated_flux = 0.;
  
  for(int k=1; k<=100; k++)
    {
      L = Lmin + (k-0.5)/Lfac;
      L2 = L*L;   
      p = pow(10,L);
      p1 = pow(10,L-0.5/Lfac);
      p2 = pow(10,L+0.5/Lfac);
      
      b0 = b0c[0] + b0c[1] * L + b0c[2]* L2;
      b1 = b1c[0] + b1c[1] * L + b1c[2]* L2;
      b2 = b2c[0] + b2c[1] * L + b2c[2]* L2;

      // cut out explicitly regions of zero flux 
      // (for low momentum and near horizontal showers)
      // since parametrization for z distribution doesn't work here
      // (can become negative) 
      
      c_cut = -0.42 + L*0.35;
      
      if (c_cut > c2) c_cut = c2; 
      
      integral_ref = b0 * (c_cut - c1) 
        + b1/2. * (c_cut*c_cut - c1*c1) 
        + b2/3. * (c_cut*c_cut*c_cut - c1*c1*c1);
      
      if (c_cut > cmax) c_cut = cmax; 
      
      integral_here = b0 * (c_cut - cmin) 
        + b1/2. * (c_cut*c_cut - cmin*cmin) 
        + b2/3. * (c_cut*c_cut*c_cut - cmin*cmin*cmin);
  
      corr[k] = integral_here/integral_ref;
      
      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
        +pe[0];
    
      integrated_flux += 1./pow(p,3) * s * corr[k] * (p2-p1);

      /*
        std::cout << k << " " 
        << corr[k] << " " 
        << p << " " 
        << s << " " 
        << p1 << " " 
        << p2 << " " 
        << integrated_flux << " " 
        << std::endl;
      */
      
      // std::cout << k << " " << corr[k] << " " << std::endl;
    }

  integrated_flux *= 1.27E3;
  std::cout << " >>> CMSCGEN.initialize <<< " <<
    " Integrated flux = " << integrated_flux << " /m**2/s " << std::endl;
  
  //  find approximate peak value, for Monte Carlo sampling
  //      peak is near L = 2 

  double ce;
  
  ce = (((((((pe[8]*2.
              +pe[7])*2.
             +pe[6])*2.
            +pe[5])*2.
           +pe[4])*2.
          +pe[3])*2.
         +pe[2])*2.
        +pe[1])*2.
    +pe[0];

  // normalize to 0.5 (not 1) to have some margin if peak is not at L=2
  //
  ce = 0.5/ce;

  for (int k=0; k<9; k++)
    {
      pe[k] = pe[k]*ce;
    }

  cemax = cemax0*corr[50];      

  return initialization;
}

double CMSCGEN::momentum_times_charge

(

)

Definition at line 426 of file CMSCGEN.cc.

References GenMuonPlsPt100GeV_cfg::cout, lat::endl(), initialization, and pq.

Referenced by CosmicMuonGenerator::nextEvent().

{
  
  if(initialization==1)
    { 
      return pq;
    }
  else
    {
      std::cout << " >>> CMSCGEN <<< warning: not initialized" << std::endl;
      return -9999.;
    }

}

---

Member Data Documentation

double CMSCGEN::b0 [private]

Definition at line 64 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::b0c[3] [private]

Definition at line 73 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::b1 [private]

Definition at line 65 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::b1c[3] [private]

Definition at line 73 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::b2 [private]

Definition at line 66 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::b2c[3] [private]

Definition at line 73 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::c [private]

Definition at line 46 of file CMSCGEN.h.

Referenced by cos_theta(), and generate().

double CMSCGEN::c1 [private]

Definition at line 61 of file CMSCGEN.h.

Referenced by initialize().

double CMSCGEN::c2 [private]

Definition at line 62 of file CMSCGEN.h.

Referenced by initialize().

double CMSCGEN::cemax [private]

Definition at line 70 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::cmax [private]

Definition at line 41 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::cmax_in [private]

Definition at line 43 of file CMSCGEN.h.

Referenced by initialize().

double CMSCGEN::cmax_max [private]

Definition at line 55 of file CMSCGEN.h.

Referenced by initialize().

double CMSCGEN::cmax_min [private]

Definition at line 54 of file CMSCGEN.h.

Referenced by initialize().

double CMSCGEN::cmin [private]

Definition at line 40 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::cmin_in [private]

Definition at line 42 of file CMSCGEN.h.

Referenced by initialize().

double CMSCGEN::corr[101] [private]

Definition at line 74 of file CMSCGEN.h.

Referenced by generate(), and initialize().

int CMSCGEN::initialization [private]

Definition at line 36 of file CMSCGEN.h.

Referenced by CMSCGEN(), cos_theta(), flux(), generate(), initialize(), momentum_times_charge(), and ~CMSCGEN().

double CMSCGEN::integrated_flux [private]

Definition at line 68 of file CMSCGEN.h.

Referenced by flux(), and initialize().

double CMSCGEN::Lfac [private]

Definition at line 59 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::Lmax [private]

Definition at line 58 of file CMSCGEN.h.

Referenced by initialize().

double CMSCGEN::Lmin [private]

Definition at line 57 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::pe[9] [private]

Definition at line 72 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::pmax [private]

Definition at line 39 of file CMSCGEN.h.

Referenced by initialize().

double CMSCGEN::pmin [private]

Definition at line 38 of file CMSCGEN.h.

Referenced by initialize().

double CMSCGEN::pmin_max [private]

Definition at line 52 of file CMSCGEN.h.

Referenced by initialize().

double CMSCGEN::pmin_min [private]

Definition at line 51 of file CMSCGEN.h.

Referenced by initialize().

double CMSCGEN::pq [private]

Definition at line 45 of file CMSCGEN.h.

Referenced by generate(), and momentum_times_charge().

TRandom2 CMSCGEN::RanGen2 [private]

Definition at line 77 of file CMSCGEN.h.

Referenced by generate(), and initialize().

bool CMSCGEN::TIFOnly_const [private]

Definition at line 80 of file CMSCGEN.h.

Referenced by generate(), and initialize().

bool CMSCGEN::TIFOnly_lin [private]

Definition at line 81 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::xemax [private]

Definition at line 49 of file CMSCGEN.h.

Referenced by generate(), and initialize().

double CMSCGEN::xemin [private]

Definition at line 48 of file CMSCGEN.h.

Referenced by generate(), and initialize().

---

The documentation for this class was generated from the following files:

• GeneratorInterface/CosmicMuonGenerator/interface/CMSCGEN.h
• GeneratorInterface/CosmicMuonGenerator/src/CMSCGEN.cc

---