#include <DsG4GdCaptureGammas.h>
Public Member Functions | |
DsG4GdCaptureGammas () | |
~DsG4GdCaptureGammas () | |
G4ReactionProductVector * | GetGammas () |
std::vector< double > | GetEnergy () |
Private Attributes | |
TFile * | specfile |
TH1F * | capspec |
Definition at line 18 of file DsG4GdCaptureGammas.h.
DsG4GdCaptureGammas::DsG4GdCaptureGammas | ( | ) |
Definition at line 25 of file DsG4GdCaptureGammas.cc.
00026 { 00027 //First back up the current TDirectory. This is hard-wired, but at 00028 //this point we assume that an output root file has been opened at the 00029 //very beginning of the run. I will try to find a prettier fix 00030 //later. jianglai 05/09/2006 00031 TDirectory *dirKeep = gDirectory; 00032 00033 // The experimental spectrum is saved in gdcap_spec.root. 00034 // This file should be placed in the data/ directory. 00035 // It will get copied to the InstallArea/${tag}/data directory when 00036 // the package is 'make'd. 00037 // The InstallArea/${tag}/data which should be linked to with the $DAYA_DATA_DIR 00038 // environment variable. 00039 00040 std::string specfilename = "gdcap_spec.root"; 00041 const char* prefix = getenv("DAYA_DATA_DIR"); 00042 if(prefix) { 00043 if(strlen(prefix)>0) { 00044 std::string newname = prefix; 00045 newname += "/"; 00046 newname += specfilename; 00047 specfilename = newname; 00048 } 00049 } 00050 00051 specfile = TFile::Open(specfilename.c_str()); 00052 00053 if(!specfile) 00054 G4cout << "Can not open gdcap_spec.root. " 00055 << "Make sure DAYA_DATA_DIR is set to it's location." << G4endl; 00056 if(!(specfile->IsOpen())) 00057 G4cout << "Can not open gdcap_spec.root. " 00058 << "Make sure DAYA_DATA_DIR is set to it's location." << G4endl; 00059 00060 capspec = (TH1F*)(specfile->Get("h1")); 00061 00062 // return to the TDirectory used before opening TFile 00063 if (dirKeep) { 00064 dirKeep->cd(); 00065 } 00066 }
DsG4GdCaptureGammas::~DsG4GdCaptureGammas | ( | ) |
Definition at line 68 of file DsG4GdCaptureGammas.cc.
00069 { 00070 /* 00071 if(specfile && specfile->IsOpen()) { 00072 specfile->Close(); 00073 G4cout <<"gdcap_spec.root closed"<<G4endl; 00074 } 00075 else 00076 G4cout <<"gdcap_spec.root is not closed."<<G4endl; 00077 */ 00078 }
G4ReactionProductVector * DsG4GdCaptureGammas::GetGammas | ( | ) |
Definition at line 80 of file DsG4GdCaptureGammas.cc.
00081 { 00082 G4ReactionProductVector* theGammas = new G4ReactionProductVector; 00083 vector<double> energy = GetEnergy(); 00084 for(unsigned int i=0; i<energy.size(); i++) { 00085 00086 G4ReactionProduct* theOne = new G4ReactionProduct; 00087 theOne->SetDefinition( G4Gamma::Gamma() ); 00088 00089 // Get the gammas direction. 00090 // Isotropic emission. 00091 G4double costheta = 2.*G4UniformRand()-1; 00092 G4double theta = acos(costheta); 00093 G4double phi = twopi*G4UniformRand(); 00094 G4double sinth = sin(theta); 00095 theOne->SetTotalEnergy( energy[i] ); 00096 G4ThreeVector temp(energy[i]*sinth*cos(phi), 00097 energy[i]*sinth*sin(phi), 00098 energy[i]*cos(theta) ); 00099 theOne->SetMomentum( temp ) ; 00100 theGammas->push_back(theOne); 00101 } 00102 return theGammas; 00103 }
vector< double > DsG4GdCaptureGammas::GetEnergy | ( | ) |
Definition at line 157 of file DsG4GdCaptureGammas.cc.
00158 { 00159 vector<double> Energy; 00160 G4double TotalEnergy; 00161 if (G4UniformRand()<0.815) 00162 TotalEnergy = 7.93; // captured by Gd157 00163 else 00164 TotalEnergy = 8.53; // captured by Gd155 00165 G4double energy_sum = 0.0; 00166 G4int gamma_num = 0; 00167 G4double gamma_energy[100]={0}; 00168 G4double energy_tmp; 00169 00170 sample: 00171 00172 00173 //Use the CLHEP random engine to sample the distribution instead 00174 //Jianglai 10/01/2006 00175 //G4double energy=capspec->GetRandom();// sample from the experimental spectrum. 00176 G4double energy = getRandomCLHEP(capspec); 00177 00178 energy_sum = energy_sum+energy; 00179 gamma_num = gamma_num+1; 00180 gamma_energy[gamma_num] = energy; 00181 00182 if(energy_sum>TotalEnergy) { 00183 if((energy_sum-TotalEnergy)<0.8) { 00184 energy_tmp = gamma_energy[gamma_num-1]+energy-(energy_sum-TotalEnergy); 00185 if(energy_tmp>7.4) { 00186 energy_sum = energy_sum-energy; 00187 gamma_num = gamma_num-1; 00188 goto sample; 00189 } 00190 gamma_num = gamma_num-1; 00191 gamma_energy[gamma_num]=energy_tmp; 00192 goto goon; 00193 } 00194 energy_sum = energy_sum-energy; 00195 gamma_num = gamma_num-1; 00196 goto sample; 00197 } 00198 goto sample; 00199 00200 goon: 00201 00202 // set the sequence 00203 for(int i=gamma_num;i>1;i--) { 00204 for(int j=1;j<i;j++) { 00205 if(gamma_energy[j]<gamma_energy[j+1]) { 00206 energy_tmp = gamma_energy[j+1]; 00207 gamma_energy[j+1] = gamma_energy[j]; 00208 gamma_energy[j] = energy_tmp; 00209 } 00210 } 00211 } 00212 if(gamma_num<=3) 00213 goto ok; 00214 00215 if(gamma_energy[gamma_num]<0.8) { 00216 energy_tmp = gamma_energy[gamma_num]+gamma_energy[gamma_num-1]; 00217 if(energy_tmp>7.4) 00218 goto ok; 00219 gamma_energy[gamma_num-1]=energy_tmp; 00220 gamma_num=gamma_num-1; 00221 } 00222 00223 ok: 00224 00225 for(int k=1;k<=gamma_num;k++) 00226 Energy.push_back(gamma_energy[k]); 00227 00228 return Energy; 00229 }
TFile* DsG4GdCaptureGammas::specfile [private] |
Definition at line 29 of file DsG4GdCaptureGammas.h.
TH1F* DsG4GdCaptureGammas::capspec [private] |
Definition at line 30 of file DsG4GdCaptureGammas.h.