00001
00002 #include <iostream>
00003 #include <iomanip>
00004 #include <math.h>
00005 #include "CosmicGenerator/CosmicGun.h"
00006
00007
00008
00009
00010
00011 extern "C"{
00012 void cosmic2_(void);
00013 void cosipr_(void);
00014 void cosgin_(void);
00015 void cosgen_(float* emin, float* emax, int* iacc);
00016 }
00017
00018
00019
00020
00021 struct genpar {
00022 float LBINWID, LEMIN, LEMAX;
00023 int NBIN;
00024 float PROBE[100];
00025 };
00026 genpar genpar_;
00027
00028 struct coscut {
00029 float ctcut;
00030 };
00031 coscut coscut_;
00032
00033 struct cosevt{
00034 float ENER, COSTH, PHI, CHRG;
00035 };
00036 cosevt cosevt_;
00037
00038 struct flxout{
00039 float FLUX, FLUX2;
00040 };
00041 flxout flxout_;
00042
00043
00044
00045
00046
00047 CosmicGun* CosmicGun::mpointer = 0;
00048
00049 CosmicGun* CosmicGun::GetCosmicGun(void){
00050 if(!mpointer) mpointer = new CosmicGun();
00051 return mpointer;
00052 }
00053
00054
00055
00056
00057
00058
00059
00060 CosmicGun::CosmicGun(void){
00061 m_event = 0;
00062 m_emin = 50;
00063 m_emax = 500;
00064 m_coscut = 0.35;
00065 m_printevt = 20;
00066 m_printmod = 50;
00067
00068 coscut_.ctcut = m_coscut;
00069 genpar_.LEMIN = log10(m_emin);
00070 genpar_.LEMAX = log10(m_emax);
00071 genpar_.NBIN = 100;
00072 genpar_.LBINWID = (genpar_.LEMAX-genpar_.LEMIN)/genpar_.NBIN;
00073
00074
00075
00076
00077 }
00078
00079
00080
00081 float CosmicGun::InitializeGenerator() {
00082 std::cout << " CosmicGun::InitializeGenerator: E(min,max)=(" << m_emin << "," << m_emax
00083 << ") GeV, and cos(ThetaCut)="<< m_coscut << std::endl;
00084 cosipr_();
00085 cosgin_();
00086 cosmic2_();
00087 return flxout_.FLUX2;
00088 }
00089
00090 void CosmicGun::PrintLevel(int printevt,int printmod){
00091 if (printevt >= 0)
00092 {
00093 m_printevt = printevt;
00094 }
00095 else
00096 {
00097 std::cerr << "CosmicGun::PrintLevel - warning ignored input printevt = " << printevt << std::endl;
00098 }
00099 if (printmod >= 1)
00100 {
00101 m_printmod = printmod;
00102 }
00103 else
00104 {
00105 std::cerr << "CosmicGun::PrintLevel - warning ignored input printmod = " << printmod << std::endl;
00106 }
00107 }
00108
00109 HepLorentzVector CosmicGun::GenerateEvent(void){
00110 int iacc = 0;
00111
00112 while(iacc == 0){
00113 cosgen_(&m_emin, &m_emax, &iacc);
00114 }
00115 m_event++;
00116
00117 float sinth = sqrt( 1-pow(cosevt_.COSTH,2) );
00118 float e = cosevt_.ENER;
00119 float px = cosevt_.ENER * sinth * sin( cosevt_.PHI);
00120 float py = cosevt_.ENER * sinth * cos( cosevt_.PHI);
00121 float pz = cosevt_.ENER * cosevt_.COSTH;
00122 HepLorentzVector p(px,py,pz,e);
00123
00124
00125 if(m_event < m_printevt)
00126 {
00127 std::cout << "CosmicGun::GenerateEvent: " << std::setw(4) << m_event
00128 << " muon charge " << std::setw(2) << cosevt_.CHRG << " with momentum : " << p << std::endl;
00129 }
00130
00131 return p;
00132 }
00133
00134 int CosmicGun::GetMuonCharge(void){
00135 return (int)cosevt_.CHRG;
00136 }
00137
00138 void CosmicGun::SetEnergyRange(float emin, float emax){
00139 if(emin >= emax || emin < 0 )
00140 {
00141 std::cout << "Error input energy range : (" << emin << " - " << emax << ") - ignored " << std::endl;
00142 return;
00143 }
00144 m_emin = emin;
00145 m_emax = emax;
00146
00147 genpar_.LEMIN = log10(m_emin);
00148 genpar_.LEMAX = log10(m_emax);
00149 genpar_.NBIN = 100;
00150 genpar_.LBINWID = (genpar_.LEMAX-genpar_.LEMIN)/genpar_.NBIN;
00151
00152 }
00153
00154 void CosmicGun::SetCosCut(float ctcut){
00155 m_coscut = ctcut;
00156
00157 coscut_.ctcut = m_coscut;
00158 }
00159