CosmicGun.cxx

Go to the documentation of this file.

 
#include <iostream>
#include <iomanip>
#include <math.h>
#include "CosmicGenerator/CosmicGun.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;
}