00001 //---------------------------------------------------------------------------// 00002 //bb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nu// 00003 // // 00004 // MAJORANA Simulation // 00005 // // 00006 // This code implementation is the intellectual property of the // 00007 // MAJORANA Collaboration. It is based on Geant4, an intellectual // 00008 // property of the RD44 GEANT4 collaboration. // 00009 // // 00010 // ********************* // 00011 // // 00012 // Neither the authors of this software system, nor their employing // 00013 // institutes, nor the agencies providing financial support for this // 00014 // work make any representation or warranty, express or implied, // 00015 // regarding this software system or assume any liability for its use. // 00016 // By copying, distributing or modifying the Program (or any work based // 00017 // on on the Program) you indicate your acceptance of this statement, // 00018 // and all its terms. // 00019 // // 00020 //bb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nu// 00021 //---------------------------------------------------------------------------// 00022 // 00023 // $Id: MJGeneratorPNNL.cc,v 1.2 2004/11/09 13:42:39 xliu Exp $ 00024 // 00025 // CLASS IMPLEMENTATION: MJGeneratorPNNL.cc 00026 // 00027 //---------------------------------------------------------------------------// 00033 // 00034 //---------------------------------------------------------------------------// 00044 //---------------------------------------------------------------------------// 00045 // 00046 00047 00048 #include <math.h> 00049 00050 #include <CLHEP/Units/PhysicalConstants.h> 00051 00052 #include "globals.hh" 00053 #include "Randomize.hh" 00054 #include "G4Event.hh" 00055 #include "G4Gamma.hh" 00056 #include "G4ParticleGun.hh" 00057 #include "G4Run.hh" 00058 #include "G4ThreeVector.hh" 00059 00060 #include "generators/MJGeneratorPNNLMessenger.hh" 00061 #include "generators/MJGeneratorPNNLRadioisotope.hh" 00062 #include "io/MJLogger.hh" 00063 00064 //---------------------------------------------------------------------------// 00065 00066 #include "generators/MJGeneratorPNNL.hh" 00067 00068 //---------------------------------------------------------------------------// 00069 00070 MJGeneratorPNNL::MJGeneratorPNNL(): 00071 fParticleGun(0) 00072 { 00073 fGeneratorName = "PNNLiso"; 00074 fG4Messenger = new MJGeneratorPNNLMessenger(this); 00075 fPosition = G4ThreeVector(0.0, 0.0, 0.0); 00076 } 00077 00078 //---------------------------------------------------------------------------// 00079 00080 MJGeneratorPNNL::MJGeneratorPNNL(const MJGeneratorPNNL & other) 00081 {;} 00082 00083 //---------------------------------------------------------------------------// 00084 00085 MJGeneratorPNNL::~MJGeneratorPNNL() 00086 { 00087 delete fParticleGun; 00088 delete fG4Messenger; 00089 } 00090 00091 //---------------------------------------------------------------------------// 00092 00093 void MJGeneratorPNNL::BeginOfRunAction(G4Run const *run) 00094 { 00095 for (G4int i=0; i<10; i++) 00096 hist_isotope[i] = 0; 00097 G4int n_particle = 1; 00098 fParticleGun = new G4ParticleGun(n_particle); 00099 fFirstCall = true; 00100 fUsePNNLGen = false; 00101 MJLog(trace) << " called." << endlog; 00102 } 00103 00104 //---------------------------------------------------------------------------// 00105 00106 void MJGeneratorPNNL::EndOfRunAction(G4Run const *run) 00107 { 00108 MJLog(routine) << "Number of times each radioisotope was sampled:"<<endlog; 00109 for (G4int i=0; i<10; i++) 00110 G4cout << i << " " << hist_isotope[i] << G4endl; 00111 MJLog(routine) << endlog; 00112 } 00113 00114 //---------------------------------------------------------------------------// 00115 00116 void MJGeneratorPNNL::GeneratePrimaryVertex(G4Event *event) 00117 { 00118 if (fFirstCall) { 00119 MJLog(routine) << "Called for first time." << endlog; 00120 MJLog(routine) << "fPNNL_source_age = " << fPNNL_source_age << endlog; 00121 MJLog(routine)<<"fPNNL_DecayChain_file = "<<fPNNL_DecayChain_file<<endlog; 00122 00123 f_pPNNLDecayChain = new MJGeneratorPNNLDecayChain(fPNNL_DecayChain_file, 00124 fPNNL_source_age); 00125 fUsePNNLGen = true; 00126 fFirstCall = false; 00127 } 00128 00129 const G4double PI = (G4double) pi; // Hack. RH. Needs improvement. 00130 00131 G4int num_events = event->GetEventID(); 00132 if ((num_events % fReportingFrequency) == 0) 00133 MJLog(routine) << "Doing event " << num_events << endlog; 00134 00135 G4int num_gammas = 0; 00136 G4int indx_isotope; 00137 while (fUsePNNLGen && num_gammas == 0) { 00138 MJLog(debugging) << "Before DoCascadeEvent." << endlog; 00139 MJGeneratorPNNLCascadeEvent CascadeEvent = 00140 f_pPNNLDecayChain->DoCascadeEvent(); 00141 MJLog(debugging) << "Before GetIndex..." << endlog; 00142 indx_isotope = f_pPNNLDecayChain->GetIndexOfLastIsotopeSampled(); 00143 MJLog(debugging) << "hist_isotopes" << endlog; 00144 if ((indx_isotope >= 0) && (indx_isotope < 10)) 00145 hist_isotope[indx_isotope]++; 00146 else { 00147 MJLog(error)<<"PrimaryGeneratorAction: Anomalous sampled isotope index:" 00148 << indx_isotope << " Aborting..." << endlog; 00149 MJLog(fatal) << endlog; 00150 } 00151 00152 G4double E_charged_particle = CascadeEvent.GetChargedParticleE(); 00153 MJLog(debugging)<<"Charged particle energy: "<<E_charged_particle<<endlog; 00154 00155 num_gammas = CascadeEvent.GetNumGammas(); 00156 MJLog(debugging) << "Number of gammas: " << num_gammas << endlog; 00157 00158 if (num_gammas > 0) { 00159 G4double* pGammaList = CascadeEvent.GetGammaList(); 00160 00161 for (G4int igamma=0; igamma<num_gammas; igamma++) { 00162 G4double Energy_gamma = pGammaList[igamma]; 00163 // G4cout << igamma << " " << Energy_gamma << G4endl; 00164 G4double cos_theta = -1.0 + 2.0*G4UniformRand(); 00165 G4double phi = 2.0*PI*G4UniformRand(); 00166 G4double theta = acos(cos_theta); 00167 G4double sin_theta = sin(theta); 00168 G4double cos_x = sin_theta*cos(phi); 00169 G4double cos_y = sin_theta*sin(phi); 00170 G4double cos_z = cos_theta; 00171 00172 fParticleGun->SetParticleDefinition(G4Gamma::GammaDefinition()); 00173 MJLog(debugging) << "Momentum." << cos_x << cos_y << cos_z << endlog; 00174 fParticleGun->SetParticleMomentumDirection( 00175 G4ThreeVector(cos_x,cos_y,cos_z)); 00176 MJLog(debugging) << "Energy." << Energy_gamma << endlog; 00177 fParticleGun->SetParticleEnergy(Energy_gamma); 00178 MJLog(debugging) << "Setting Particle Position." << endlog; 00179 fParticleGun->SetParticlePosition(fPosition); 00180 MJLog(debugging) << "Position : " << fPosition << endlog; 00181 fParticleGun->GeneratePrimaryVertex(event); 00182 } // for (G4int igamma=0; igamma<num_gammas; igamma++) 00183 } // if(num_gammas > 0) 00184 } // while (fUsePNNLGen && num_gammas == 0) 00185 } 00186 00187 //---------------------------------------------------------------------------// 00188 //---------------------------------------------------------------------------//
1.3.9.1