G4DNAChampionElasticModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 
00029 #include "G4DNAChampionElasticModel.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "G4DNAMolecularMaterial.hh"
00033 
00034 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00035 
00036 using namespace std;
00037 
00038 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00039 
00040 G4DNAChampionElasticModel::G4DNAChampionElasticModel(const G4ParticleDefinition*,
00041                                              const G4String& nam)
00042 :G4VEmModel(nam),isInitialised(false)
00043 {
00044 //  nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
00045 
00046   killBelowEnergy = 7.4*eV; 
00047   lowEnergyLimit = 0 * eV; 
00048   highEnergyLimit = 1. * MeV;
00049   SetLowEnergyLimit(lowEnergyLimit);
00050   SetHighEnergyLimit(highEnergyLimit);
00051 
00052   verboseLevel= 0;
00053   // Verbosity scale:
00054   // 0 = nothing 
00055   // 1 = warning for energy non-conservation 
00056   // 2 = details of energy budget
00057   // 3 = calculation of cross sections, file openings, sampling of atoms
00058   // 4 = entering in methods
00059   
00060   if( verboseLevel>0 ) 
00061   { 
00062     G4cout << "Champion Elastic model is constructed " << G4endl
00063            << "Energy range: "
00064            << lowEnergyLimit / eV << " eV - "
00065            << highEnergyLimit / MeV << " MeV"
00066            << G4endl;
00067   }
00068   fParticleChangeForGamma = 0;
00069   fpMolWaterDensity = 0;
00070 }
00071 
00072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00073 
00074 G4DNAChampionElasticModel::~G4DNAChampionElasticModel()
00075 {  
00076   // For total cross section
00077   
00078   std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00079   for (pos = tableData.begin(); pos != tableData.end(); ++pos)
00080   {
00081     G4DNACrossSectionDataSet* table = pos->second;
00082     delete table;
00083   }
00084 
00085    // For final state
00086    
00087    eVecm.clear();
00088 
00089 }
00090 
00091 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00092 
00093 void G4DNAChampionElasticModel::Initialise(const G4ParticleDefinition* /*particle*/,
00094                                        const G4DataVector& /*cuts*/)
00095 {
00096 
00097   if (verboseLevel > 3)
00098     G4cout << "Calling G4DNAChampionElasticModel::Initialise()" << G4endl;
00099 
00100   // Energy limits
00101   
00102   if (LowEnergyLimit() < lowEnergyLimit)
00103   {
00104     G4cout << "G4DNAChampionElasticModel: low energy limit increased from " << 
00105         LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
00106     SetLowEnergyLimit(lowEnergyLimit);
00107     }
00108 
00109   if (HighEnergyLimit() > highEnergyLimit)
00110   {
00111     G4cout << "G4DNAChampionElasticModel: high energy limit decreased from " << 
00112         HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
00113     SetHighEnergyLimit(highEnergyLimit);
00114   }
00115 
00116   // Reading of data files 
00117   
00118   G4double scaleFactor = 1e-16*cm*cm;
00119 
00120   G4String fileElectron("dna/sigma_elastic_e_champion");
00121 
00122   G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
00123   G4String electron;
00124  
00125   // *** ELECTRON
00126   
00127     // For total cross section
00128     
00129     electron = electronDef->GetParticleName();
00130 
00131     tableFile[electron] = fileElectron;
00132 
00133     G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
00134     tableE->LoadData(fileElectron);
00135     tableData[electron] = tableE;
00136     
00137     // For final state
00138 
00139     char *path = getenv("G4LEDATA");
00140  
00141     if (!path)
00142     {
00143       G4Exception("G4ChampionElasticModel::Initialise","em0006",
00144                   FatalException,"G4LEDATA environment variable not set.");
00145       return;
00146     }
00147        
00148     std::ostringstream eFullFileName;
00149     eFullFileName << path << "/dna/sigmadiff_cumulated_elastic_e_champion.dat";
00150     std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
00151      
00152     if (!eDiffCrossSection) 
00153       G4Exception("G4DNAChampionElasticModel::Initialise","em0003",
00154                   FatalException,"Missing data file:/dna/sigmadiff_cumulated_elastic_e_champion.dat");
00155     
00156     eTdummyVec.push_back(0.);
00157 
00158     while(!eDiffCrossSection.eof())
00159     {
00160         double tDummy;
00161         double eDummy;
00162         eDiffCrossSection>>tDummy>>eDummy;
00163         
00164         // SI : mandatory eVecm initialization
00165         
00166         if (tDummy != eTdummyVec.back()) 
00167         { 
00168           eTdummyVec.push_back(tDummy); 
00169           eVecm[tDummy].push_back(0.);
00170         }
00171           
00172         eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
00173 
00174         if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
00175         
00176     }
00177 
00178     // End final state
00179     
00180   if (verboseLevel > 2) 
00181     G4cout << "Loaded cross section files for Champion Elastic model" << G4endl;
00182 
00183   if( verboseLevel>0 ) 
00184   { 
00185     G4cout << "Champion Elastic model is initialized " << G4endl
00186            << "Energy range: "
00187            << LowEnergyLimit() / eV << " eV - "
00188            << HighEnergyLimit() / MeV << " MeV"
00189            << G4endl;
00190   }
00191 
00192   // Initialize water density pointer
00193   fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00194 
00195   if (isInitialised) { return; }
00196   fParticleChangeForGamma = GetParticleChangeForGamma();
00197   isInitialised = true;
00198 
00199 }
00200 
00201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00202 
00203 G4double G4DNAChampionElasticModel::CrossSectionPerVolume(const G4Material* material,
00204                                            const G4ParticleDefinition* p,
00205                                            G4double ekin,
00206                                            G4double,
00207                                            G4double)
00208 {
00209   if (verboseLevel > 3)
00210     G4cout << "Calling CrossSectionPerVolume() of G4DNAChampionElasticModel" << G4endl;
00211 
00212  // Calculate total cross section for model
00213 
00214  G4double sigma=0;
00215 
00216  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
00217 
00218  if(waterDensity!= 0.0)
00219 //  if (material == nistwater || material->GetBaseMaterial() == nistwater)
00220   {
00221   const G4String& particleName = p->GetParticleName();
00222 
00223   if (ekin < highEnergyLimit)
00224   {
00225       //SI : XS must not be zero otherwise sampling of secondaries method ignored
00226       if (ekin < killBelowEnergy) return DBL_MAX;
00227       //      
00228       
00229         std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00230         pos = tableData.find(particleName);
00231         
00232         if (pos != tableData.end())
00233         {
00234           G4DNACrossSectionDataSet* table = pos->second;
00235           if (table != 0)
00236           {
00237             sigma = table->FindValue(ekin);
00238           }
00239         }
00240         else
00241         {
00242            G4Exception("G4DNAChampionElasticModel::ComputeCrossSectionPerVolume","em0002",
00243                       FatalException,"Model not applicable to particle type.");
00244         }
00245   }
00246 
00247   if (verboseLevel > 2)
00248   {
00249     G4cout << "__________________________________" << G4endl;
00250     G4cout << "°°° G4DNAChampionElasticModel - XS INFO START" << G4endl;
00251     G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
00252     G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
00253     G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
00254     //      G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
00255     G4cout << "°°° G4DNAChampionElasticModel - XS INFO END" << G4endl;
00256   } 
00257 
00258  } 
00259          
00260  return sigma*waterDensity;
00261 // return sigma*material->GetAtomicNumDensityVector()[1];
00262 }
00263 
00264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00265 
00266 void G4DNAChampionElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
00267                                               const G4MaterialCutsCouple* /*couple*/,
00268                                               const G4DynamicParticle* aDynamicElectron,
00269                                               G4double,
00270                                               G4double)
00271 {
00272 
00273   if (verboseLevel > 3)
00274     G4cout << "Calling SampleSecondaries() of G4DNAChampionElasticModel" << G4endl;
00275 
00276   G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
00277   
00278   if (electronEnergy0 < killBelowEnergy)
00279   {
00280     fParticleChangeForGamma->SetProposedKineticEnergy(0.);
00281     fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
00282     fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
00283     return ;
00284   }
00285 
00286   if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
00287   {  
00288   
00289     G4double cosTheta = RandomizeCosTheta(electronEnergy0);
00290   
00291     G4double phi = 2. * pi * G4UniformRand();
00292 
00293     G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
00294     G4ThreeVector xVers = zVers.orthogonal();
00295     G4ThreeVector yVers = zVers.cross(xVers);
00296 
00297     G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
00298     G4double yDir = xDir;
00299     xDir *= std::cos(phi);
00300     yDir *= std::sin(phi);
00301 
00302     G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
00303 
00304     fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
00305 
00306     fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
00307   }
00308 
00309 }
00310 
00311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00312 
00313 G4double G4DNAChampionElasticModel::Theta
00314   (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff)                                                    
00315 {
00316   G4double theta = 0.;
00317   G4double valueT1 = 0;
00318   G4double valueT2 = 0;
00319   G4double valueE21 = 0;
00320   G4double valueE22 = 0;
00321   G4double valueE12 = 0;
00322   G4double valueE11 = 0;
00323   G4double xs11 = 0;   
00324   G4double xs12 = 0; 
00325   G4double xs21 = 0; 
00326   G4double xs22 = 0; 
00327 
00328   if (particleDefinition == G4Electron::ElectronDefinition()) 
00329   {
00330     std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
00331     std::vector<double>::iterator t1 = t2-1;
00332  
00333     std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), integrDiff);
00334     std::vector<double>::iterator e11 = e12-1;
00335           
00336     std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), integrDiff);    
00337     std::vector<double>::iterator e21 = e22-1;
00338                 
00339     valueT1  =*t1;
00340     valueT2  =*t2;
00341     valueE21 =*e21;
00342     valueE22 =*e22;
00343     valueE12 =*e12;
00344     valueE11 =*e11;
00345 
00346     xs11 = eDiffCrossSectionData[valueT1][valueE11];
00347     xs12 = eDiffCrossSectionData[valueT1][valueE12];
00348     xs21 = eDiffCrossSectionData[valueT2][valueE21];
00349     xs22 = eDiffCrossSectionData[valueT2][valueE22];
00350   }
00351   
00352   if (xs11==0 && xs12==0 && xs21==0 && xs22==0) return (0.);   
00353 
00354   theta = QuadInterpolator   ( valueE11, valueE12, 
00355                                valueE21, valueE22, 
00356                                xs11, xs12, 
00357                                xs21, xs22, 
00358                                valueT1, valueT2, 
00359                                k, integrDiff );
00360                                
00361   return theta;
00362 }
00363 
00364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00365 
00366 G4double G4DNAChampionElasticModel::LinLogInterpolate(G4double e1, 
00367                                                         G4double e2, 
00368                                                         G4double e, 
00369                                                         G4double xs1, 
00370                                                         G4double xs2)
00371 {
00372   G4double d1 = std::log(xs1);
00373   G4double d2 = std::log(xs2);
00374   G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
00375   return value;
00376 }
00377 
00378 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00379 
00380 G4double G4DNAChampionElasticModel::LinLinInterpolate(G4double e1, 
00381                                                         G4double e2, 
00382                                                         G4double e, 
00383                                                         G4double xs1, 
00384                                                         G4double xs2)
00385 {
00386   G4double d1 = xs1;
00387   G4double d2 = xs2;
00388   G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
00389   return value;
00390 }
00391 
00392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00393 
00394 G4double G4DNAChampionElasticModel::LogLogInterpolate(G4double e1, 
00395                                                         G4double e2, 
00396                                                         G4double e, 
00397                                                         G4double xs1, 
00398                                                         G4double xs2)
00399 {
00400   G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
00401   G4double b = std::log10(xs2) - a*std::log10(e2);
00402   G4double sigma = a*std::log10(e) + b;
00403   G4double value = (std::pow(10.,sigma));
00404   return value;
00405 }
00406 
00407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00408 
00409 
00410 G4double G4DNAChampionElasticModel::QuadInterpolator(G4double e11, G4double e12, 
00411                                                        G4double e21, G4double e22, 
00412                                                        G4double xs11, G4double xs12, 
00413                                                        G4double xs21, G4double xs22, 
00414                                                        G4double t1, G4double t2, 
00415                                                        G4double t, G4double e)
00416 {
00417   // Log-Log
00418 /*
00419   G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
00420   G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
00421   G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
00422 
00423 
00424   // Lin-Log
00425   G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
00426   G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
00427   G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
00428 */
00429 
00430   // Lin-Lin
00431   G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
00432   G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
00433   G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
00434 
00435   return value;
00436 }
00437 
00438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00439 
00440 G4double G4DNAChampionElasticModel::RandomizeCosTheta(G4double k) 
00441 {
00442  
00443  G4double integrdiff=0;
00444  G4double uniformRand=G4UniformRand();
00445  integrdiff = uniformRand;
00446  
00447  G4double theta=0.;
00448  G4double cosTheta=0.;
00449  theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
00450 
00451  cosTheta= std::cos(theta*pi/180);
00452 
00453  return cosTheta; 
00454 }

Generated on Mon May 27 17:48:02 2013 for Geant4 by  doxygen 1.4.7