CosmicGun.cxx

Go to the documentation of this file.

 
#include "CosmicGun.h"
#include <iomanip>
#include <iostream>
#include <math.h>
 
//--------------------------------------------------------
//
// These are the fortran subroutines
//
extern "C" {
void cosmic2_( void );
void cosipr_( void );
void cosgin_( void );
void cosgen_( float* emin, float* emax, int* iacc );
}
//--------------------------------------------------------
//
// and the common blocks
//
struct genpar {
  float LBINWID, LEMIN, LEMAX;
  int   NBIN;
  float PROBE[100];
};
genpar genpar_;
 
struct coscut {
  float ctcut;
};
coscut coscut_;
 
struct cosevt {
  float ENER, COSTH, PHI, CHRG;
};
cosevt cosevt_;
 
struct flxout {
  float FLUX, FLUX2;
};
flxout flxout_;
 
//--------------------------------------------------------
//
// singleton pattern
//
CosmicGun* CosmicGun::mpointer = 0;
 
CosmicGun* CosmicGun::GetCosmicGun( void ) {
  if ( !mpointer ) mpointer = new CosmicGun();
  return mpointer;
}
 
//--------------------------------------------------------
//
// constructor with some default settings
//
// 9 August 2005, RMcP: remove initialization from constructor.
CosmicGun::CosmicGun( void ) {
  m_event    = 0;
  m_emin     = 50;
  m_emax     = 500;
  m_coscut   = 0.35;
  m_printevt = 20;
  m_printmod = 50;
 
  coscut_.ctcut   = m_coscut;
  genpar_.LEMIN   = log10( m_emin );
  genpar_.LEMAX   = log10( m_emax );
  genpar_.NBIN    = 100;
  genpar_.LBINWID = ( genpar_.LEMAX - genpar_.LEMIN ) / genpar_.NBIN;
 
  // cosipr_();
  // cosgin_();
}
 
// Add separate generator initialization routine
// to avoid forced default initialization, RMcP 9 Aug 05.
float CosmicGun::InitializeGenerator() {
  std::cout << " CosmicGun::InitializeGenerator: E(min,max)=(" << m_emin << "," << m_emax
            << ") GeV, and cos(ThetaCut)=" << m_coscut << std::endl;
  cosipr_();
  cosgin_();
  cosmic2_();
  return flxout_.FLUX2;
}
 
void CosmicGun::PrintLevel( int printevt, int printmod ) {
  if ( printevt >= 0 ) { m_printevt = printevt; }
  else
  {
    std::cerr << "CosmicGun::PrintLevel - warning ignored input printevt = " << printevt
              << std::endl;
  }
  if ( printmod >= 1 ) { m_printmod = printmod; }
  else
  {
    std::cerr << "CosmicGun::PrintLevel - warning ignored input printmod = " << printmod
              << std::endl;
  }
}
 
HepLorentzVector CosmicGun::GenerateEvent( void ) {
  int iacc = 0;
 
  while ( iacc == 0 ) { cosgen_( &m_emin, &m_emax, &iacc ); }
  m_event++;
 
  float            sinth = sqrt( 1 - pow( cosevt_.COSTH, 2 ) );
  float            e     = cosevt_.ENER;
  float            px    = cosevt_.ENER * sinth * sin( cosevt_.PHI );
  float            py    = cosevt_.ENER * sinth * cos( cosevt_.PHI );
  float            pz    = cosevt_.ENER * cosevt_.COSTH;
  HepLorentzVector p( px, py, pz, e );
 
  //  if(m_event < m_printevt || m_event%m_printmod == 0)
  if ( m_event < m_printevt )
  {
    std::cout << "CosmicGun::GenerateEvent: " << std::setw( 4 ) << m_event << " muon charge "
              << std::setw( 2 ) << cosevt_.CHRG << " with momentum : " << p << std::endl;
  }
 
  return p;
}
 
int CosmicGun::GetMuonCharge( void ) { return (int)cosevt_.CHRG; }
 
void CosmicGun::SetEnergyRange( float emin, float emax ) {
  if ( emin >= emax || emin < 0 )
  {
    std::cout << "Error input energy range : (" << emin << " - " << emax << ") - ignored "
              << std::endl;
    return;
  }
  m_emin = emin;
  m_emax = emax;
 
  genpar_.LEMIN   = log10( m_emin );
  genpar_.LEMAX   = log10( m_emax );
  genpar_.NBIN    = 100;
  genpar_.LBINWID = ( genpar_.LEMAX - genpar_.LEMIN ) / genpar_.NBIN;
}
 
void CosmicGun::SetCosCut( float ctcut ) {
  m_coscut = ctcut;
 
  coscut_.ctcut = m_coscut;
}